Investigating the Morphology of a Free-Falling Jet with an Accurate Finite Element and Level Set Modeling

Abstract

This study investigates the morphology of a free-falling liquid jet by using a computational approach with an experimental validation. Numerical simulations are developed by means of the Finite Element Method (FEM) for solving the viscous fluid flow and the level set method in order to track the interface between the fluid and air. Experiments are conducted in order to capture the shape of a free-falling jet of viscous fluid via circular orifice, where the shape is measured optically. The numerical results are found to be in agreement with the experimental data, demonstrating the validity of the proposed approach. Furthermore, we analyze the role of the surface tension by implementing linear as well as nonlinear surface energy models. All computational codes are developed with the aid of open-source packages from FEniCS and made publicly available. The combination of experimental and numerical techniques provides a comprehensive understanding of the morphology of free-falling jets and may be extended to multiphysics problems rather in a straightforward manner.

1. Introduction

Multi-phase flow is a wide field of research. In this paper, we are concerned with immiscible two-phase flows that involves moving interfaces. A wide range of technical applications exists such as fuel injection, fiber spinning, ink-jet printing, etc. [1], where two-phase flows are in use. The analytical analysis of this two-phase flow problem goes back to the 19th century. The key challenge is tracking the interface between the phases in order to obtain an accurate solution. The interface is modeled by a so-called surface tension that is measured in hydrostatics but used in fluid dynamics. Therefore, we need computational techniques for a better understanding of the role of surface tension in fluid flow. For modeling incompressible two-phase flows with surface tension, commonly used methods are the Volume Of Fluid (VOF) method [2,3,4], the Level Set (LS) method [5,6,7,8], the front tracking method [9], phase-field methods [10], and Smoothed Particle Hydrodynamics (SPH) [11,12,13,14]. For the scope of this paper, we focus on the LS method.

Originally, the LS method has been developed for the two-phase flow problem [15,16]. A level set function is a signed-distance function that represents the shortest distance to the interface. Comparing to VOF, the LS method has the advantage of accurate computation of interface normal and curvature. However, the biggest drawback of the LS method is that it causes mass loss, in other words, the balance of mass in the LS method is not conserved. As a remedy, Enright et al. [17] proposed a hybrid particle level set method in which they populate the interface with particles and use algorithms to attract them toward the interface and correct level set function as the particles cross the interface. Various researchers have analyzed that this mass loss may also be reduced by refining the mesh. For example, Herrmann [18] solves the level set equations on an auxiliary high-resolution grid and Gibou et al. [19] uses adaptive mesh refinement to refine the regions near the interface.

The Conservative Level Set (CLS) method is proposed by Olsson and Kreiss [20]. This model improves the mass conservation while maintaining the simplicity of the original method. We follow [21] and implement the CLS method herein. The main idea of the CLS method is to replace the signed distance function of the standard method with a regularized characteristic function. In this way, the level set function becomes a smeared Heaviside function, which is 0 within fluid and 1 in air; it varies smoothly from zero to one across the interface. This order parameter achieves the expected mass conservation. Zhao et al. [22] proposed an improved conservative level set method. In their work, the CLS method is used to capture interfaces, which has adequate mass conservation properties, and the standard LS method is employed to obtain accurate surface information.

Both the Finite Element Method (FEM) and the Finite Volume Method (FVM) are widely used numerical methods for solving fluid mechanics problems. FVM is adopted especially by commercial software solutions. This method is robust and stable, still developing in recent years [23,24,25,26]. However, it is challenging to increase the accuracy and use it for solving multiphysics problems. The Finite Element Method (FEM) is an accurate approximative technique with monotonous convergence properties. It is dominantly used in solid mechanics, and thus, an adequate candidate for fluid flow problems with a rather straightforward extension in fluid-structure interaction problems. In recent years, it has also been widely used in the research of fluid [27,28,29,30]. For example, FEM has delivered results in thermodynamics and electromagnetism by using standard form functions [31] accurately in coupled problems [32]. However, FEM has stability issues in solving fluid flow problems.

Surface tension models are crucial for multiphase fluid dynamics. Popinet (2018) reviewed key developments in Eulerian numerical models for simulating surface tension, emphasizing the stability of implicit time stepping and curvature estimation [33]. Hashemi et al. (2020) and Maarouf et al. (2022) proposed enriched finite element and characteristic-based finite element level set methods, respectively, to address strong pressure discontinuities and free interface problems, particularly focusing on the impact of surface tension [29,34]. Khalloufi et al. (2016) developed advanced numerical frameworks employing Arbitrary Lagrangian–Eulerian (ALE) methods and anisotropic mesh adaptation to enhance the computational efficiency and stability of free surface flows and two-phase fluids [35]. In terms of theoretical modeling of surface tension, Brackbill et al. (1992) introduced the Continuum Surface Force (CSF) model, which treats surface tension as a continuous three-dimensional effect across the interface, effectively handling complex interface topologies [36]. Recently, Kleinheins et al. (2023, 2024) further advanced surface tension models for multicomponent solutions, including the sigmoid model and a unified surface tension model, focusing on the effects of non-idealities and synergistic behaviors [37,38]. Additionally, Xie et al. (2022) extended the Gurtin–Murdoch surface elasticity theory and proposed a viscous surface hydrodynamics theory for modeling droplet wettability and moving contact line problems on soft substrates [39]. These studies collectively enhance our understanding and numerical simulation capabilities regarding the complex role of surface tension in multiphase flows and free surface dynamics. In this work, we use Xie’s model to investigate the effects of various factors within a complex surface tension coefficient model on the jet surface through numerical methods.

A stabilization method is often used for the FEM. For circumventing numerical problems, the inf-sup or so-called Ladyzhenskaya–Babuska–Brezzi condition needs to be fulfilled. One possibility is to use quadratic polynomials for the velocity and linear polynomials for the pressure in the element formulation. They are called P2/P1 or Taylor–Hood elements [40]. There is another way to overcome the numerical problems by using small elements that are smaller than the characteristic length scale, for example, Kolmogorov scale. This approach is called Direct Numerical Simulation (DNS) [41,42,43,44]. However, even though this method sounds simple, its implementation usually requires a lot of computing resources, such as access to a super-computer. There is another class of so-called projection methods originally introduced by Chorin [45] and its derivatives [46,47,48,49]. The staggered scheme bears difficulties in multiphysics extension. One common approach is to use stabilization terms [50,51] or multiscale methods [52]. Their extension to interface tracking is possible [53,54]. Without stabilization and by using P1/P1 elements, an accurate FEM is proposed by Abali [55], this approach is validated [56]. In recent years, there has been a widespread research interest in the integration of the level set method and the finite element method for studying multiphase flow [28,34,57,58]. This combined approach has gained significant attention and has shown promising applications, particularly in the field of additive manufacturing [59]. In this paper, we employ the latter approach for capturing two-phase flow interfaces. We study the surface shape of a free-falling fluid and compare it with experimental results to validate the implementation. Finally, in contrast to existing work, we consider a complex surface tension model proposed by Xie et al. [39], which comprises both linear and nonlinear components, and study the influence of each component on the results from a numerical perspective.

2. Numerical Method

2.1. Conservative Level Set Method

In the Level Set (LS) method, the level set function, $\varphi$, is defined as a signed distance function, ${\varphi }_d\left(\mathit{x}\right)$, as follows:

$${\varphi }_d\left(\mathit{x}\right)=\left\{\begin{array}{cc}\underset{{\mathit{x}}_{\mathrm{int}}\in \Gamma }{min}\left|\mathit{x}-{\mathit{x}}_{\mathrm{int}}\right|,& \mathit{x}\in {\Omega }_1,\\ -\underset{{\mathit{x}}_{\mathrm{int}}\in \Gamma }{min}\left|\mathit{x}-{\mathit{x}}_{\mathrm{int}}\right|,& \mathit{x}\in {\Omega }_2,\end{array}\right.$$

(1)

where ${\Omega }_1$ and ${\Omega }_2$ are subdomains that represent two different phases separated by an interface, $\Gamma ={\Omega }_1\bigcap {\Omega }_2$. The location of a fluid particle, $\mathit{x}$, is within the domain ${\Omega }_1\bigcup {\Omega }_2$; whereas ${\mathit{x}}_{\mathrm{int}}$ denotes the location of a particle on the interface, ${\mathit{x}}_{\mathrm{int}}\in \Gamma$. The meaning of the function, ${\varphi }_d\left(\mathit{x}\right)$, is the minimal distance from a location to the interface, $\Gamma$. In other words, ${\varphi }_d=0$ characterizes the interface, $\Gamma$.

According to [21], in the Conservative Level Set (CLS) method, the interface is tracked with a so-called regularized characteristic function:

$${\varphi }_c\left(\mathit{x}\right)={\displaystyle \frac{1}{1+exp\left({\varphi }_d\left(\mathit{x}\right)/\epsilon \right)}}.$$

(2)

This regularization creates a finite region, where ${\varphi }_c$ varies from 0 to 1 controlled by the thickness parameter, $\epsilon$, which corresponds to the half of the interface thickness. The thickness parameter, $\epsilon$, is determined by the element size in the mesh and a finer mesh implies a smaller value of $\epsilon$. In this way, at the interface, ${\varphi }_d=0$ results in ${\varphi }_c=0.5$ whereas ${\varphi }_c$ is a diffuse distribution from 0 to 1 along the distance of $2\epsilon$. After a simple initialization, the result of ${\varphi }_c$ is shown in Figure 1. In order to prevent numerical issues, the value of the thickness parameter should be slightly larger than the element size. In simpler words, sufficiently many elements need to be used within the transition region of length $2\epsilon$. As the mesh becomes finer, the thickness parameter may be decreased proportionally to enhance the accuracy of the interface. In simpler words, the accuracy is mesh-dependent and we aim for a monotonous convergence in ${\varphi }_c$ in order to ensure the reliability of the implementation.

2.2. Level Set Advection and Reinitialization

The interface represents a fictitious surface between the two phases. Since the interface is attached to the massive particles of one of the phases, its motion is governed by the motion of the massive particles. Hence, in order to achieve the evolution of the interface, the advection equation is solved–stemming from the balance of mass—as follows:

$${\displaystyle \frac{\partial {\varphi }_c}{\partial t}}+\mathbf{\nabla }·\left({\varphi }_c\mathit{v}\right)=0.$$

(3)

This equation indicates that the function ${\varphi }_c(\mathit{x},t)$ is advected under the velocity field, $\mathit{v}$, which is indeed the velocity of the massive particles of the fluid. Therefore, $\mathit{v}$ is solved from the balance of momentum with (linear) constitutive equations resulting in Navier–Stokes equations.

Reinitialization is a technique used in the level set method in order to ensure that the characteristic function, ${\varphi }_c$, remains accurate over time. This characteristic function represents the interface between the two fluids or materials. In the level set method, we only define the signed distance function, ${\varphi }_d$, and the characteristic function, ${\varphi }_c$, at the initial time step to define the initial interface. Transiently, we solve the advection Equation (3), where ${\varphi }_c$ may become distorted due to numerical errors or inaccuracies in the numerical solution. Reinitialization is a technique used to rectify these errors and reinitialize ${\varphi }_c$ at every time step. The basic idea behind reinitialization is solving some advection-like equation, for acquiring ${\varphi }_r(\mathit{x},\tau )$, depending on a sub-time variable, $\tau$. A governing equation for ${\varphi }_r(\mathit{x},\tau )$ is solved several times until convergence is achieved within the time step, $t_i$. By following [21], we choose the following reinitialization equation:

$${\displaystyle \frac{\partial {\varphi }_r}{\partial \tau }}+\mathbf{\nabla }·\left[{\varphi }_r\left(1-{\varphi }_r\right){\mathbf{n}}_{\Gamma }\right]=\epsilon \mathbf{\nabla }·\left[{\mathbf{n}}_{\Gamma }\left(\mathbf{\nabla }{\varphi }_r·{\mathbf{n}}_{\Gamma }\right)\right],$$

(4)

where ${\mathbf{n}}_{\Gamma }\left(\mathit{x},t_i\right)$ is the surface normal on the interface, $\epsilon$ is the thickness parameter in (2). The second part on the left-hand side of (4), ${\varphi }_r\left(1-{\varphi }_r\right){\mathbf{n}}_{\Gamma }$, corresponds to a so-called compression flux. This artificial compression flux acts in regions where $0<{\varphi }_r<1$ and along the surface normal at the interface. On the right-hand side of Equation (4), a small amount of artificial viscosity term is introduced to avoid numerical issues caused by discontinuities at the interface. The surface normal at the interface, ${\mathbf{n}}_{\Gamma }\left(\mathit{x},t_i\right)$, is given by

$${\mathbf{n}}_{\Gamma }={\displaystyle \frac{\mathbf{\nabla }{\varphi }_c(\mathit{x},t)}{|\mathbf{\nabla }{\varphi }_c(\mathit{x},t)|}}.$$

(5)

For the initial condition, $\tau =0$, we choose

$${\varphi }_r(\mathit{x},0)={\varphi }_c\left(\mathit{x},t_i\right).$$

(6)

2.3. Governing Equations

For an incompressible flow, the mass density is assumed to be constant in time, then continuity equation reads

$$\mathbf{\nabla }·\mathit{v}=0.$$

(7)

Balance of momentum,

$$\rho {\displaystyle \frac{\partial \mathit{v}}{\partial t}}+\rho (\mathit{v}·\mathbf{\nabla })\mathit{v}=\mathbf{\nabla }·\mathit{\sigma }+\mathit{f},$$

(8)

incorporates the Cauchy stress tensor, $\mathit{\sigma }$, to be defined as well as volumetric forces, $\mathit{f}$, to be given. For the stress tensor, we use the well-known linear material model,

$$\mathit{\sigma }=-p\mathbf{I}+\mu \left(\mathbf{\nabla }\mathit{v}+{(\mathbf{\nabla }\mathit{v})}^{⊺}\right),$$

(9)

where $\mu$ is viscosity in $\mathrm{Pa}·\mathrm{s}$. As external forces, we employ gravitational body forces, ${\mathit{f}}_g=\rho g{\mathit{e}}_g$ as well as the forces due to surface tension ${\mathit{f}}_{\mathrm{s}}$.

Surface tension is modeled at the interface by forces per area as given by ${\varphi }_c$. In two-phase flow, free surfaces has garnered significant scholarly interest about modeling the surface tension [29,34,35,60,61]. The usual expression from the literature [16] reads

$${\mathit{f}}_{\mathrm{s}}={\mathit{\sigma }}_{\mathrm{s}}\kappa {\mathbf{n}}_{\Gamma }\delta ,$$

(10)

and we refer to [36] for a derivation. The scalar ${\mathit{\sigma }}_{\mathrm{s}}$ is the surface tension coefficient in N/mm, ${\mathbf{n}}_{\Gamma }$ is the unit surface normal on the interface given by (5), and $\delta$ represents an approximation of the Dirac distribution, which is situated at the interface denoted by $\Gamma$. The curvature of the interface, $\kappa$, is geometrically given by

$$\kappa =-\mathbf{\nabla }·{\mathbf{n}}_{\Gamma }=-\mathbf{\nabla }·{\displaystyle \frac{\mathbf{\nabla }{\varphi }_c}{|\mathbf{\nabla }{\varphi }_c|}}.$$

(11)

By means of

$$\delta =\left|\mathbf{\nabla }{\varphi }_c\right|,$$

(12)

we obtain from Equation (10)

$${\mathit{f}}_{\mathrm{s}}={\mathit{\sigma }}_{\mathrm{s}}\kappa \mathbf{\nabla }{\varphi }_c.$$

(13)

Alternatively, by expressing the surface tension force as the divergence of a surface stress tensor in an analogous manner to the bulk properties [62], as follows:

$$\begin{array}{c}{\mathbf{T}}_{\mathrm{s}}={\mathit{\sigma }}_{\mathrm{s}}\left(\mathbf{I}-{\mathbf{n}}_{\Gamma }\otimes {\mathbf{n}}_{\Gamma }\right)\delta ,\\ {\mathit{f}}_{\mathrm{s}}=\mathbf{\nabla }·{\mathbf{T}}_{\mathrm{s}},\end{array}$$

(14)

where $\delta$ is computed by Equation (12). An advantage with Equation (14) over Equation (13) is the absence of curvature computation.

Furthermore, Xie et al. [39] proposed the adhesive surface hydrodynamics theory that takes into account the influence of velocity gradients on interface in a surface tension model. In other words, viscous friction force during the dynamic flow process is taken into account, as shown in the following equation:

$$\begin{array}{c}{\mathbf{T}}_{\mathrm{hydro}}=\left({\mathit{\sigma }}_{\mathrm{s}}{\mathbf{I}}_{\mathrm{s}}+{\mu }_{\mathrm{s}}\mathbf{P}\mathrm{Sym}\left(\mathbf{\nabla }\otimes \mathit{v}\right)\mathbf{P}\right)\delta ,\\ {\mathit{f}}_{\mathrm{hydro}}=\mathbf{\nabla }·{\mathbf{T}}_{\mathrm{hydro}}.\end{array}$$

(15)

We have extended the tensor expression of this model to the form suitable for the level set method. Here ${\sigma }_{\mathrm{s}}$ represents the same surface tension coefficient as in Equation (14). The second part considers the contribution of viscous friction with a surface viscosity, ${\mu }_{\mathrm{s}}$, at the interface between the two media such as viscous fluid and air. In this case, we consider the velocity of viscous fluid particles at the interface, in other words, air stands still. In general, this term adds an additional term analogous from the Euler to Newton fluids in the bulk. There are no well established experimental procedures for determining this parameter. In this work, we investigate the influence of this parameter on the dynamic flow process qualitatively by assigning ${\mu }_{\mathrm{s}}$ different numerical values. The terms, ${\mathbf{I}}_{\mathrm{s}}$ and $\mathbf{P}$, are shown as

$$\begin{array}{c}{\mathbf{I}}_{\mathrm{s}}=\mathbf{P}\mathbf{I}=\mathbf{P},\\ \mathbf{P}=\mathbf{I}-{\mathbf{n}}_{\Gamma }\otimes {\mathbf{n}}_{\Gamma }.\end{array}$$

(16)

Consequently, a diffuse interface model is used in the surface tension force expression, which means that the surface tension is transformed into a volume force in the interface area and this area is composed of a few layers of elements related to the thickness parameter, $\epsilon$, used in Equation (2). Hereby, we use the diffuse interface model to make these quantities vary smoothly over the interface, by adopting the rule of mixtures,

$$\begin{array}{cc}\hfill \rho & ={\rho }_{\mathrm{fluid}}{\varphi }_c+{\rho }_{\mathrm{air}}\left(1-{\varphi }_c\right),\hfill \\ \hfill \mu & ={\mu }_{\mathrm{fluid}}{\varphi }_c+{\mu }_{\mathrm{air}}\left(1-{\varphi }_c\right).\hfill \end{array}$$

(17)

We use the following conversion with characteristic length and time-scales in order to acquire dimensionless quantities,

$$\overline{\mathit{x}}={\displaystyle \frac{\mathit{x}}{x_0}},\overline{\mathit{v}}={\displaystyle \frac{\mathit{v}}{v_0}},\overline{t}={\displaystyle \frac{t}{x_0/v_0}},\overline{\rho }={\displaystyle \frac{\rho }{{\rho }_0}},\overline{\mu }={\displaystyle \frac{\mu }{{\mu }_0}}.$$

(18)

Hence, we obtain the dimensionless form of Equations (7) and (8), as follows:

$$\mathbf{\nabla }·\left(\overline{\rho }\mathit{v}\right)=0.$$

(19)

$$\begin{array}{cc}\hfill \overline{\rho }{\displaystyle \frac{\partial \overline{\mathit{v}}}{\partial \overline{t}}}+\overline{\rho }(\overline{\mathit{v}}·\mathbf{\nabla })\overline{\mathit{v}}=& -\mathbf{\nabla }\overline{p}+{\displaystyle \frac{\overline{\mu }}{Re}}\mathbf{\nabla }·\left(\mathbf{\nabla }\overline{\mathit{v}}+{(\mathbf{\nabla }\overline{\mathit{v}})}^{⊺}\right)+{\displaystyle \frac{\overline{\rho }}{F{r}^2}}{\mathit{e}}_g+{\displaystyle \frac{1}{We}}\overline{{\mathit{f}}_{\mathrm{s}}},\hfill \end{array}$$

(20)

where ${\mathit{e}}_g$ is the unit vector representing the direction of gravity, and $\overline{{\mathit{f}}_{\mathrm{s}}}$ is the nondimensionalized surface tension. Reynolds number, $Re$, Weber number, $We$, and Froude number, $Fr$, read

$$Re={\displaystyle \frac{{\rho }_0{v}_0{x}_0}{{\mu }_0}},Fr={\displaystyle \frac{v_0}{\sqrt{x_{0}g}}},We={\displaystyle \frac{{\rho }_0{v}_0^2{x}_0}{{\mathit{\sigma }}_{\mathrm{s}}}}.$$

(21)

For time discretization, we use the Finite Difference Method (FDM) also called the Euler backward method:

$${\displaystyle \frac{\partial \overline{\mathit{v}}}{\partial t}}={\displaystyle \frac{\overline{\mathit{v}}-{\overline{\mathit{v}}}_0}{\Delta \overline{t}}},$$

(22)

where $\Delta \overline{t}$ is time step and $\overline{\mathit{v}}$ is an unknown at the current time step with the already computed (known) ${\overline{\mathit{v}}}_0$ in the last time step.

By following a standard variational formulation, for the unknowns velocity and pressure with their variations, $\delta \overline{\mathit{v}}$ and $\delta \overline{p}$, and performing partial integration on the surface tension part, we obtain the weak forms:

$$F_1={\int }_{\Omega }\overline{\rho }\mathbf{\nabla }·\left(\overline{\mathit{v}}\right)\delta \overline{p}\mathrm{d}x.$$

(23)

$$\begin{array}{cc}\hfill F_2=& {\displaystyle \frac{1}{\Delta \overline{t}}}{\int }_{\Omega }\overline{\rho }\left(\overline{\mathit{v}}-{\overline{\mathit{v}}}_0\right)·\delta \overline{\mathit{v}}\mathrm{d}x+{\int }_{\Omega }\overline{\rho }(\overline{\mathit{v}}·\mathbf{\nabla })\overline{\mathit{v}}·\delta \overline{\mathit{v}}\mathrm{d}\mathit{x}\hfill \\ & +{\int }_{\Omega }\mathbf{\nabla }\overline{p}·\delta \overline{\mathit{v}}\mathrm{d}x+{\displaystyle \frac{1}{Re}}{\int }_{\Omega }\overline{\mu }\mathbf{\nabla }\delta \overline{\mathit{v}}·\left(\mathbf{\nabla }\overline{\mathit{v}}+{(\mathbf{\nabla }\overline{\mathit{v}})}^{\mathrm{T}}\right)\mathrm{d}x\hfill \\ & -{\int }_{\Omega }\delta \overline{\mathit{v}}·{\displaystyle \frac{\overline{\rho }}{F{r}^2}}{\mathit{e}}_g\mathrm{d}x+{\int }_{\Omega }{\displaystyle \frac{1}{We}}{\mathbf{T}}_{\mathrm{s}}·\mathbf{\nabla }\delta \overline{\mathit{v}}\mathrm{d}x.\hfill \end{array}$$

(24)

We emphasize the choice of integration by parts only on the terms where needed. This choice eliminates some numerical problems despite the fact that it is not a common practice. The use of P2/P1 elements for velocity and pressure is adequate but computationally costly. Especially in the case of using fine meshes because of the interface, we aim to use a P1/P1 scheme. Therefore, we follow [55] and solve an additional weak form stemming from the balance of angular momentum and test this one by $\mathbf{\nabla }\delta \overline{p}$ in order to prevent numerical instabilities. In short, this term acts as a stabilization term without adding an artificial viscosity or limiting procedure. Hence, we add

$$\begin{array}{cc}\hfill F_3& ={\int }_{\Omega }(\left(\overline{\mathit{v}}-{\overline{\mathit{v}}}_0\right)+\Delta \overline{t}(\overline{\mathit{v}}·\mathbf{\nabla })\overline{\mathit{v}}+{\displaystyle \frac{\Delta \overline{t}}{\overline{\rho }}}\mathbf{\nabla }\overline{p}-{\displaystyle \frac{\Delta \overline{t}}{F{r}^2}}{\mathit{e}}_g\hfill \\ & -{\displaystyle \frac{\Delta \overline{t}}{\overline{\rho }We}}{\mathit{f}}_{\mathrm{s}}-{\displaystyle \frac{\Delta \overline{t}\overline{\mu }}{\overline{\rho }Re}}\mathbf{\nabla }·\left(\mathbf{\nabla }\overline{\mathit{v}}+{(\mathbf{\nabla }\overline{\mathit{v}})}^{\mathrm{T}}\right))·\mathbf{\nabla }\delta \overline{p}\mathrm{d}x,\hfill \end{array}$$

(25)

and solve

$${\left.F\right|}_{\Omega }=F_1+F_2+F_3.$$

(26)

We use standard Lagrange elements with a linear interpolation function for velocity and pressure. The standard FEM formulation with piecewise continuity ensures a monotonous convergence and the unknown fields are solved at once providing accurate solution in coupled terms. Shape functions are piecewise (elementwise) continuous over the whole domain with a compact support on the element such that nodal values have influence locally. This representation is adequate for approximation in Hilbertian Sobolev space ${\mathcal{H}}^1$ [63]. Elementwise interpolation is of order q and span ${\mathcal{P}}_q$ on a triangulation, $\mathcal{T}$, consisting of elements, ${\Omega }^{\mathrm{E}}$. Linear Lagrange elements are used for pressure and velocity, which means we use interpolation polynomial of first order as shape function for each element in both scalar space and vector space,

$$\begin{array}{c}\mathcal{V}=\left\{\{\overline{p},\overline{\mathit{v}}\}\in {\mathcal{H}}^1(\Omega ):\{\overline{p},\overline{\mathit{v}}\}{|}_{{\Omega }^{\mathrm{E}}}\in {\mathcal{P}}_1\left({\Omega }^{\mathrm{E}}\right)\forall {\Omega }^{\mathrm{E}}\in \mathcal{T}\right\}.\end{array}$$

(27)

For the advection and reinitialization, the weak forms are

$$F_{{\varphi }_c}={\int }_{\Omega }\delta {\varphi }_c{\displaystyle \frac{{\varphi }_c-{\varphi }_{c0}}{\Delta \overline{t}}}\mathrm{d}x-{\int }_{\Omega }\mathbf{\nabla }\delta {\varphi }_c·\left({\varphi }_c\overline{\mathit{v}}\right)\mathrm{d}x$$

(28)

$$\begin{array}{cc}\hfill F_{{\varphi }_r}=& {\int }_{\Omega }\delta {\varphi }_c{\displaystyle \frac{{\varphi }_r-{\varphi }_{r0}}{\Delta \overline{\tau }}}\mathrm{d}x-{\int }_{\Omega }({\varphi }_r\left(1-{\varphi }_r\right)\left(\mathbf{\nabla }\delta {\varphi }_c·{\mathbf{n}}_{\Gamma }\right)\hfill \\ & +\epsilon \left({\mathbf{n}}_{\Gamma }·\mathbf{\nabla }{\varphi }_r\right)\left({\mathbf{n}}_{\Gamma }·\mathbf{\nabla }\delta {\varphi }_c\right))\mathrm{d}x.\hfill \end{array}$$

(29)

Analogously, for solving ${\varphi }_c$ and ${\varphi }_r$, we use linear Lagrange elements on the same mesh, spanning (identical) scalar spaces

$$\begin{array}{c}\mathcal{S}=\left\{\left\{{\varphi }_{r/c}\right\}\in {\mathcal{H}}^1(\Omega ):\left\{{\varphi }_{c/r}\right\}{|}_{{\Omega }^{\mathrm{E}}}\in {\mathcal{P}}_1\left({\Omega }^{\mathrm{E}}\right)\forall {\Omega }^{\mathrm{E}}\in \mathcal{T}\right\}.\end{array}$$

(30)

2.4. Algorithm Scheme

The algorithm scheme is illustrated in Figure 2. The process involves two iterations, both pertaining to the time step. In the first iteration, we solve the balance equations, and in the second iteration, we solve the advection equation based on the velocity field results. Subsequently, in the third step, which is a sub-iteration at each time step, we carry out the reinitialization process.

3. Free Falling Jet Problem

The behavior of a fluid as it pours from an outlet into the air, forming a free-falling stable jet that accelerates, stretches, and narrows under the influence of gravity, is of significant interest in the fields of fluid mechanics and engineering practice [64,65,66]. The study of this jet flow behavior has found diverse applications, such as in the production of small fluid particles using the sol-gel process, fabrication of polymer fibers using spinning processes, and development of biomedical devices.

In this paper, we establish a 2D finite element model based on the free-falling jet subjected to gravity and validate the accuracy of the proposed model for interface simulation by comparing it with experimental results. In Figure 3 boundary regions are separated by marking them in different colors. The fluid enters across the blue boundary, passes through a $1\times 1$ pipe (a dimensionless geometric parameter), and falls freely under the influence of gravity (with an initial velocity given at the blue boundary, inlet). The no-slip boundary condition is set on yellow boundaries, while the pressure is set to 0 on the green boundary, outlet. Since ${\varphi }_c=1$ denotes fluid, we emphasize that for simulating the continuous inflow of fluid, we set the condition of ${\varphi }_c=1$ on the blue boundary for all times. The boundary conditions are summarized as follows:

$$\begin{array}{ccc}\hfill \mathit{u}=& \left(0,-v_0\right),\hfill & \hfill x\in {\Gamma }_1,\\ \hfill {\varphi }_c=& 1,\hfill & \hfill x\in {\Gamma }_1,\\ \hfill p=& 0,\hfill & \hfill x\in {\Gamma }_2,\\ \hfill \mathit{u}=& (0,0),\hfill & \hfill x\in {\Gamma }_3.\end{array}$$

(31)

We model the ejection by setting velocity boundary condition, $\mathit{u}$, and level set function value boundary condition, ${\varphi }_c$, on ${\Gamma }_1$ (the blue boundary) in Figure 3. A constant velocity $v_0$ directed downwards along the vertical axis is established on this boundary, and the level set function value is consistently set to 1, facilitating the simulation of continuous flow at this boundary. Once this flow stabilizes, a steady jet is successfully simulated.

4. Experiment

In order to investigate the same problem, we conduct a simple experiment and measure the shape of a liquid jet falling freely under the gravitational forces. In Figure 4 we set up a vertical apparatus consisting of a reservoir, a valve to regulate flow rate, a pump to generate the flow, and a nozzle through which the liquid is delivered. A flow regulating valve allows for the adjustment of flow rate, enabling the formation of a stable jet by selecting an appropriate flow rate manually. The current velocity is determined by considering both the flow rate and the diameter of the pip. The velocities of different fluids used in the experiment are listed in

Table 1. Parameters setting.

Fluid	Viscosity (Pa·s)	Surface Tension (μN/mm)	Density (g/mm3)	Velocity (mm/s)	Re
Peanut oil	0.061	34	0.91 $\times {10}^{-3}$	500	59.67
Castor oil	0.57	35.9	0.956 $\times {10}^{-3}$	400	5.38
Glycerol (18 wt%)	1.3	72.7	1.26$\times {10}^{-3}$	200	1.55

. A container is placed beneath the nozzle to collect the liquid that fell from the jet, and a grid background is positioned behind the jet to measure its diameter by a simple image analysis. The grid is placed close enough to the jet to obtain accurate measurements without interfering with its trajectory.

We adjusted the flow rate until the jet was stable and straight, falling vertically downward without deviation. Using a video camera, we captured the motion of the falling jet to record its shape. The specific model of the camera system was the Hitachi KP-FBR30PCL. (Hitachi KP-FBR30PCL 1/3” CCD, Remote Head Camera, NTSC, Color Raw, 659H × 494V, Camera Link, Lens type: C-Mount, Output: Dual 26-Pin SDR, 60fps. This camera system satisfied requirements for high-speed shutter photography and output of raw files, which are advantageous for post-processing needs. The reservoir had to be positioned above the outlet, with the jet primarily driven by gravity. The valve was used to restrict the flow so that the gravity-driven jet did flow steadily at a fixed outlet rate for a period of time, allowing for multiple photographs to be taken. For ensuring accurate measurements, we used a camera with a long focal length and a small aperture. The long focal length minimized the effects of perspective distortion and lens aberrations, while the small aperture increased the depth of field and allowed both the foreground liquid jet and the background grid to be in focus. Additionally, the jet and the background grid had to be as close as possible to minimize the difference between the statistical results and the actual dimensions.

In this experiment, we investigated three types of fluids: peanut oil, castor oil, and glycerin. The specific material properties of these substances are found in the literature and presented in Table 1.

In order to process the results shown in Figure 5, images of the stable liquid jet were captured. Several points were then selected along the flow direction, and the diameter of the liquid jet was calculated and recorded at each position using the scale provided on the left side of Figure 5 and the background grid. For each stable jet, 4–5 sets of data were taken. The remaining data were averaged to obtain the surface profile result for the current jet.

5. Numerical Results

The numerical simulation is carried out using the open-source packages developed by the FEniCS platform. (The FEniCS computing platform, https://fenicsproject.org/) In the simulation, we utilize the dimensionless parameters after Equation (18), where the reference length $x_0$ is set to 8 mm (the fluid outlet width), and the reference velocity $v_0$ and other reference material parameters are compiled in Table 1 based on the properties of each material. All codes are made public on Github. (Github: https://github.com/liuyiming0507/FreeFallingJet).

Initially we conduct a convergence test on the results of the surface tension model. Figure 6b presents the normal stress, ${\mathrm{T}}_{\mathrm{yy}}^s$, as it results from surface tension tensor, ${\mathbf{T}}_{\mathrm{s}}$, computed by different Degrees Of Freedom (DOFs). In FEniCS, the Degrees Of Freedom (DOFs) are determined based on the mixed variable space that contains all the unknowns in the fluid governing equations, specifically the velocity and pressure. The total number of DOFs corresponds to the sum of the independent variables required to describe the solution, accounting for all components of velocity and the pressure across the entire computational domain. We select the extreme value of ${\mathrm{T}}_{\mathrm{yy}}^s$ at the boundary interface where $y=9.5$ when the initial part of the jet stream is stable at $\overline{t}=5$, as seen in Figure 6a. The component ${\mathrm{T}}_{\mathrm{yy}}^s$ of the surface tension tensor has been selected for convergence results. In the level set method, the transition region, representing the interface, becomes narrower as the mesh is refined, leading to a more accurate representation of the interface. Since the surface tension is influenced by this transition region (based on Equation (14), Equation (12), Equation (5) and Equation (2)), choosing a representative component like ${\mathrm{T}}_{\mathrm{yy}}^s$ is crucial for assessing convergence. Furthermore, ${\mathrm{T}}_{\mathrm{yy}}^s$ is selected specifically because normal stress generally has larger numerical values in $yy$ direction.

The convergence of ${\mathrm{T}}_{\mathrm{yy}}^s$ is relatively slow and shows that it has strong mesh dependence. We used a uniform triangular mesh, shown in Figure 7, to eliminate the influence of irregular mesh on the results.This phenomenon relies on the surface tension model Equation (14) in this study, based on the level set method. In Equation (2), the choice of thickness parameter, $\epsilon$, depends on the size of the element. A finer mesh leads to a smaller $\epsilon$ value, resulting in a thinner transition interface region. Obviously, with a thinner interface, the actual “interface” is more accurately modeled, and thus, the area, where ${\mathrm{T}}_{\mathrm{yy}}^s$ acts, becomes smaller, see Figure 8. The results of ${\varphi }_c=0.5$ are maintaining a consistent interface contour as in Figure 9 at different resolutions. As contours at resolutions 40 and 50 are indistinguishable, the system has converged.

The results of numerical simulations and experimental comparisons are shown in Figure 10. We focus on the shape of free-falling jet and investigate the effects of Reynolds number on the jet shape. We transform the dimensionless outcomes of the numerical simulation into physical values for comparing the interface between the liquid jet and air (the contour of ${\varphi }_c=0.5$) with experimental data. In Figure 10, the solid line represents the numerical results, while the experimental results are depicted as discrete data points. The comparison between experimental and numerical results reveal that, at high Reynolds numbers, the inertial forces dominate the flow behavior, and the numerical simulations match the experimental jet shape well. However, for small Reynolds numbers, the numerical simulations fail to capture the initial contraction of the jet with the same accuracy. A possible explanation for this discrepancy is that the employed surface tension model is based on an interface area approximation derived from the level set method, Equation (13) and Equation (11), which introduces errors. At low Reynolds numbers, the dominance of the inertial term diminishes, leading to the emergence of these errors in the surface tension model, especially in areas of large curvature. Additionally, for low Reynolds numbers, viscous force have a greater influence on the fluid compared to inertial force. In such cases, using a linear viscosity model may not accurately capture the behavior. In this study, we further explored the effects of a complex surface tension models, and the impact of other viscosity models is not further considered. Despite this, for the jet tail, the maximum diameter among these three fluids is around 4∼5 mm, while the minimum is around 2 mm, showing good agreement between the experimental and numerical results. Based on the experimental data collection points in Figure 5, we calculated the mean relative error (MRE) of the numerical results at these points compared to the experimental results for each case. The results are as follows: the MRE for Peanut oil is $3.59\%$, for Castor oil is $5.98\%$, and for Glycerol is $9.72\%$. Overall, the mean relative errors are acceptable.

The analyzed 2D model is a reduction in space and needs to be validated in a 3D simulation. We stress that the implementation is working in 2D as well as 3D, indeed, the degrees of freedom and thus the computation time increases significantly. We have used a computing node using Intel Xeon E7-4850, in total 64 cores each with the 40 MB cache, equipped with 256 GB memory in total, running Linux Kernel 5 Ubuntu 20.04. The parallelization is done by mpirun and the data is merged in Paraview. The conducted 3D simulation is dempnstrated in Figure 11, revealing that the finite element framework based on FEniCS established in this study can effectively simulate 3D cases in case of having enough computational power. In Figure 12, we used Glycerol material data as an example to compare the results of 3D simulation with those of 2D simulation. At the regions of interest, namely the initial contraction of the jet and the droplet at the jet tail, it is observed that the 2D and 3D results exhibit good agreement, indicating that simplifying the problem to 2D has minimal impact on the cases investigated in this research. Considering time efficiency and improved productivity, the investigation of 2D cases suffices for the scope of this study.

Next, we compare the impact of different Weber numbers on the results. The Weber number is a dimensionless number in fluid mechanics that compares the inertia to the surface tension. A smaller Weber number indicates that the surface tension has a larger influence on the fluid behavior compared to the inertia. We keep the Reynolds number and the Froude number constant, selecting $We=0.5,1,5,10$. The comparison results of ${\varphi }_c$ are depicted in Figure 13. As seen in Equation (21), a smaller Weber number implies that the influence of surface tension on the fluid, compared to the inertia, is greater. From these results, it is evident that at lower Weber numbers the fluid flow is slower and the droplets at the jet tail are larger. Conversely, a higher Weber number indicates a larger influence of the inertia over the surface tension. As the Weber number increases, the fluid flow becomes faster and the droplets at the jet tail become smaller.

Furthermore, we study the influence of different components in Equation (15). In Figure 14a, by keeping ${\mathit{\sigma }}_{\mathrm{s}}$ constant and comparing the results for different values of ${\mu }_{\mathrm{s}}$, we observe an analogous trend as seen in the results for different Weber numbers. In this surface hydrodynamics model, the variation in ${\mu }_{\mathrm{s}}$ has a small effect on the droplet size at the jet tail and also the curvature at the droplet head. Overall, the shape changes are minimal. However, as shown in Figure 14b, the choice of ${\mu }_{\mathrm{s}}$ significantly affects the contour of the initial contraction section of the jet. In Figure 15, we investigate the influence of different ${\sigma }_{\mathrm{s}}$ values on the results. In Figure 15a, comparing the results with those in Figure 14a, it is evident that there is a significant variation in the size of the droplets at the jet tail. Under lower ${\sigma }_{\mathrm{s}}$, the droplets appear narrow and small, whereas as the ${\sigma }_{\mathrm{s}}$ increases, the droplets become more rounded, indicating a pronounced effect of the surface tension coefficient on the droplets’ size at the jet tail. In contrast, Figure 15b demonstrates that different ${\sigma }_{\mathrm{s}}$ have minimal impact on the contour of the initial contraction section of the jet. We have discussed the impact of the surface tension term and the viscous friction term in Equation (15) based on the results. In summary, the surface tension term affects the size and shape of the droplets at the jet tail, while the viscous friction term influences the contour shape of the initial contraction section of the jet. It is important to highlight that the viscous friction term is expressed through the velocity gradient, which is accompanied by a coefficient ${\mu }_{\mathrm{s}}$. Although the present study did not specifically determine this coefficient, it is noteworthy that variations in this coefficient have been reported to significantly influence the initial contraction of the jet in numerical simulations. Thus, the role of this coefficient warrants further attention in future investigations to better understand its impact on jet dynamics.

6. Conclusions

In this paper, we investigate the morphology of a free-falling liquid jet by comparing experimental outcomes with numerical simulations and discussing the impact of Reynolds number on jet shape. At high Reynolds numbers, inertial forces prevail in the flow behavior, resulting in a strong correspondence between numerical simulations and experimental jet shapes. On the other hand, at low Reynolds numbers, numerical simulations exhibit reduced accuracy in capturing the initial contraction of the jet. Nonetheless, at the jet tail, a notable agreement is observed between the experimental and numerical findings. In the numerical simulation, we employ an accurate finite element method, which effectively reduces computation time. The concurrence with experimental results demonstrates that our numerical model is both accurate and viable. Next, we investigate the effect of the Weber number. The Weber number provides an important method for characterizing the effects of surface tension and inertia on the behavior of the fluid in our model.

Finally, we investigate the impact of the surface tension term and the viscous friction term on the results within the applied surface hydrodynamics theory. From the conclusions, it is observed that these two terms exhibit a certain level of independence in their influence on the numerical computation results. From Figure 14b, we observe the influence of the viscous friction component in the new surface tension model on the initial contraction of the jet. In addition, Figure 16 also demonstrates the distribution of the x-direction velocity component (horizontal velocity component) in the initial contraction part. This solution indicates that the velocity gradient component in the selected adhesive surface hydrodynamics model is worthy of consideration. The ability to control and understand these effects is crucial for predicting and manipulating the behavior of fluid flows in a variety of applications such as solder ball jetting, 3D-printing, and ink-jet printing [67,68,69].