Numerical Study of an Indicator Function for Front-Tracking Methods

Department of Mathematics, Kangwon National University, Chuncheon 24341, Gangwon-do, Republic of Korea Department of Mathematics, Korea University, Seoul 02841, Republic of Korea School of Computer Science and Engineering, Sun Yat-Sen University, Guangzhou 510275, China Department of Mathematics, Zhengzhou University of Aeronautics, Zhengzhou 450046, China Department of Mathematics, Harbin Institute of Technology at Weihai, Weihai 264209, China


Introduction
Dirac delta function plays an important role in various scienti c and industrial problems. One of them is on the numerical simulations of uid ows including multiphase ows and uid solid interaction (FSI) problems. e immersed boundary methods (IBM) are on the front burner for the numerical simulations of these problems and one of the key points lay on the methods to accurately separate the moving multiphase or solid boundaries. e phase eld method (PHM) [1][2][3][4], level set method (LSM) [5][6][7], and volume of uid methods (VOF) [8][9][10], express the di erent phases using particular scalar functions. ese functions are designed to be advected by uid ows and one phase changes to another one continuously to avoid numerical stability problem. It should be noted that the interphase is di used and implicitly captured in the three methods. However, the immersed boundary method is di erent for its explicit boundary expression. Unlike PHM, LSM, and VOF, distinguish of di erent phase can be denoted directly by scalar phase functions; in IBM, a so-called indicator function needs to be constructed for these phase distinctions indirectly. In this process, if a Dirac delta function is used for the construction of the indicator function, the di used interaction forms as PHM, LSM, and VOF. From this point of view, these four methods are all belonged to di used models in simulations of multiphase ows. In IBM, the uid interface is explicitly represented by several Lagrangian points and the uid ow is solved on a stationary Eulerian mesh. Furthermore, when solving Richards' equation [11,12] for layered soils, the layers of soils can be expressed using an indicator function.
Multiphase uid ows are important in various scienti c and industrial problems.
In the front-tracking method [21] such as IBM, the uid-uid interface is represented by a Lagrangian grid and the flow is solved on a stationary Eulerian mesh. e Lagrangian points are represented by a set of marker points and move with the fluid flow defined at the Eulerian grid. To represent viscosity or density difference between different fluids, and indicator function is calculated based on the marker points [21][22][23]. erefore, the indicator function has an important role in not only theoretical issues coupling with the hydrodynamics equations [24], but also in practical approaches in recent studies. e density field was determined by the indicator function derived from an irrotational discrete delta vector in the droplet simulations with a large density ratio in [25]. To incorporate the viscoelasticity effect of Oldroyd-B fluid, the smoothed Heaviside indicator function was used for simulating the dynamics of Newtonian vesicle in [26]. e role of the indicator function was represented by the solid fraction in the diffuse-interface IB framework for conjugate-heattransfer problems proposed by [27]. Moreover, the indicator function was used to identify the spatial distribution of physical quantities to predict the structure and refractive index profile of fiber-optic components in [28]. In [21], the authors presented the indicator function I(x, y) based on a continuous version as follows: where A is an area. en, symbolically, we have the following equation: where n ′ is the outward unit normal vector at the domain boundary zA. However, it is not straightforward to connect the relationship between the numerical scheme and the continuous version of the indicator function. erefore, the objective of this paper is to present a detailed derivation and numerical investigation of the indicator function for front-tracking methods. e paper is organized in the following manner. In Section 2, we present the detailed derivation of the indicator function. In Section 3, various numerical experiments are performed to investigate the effect of parameters such as distance between points, uniformity of the distance, and types of the Dirac delta functions on the indicator function. In Section 5, conclusions are drawn.

Derivation of the Indicator Function
Let δ(x) be a smoothed Dirac delta function and satisfy the following equation: For example, a 4-point delta function [29][30][31][32] is given as follows: where h is the space step size in the discrete equation and given by the following equation: Now, we consider a series of the shifted delta function δ(x − X 0 ) for some − h < X 0 < h and it is shown in Figure 1 for some X 0 � 0.7 and h � 1, here we take h � 1 for better visualization. We describe the h effect in Section 3.1.
Let x i � ih for i � ± 1, ± 2, . . .. From the definition (4), we have the following equation: For Let us consider a point x � X 0 at the interface in the onedimensional space. We want to construct an indicator function which increases monotonically from zero to one across the point X 0 as x changes its position from left to right of the point, i.e., (7) Figure 2 shows the delta function δ(x) (solid line) and the indicator function Taking a second derivative of equation (7) with respect to x, then we have the following equation: It implies that the indicator function can be found by solving the Poisson equation (8). Let us consider a 4-point delta function with X 0 � 0 and a discretization for equation (8) as follows: where we have used the unit space step size h � 1 and By multiplying both sides by the inverse of the coefficient matrix, we have the following equation: From equation (11), we obtain the following results: where we have used  (12). Figure 3(b) shows discrete delta and indicator functions with various point positions, − 1 < X 0 < 1.
Next, let us consider a hat shaped indicator function as shown in Figure 4 and the sum of the corresponding Dirac delta functions is δ(x + 4) − δ(x − 4). If we solve the following equation for the indication function with zero Dirichlet boundary condition, then we have the result in Figure 4.
Now, we consider this in the two-dimensional space. Let a computational domain Ω � (a, b) × (c, d) be partitioned in Cartesian geometry. Let N x and N y be the number of cells in the xand y-directions, respectively. We assume a uniform mesh with mesh spacing h � where n(x) is the outward normal vector at x. Next, let us rewrite equation (14) as follows: We replace n(x)|∇I(x)| in equation (14). by M l�1 n(X l )δ (x − X l )Δs l , then we have the following equation: where After taking the divergence operator to equation (16), we have the following Poisson's equation: We solve equation (17). with appropriate Dirichlet boundary condition. In this study, we use the Gauss-Seidel method.

Effects of h.
In this section, we describe the effects of spatial grid size h. Figure 10 shows 4-point delta functions and corresponding indicator functions for different spatial step h � 1, 0.5, and 0.25. We observe that as h is smaller, the delta function becomes sharper and the interface transition layer of the indicator function becomes shorter.
In computational domain Ω � (− 15, 15) × (− 15, 15), we consider uniform immersed boundary points for a circle  radius of 10 with different h � 1, 0.5, and 0.25 as shown in Figure 11(a). In Figure 11(b), we solve the indicator function, accordingly the points in Figure 11(a). It has a more stiff interface as h is smaller. Figure 11(c) shows ||I (k) − I (k− 1) || 2 for the Gauss-Seidel iteration with a tolerance tol � 1.0e − 6, where k is the number of iterations and I (0) is a zero vector. We can observe that the number of iterations is getting larger as h approaches zero.

Fluid Application
In this section, we give two examples for the applications of indicator function in multiphase flows and fluid-structure interaction (FSI) problem. First, we consider the two-phase flows problem as the dynamics of a droplet suspended in the simple shear flow. e droplet is between two parallel plates as shown in Figure 12. Let u(x, t) � (u(x, t), v(x, t)) be fluid velocity at x � (x, y) and time t. e Lagrangian variable, denoted as X(s, t), where s is arc length parameter, expresses the immersed boundary. e governing equations on the domain Ω are as follows: where Γ is the interface between the two fluids. We define the Reynolds number in equation (22) as Re � ρ 2 R 2 _ c, We � Re · Ca, and Ca � ηR _ c/σ, where ρ is fluid density, η fluid viscosity, R droplet radius, _ c the shear rate, and σ the surface tension coefficient. We define a gradient field, To compute the indicator function, the following equation is solved: ΔI(x, t) � ∇ · Γ n(X(s, t))δ 2 (x − X(s, t))ds . (25) en, the variable density ρ and viscosity μ are defined as follows: where ρ i and μ i for i � 1, 2 are density and viscosity of fluid i, respectively. For simplicity, the viscosity and density of the droplet is set equal to that of the matrix phase. Our focus is on the generations of indicator functions from the unit normal of the interface in the drop deformation process. We refer to [37] for the details of numerical algorithms. We consider the simulations of drop deformation on Ω � [0, 2] × [0, 2] with a space mesh size 128 × 128. A droplet with a radius 0.5 is put into the center of the domain with initially three sets of different Lagrangian points for comparisons. Delta function  Figure 13. While the drop deforms, the interface is remeshed by adding and deleting points with the condition Δs max /Δs min ≤ 2.0 for different Δs ini leading to the boundary points being more evenly distributed. As the time evolves, the indicator function could be well resolved for Δs ini � 1.0 h as shown in the second row of Figure 13. However, increasing Δs will lead the indicator function a slightly serrated border as shown in the second and sixth rows of Figure 13 (case Δs ini � 5.0 h and 10.0 h), which is unacceptable to separate two immiscible fluids. Moreover, we could see the drop interface shrink in the normal direction quickly by setting the Δs larger while drop deforms in the shear flow, as shown in the odd row of Figure 13. Table 1 lists the area losses as the drop evolves at the time t � 2.5. It shows that the larger Δs, the more the drop area loss.
Next, we take the direct-forcing immersed boundary method as an example application in the FSI problem of the indicator function. e nondimensional governing equations are similar to the multiphase flows and read as follows: where u is the velocity, p is the pressure, u is the body force, and Re is the Reynolds number. From equation (27), we can get the following equation:  to handle the multiphase flows and FSI lies on the how to evaluate the boundary force. For a moving rigid body immersed on the fluid, the imposed force term, the authors in [39] pointed out it consists of two parts that can be written as the following equations (30) and (31): where zΩ is the solid boundary, Ω the solid phase, Δt the time step, I ij the indicator function, u ij , v ij the velocity in Eulerian gird, N x , N y the mesh grid numbers, and V ij the mesh volume in ij grid point. One force in the first terms of right hand sides of equations (30) and (31) is that of a submersed body acting on the external fluid, which can be determined simply by integrating the immersed boundary force densities by equation (29) on the solid boundary that is already evaluated by the direct forcing. Another one that contributes to the unsteady flow inside the solid phase, is the so-called internal or virtual flow force [39] in the second terms of the right hand side of equations (30) and (31). Usually, if the object is fixed in the fluid, we impose the f � (f x , f y ) term to satisfy the no-slip boundary condition on the solid surface leading to the second term equals zero. However, when the solid boundary is moving, we should not neglect the effect of the virtual flow, even though its contribution is small [39]. For a complexly enclosed geometry, the indicator function can be used to calculate to its inside virtual flow by letting its value 1, and outside 0 of the moving solid object.
For this moving solid problem in fluid, we consider a 2D NACA0012 airfoil as the representative of a fish body periodically undulating in the incoming flow. e study of this kind of fish kinematics can be found in [40,41]. Our focus is on the fish shape expressed by the indicator function. As in the first example of drop deformation in the shear flow, we investigate the effect of ds on the indicator function. e computational domain is [0, 4] × [0, 2], with mesh size 256 for y direction. e airfoil is initially represented by 202, 102, and 42 Lagrangian points leading to the point distance Δs is about 0.918h, 2.117h, and 6.134h, where h is the Eulerian grid size. e odd rows of Figure 14 shows the fish undulations with outer normal vectors attached on the boundary points for different Δs settings. Because the boundary deforms not too much, we do not need to remesh it. As shown in Figure 14, by increasing Δs, the indicator functions get zigzag distribution on their boundaries between 0 and 1, which could lead to the calculations of f x and f y in equations (30) and (31) a little oscillations. However, the Δs ≈ 1.0h is good enough to resolve the smooth distribution of indicator function as in the first drop deformation example. It should be noted that the diffused fish interface is not well resolved at the high curvature section of the fishtail for larger Δs as shown in the third case of Figure 14. Letting Δs ≈ 1.0h resolves it quite well.

Conclusion
In this article, we presented a detailed derivation and numerical investigation of an indicator function together with its fluid applications in the immersed boundary method. We used the discrete Dirac delta function to construct an indicator function from a set of Lagrangian points. We solved the resulting discrete Poisson equation with the zero Dirichlet boundary condition using an iterative method. We presented several computational tests to study the effect of parameters such as distance between points, uniformity of the distance, and types of the Dirac delta functions on the indicator function. Two applications indicate that the construction of the indicator function are effective, multipurpose, and easy-to-use in separating the changing phases. We note that the indicator function can be extended directly and readily to the three-dimensional space. Based on the above-given numerical investigations, as future works, we

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that there are no conflicts of interest.