Immersed Boundary Method for Simulating Interfacial Problems

: We review the immersed boundary (IB) method in order to investigate the ﬂuid-structure interaction problems governed by the Navier–Stokes equation. The conﬁguration is described by the Lagrangian variables, and the velocity and pressure of the ﬂuid are deﬁned in Cartesian coordinates. The interaction between two different coordinates is involved in a discrete Dirac-delta function. We describe the IB method and its numerical implementation. Standard numerical simulations are performed in order to show the effect of the parameters and discrete Dirac-delta functions. Simulations of ﬂow around a cylinder and movement of Caenorhabditis elegans are introduced as rigid and ﬂexible boundary problems, respectively. Furthermore, we provide the MATLAB codes for our simulation. choice of the discrete delta function affects the dynamics. The error estimation of regular functions is a topic that has been actively discussed. We also showed that the use of the conventional IB method is sufﬁcient to identify the behavioral characteristics of ﬂows around complex boundaries by considering the effect of the Reynolds number. Nevertheless, it remains an open problem to simulate ﬂows around complex boundaries in the case of high Reynolds numbers.


Introduction
Fluid-structure interaction is widely observed in various fields [1][2][3][4][5]. However, simulating fluid-structure interaction is a challenge problem in mathematics, because of its boundary layer transition. The immersed boundary (IB) method was originally proposed by Peskin in 1977 in order to model flow patterns around heart valves [6,7]. It is widely used for simulating fluid-structure interactions, such as turbulent flow [8,9], biomembrane [10,11], glioma invasion [12], and swimming eel [13]. In the IB method, the Eulerian variables are defined based on a fixed Cartesian coordinate system, and the Lagrangian variables move freely through the fixed Cartesian coordinates. Their interaction is expressed by a discrete Dirac-delta function. The advantages of the IB method are its simplicity and efficiency in handling complicated geometries without mesh regeneration.
In the original IB method, the massless boundary moves with the velocity of the ambient flow by the fluid force, and the density of the structure is spread out to the ambient flow in order to consider the inertia of the elastic structure [6]. The movement of the elastic boundary with mass can be easily simulated by the extension of the original IB method by Kim and Peskin, which is called the penalty IB method [14]. The penalty IB method has recently been extended to simulate the elastic boundary with mass [14] and the interaction between a rigid body and a surrounding fluid [15]. The numerical simulation with high Reynolds numbers is strictly limited when using the conventional IB method, because it is needed to capture the flow transition in a thin boundary layer at high Reynolds numbers. It can be resolved by refining the mesh grid locally [16,17] or imposing a velocity profile based on a wall model [18,19] at the cutting egde. This paper briefly reviews the IB method and its numerical implementation. The derivation of the equation of motion for the IB method is quite intuitive; however, some issues remain regarding its implementation for beginners, such as the choice of delta functions and the effect of parameters. To resolve this, the actual program code is provided with the corresponding numerical simulations. Moreover, a practical example is introduced in order to explain the modeling of biological phenomena using the IB method.
The remainder of this paper is organized, as follows. We describe the mathematical formulations for the IB method in Section 2. Section 3 presents the discretization of the IB method. The numerical results are presented to show the effect of the Reynolds number and discrete Dirac-delta function in Section 4. Finally, Section 5 presents the conclusion. Furthermore, we provide the MATLAB code for the numerical implementation from the corresponding author's webpage and its structure can be found in Appendix A.

Mathematical Formulations
Let X(s, t) = (X(s, t), Y(s, y)) be the position of the Lagrangian points with label s, describing the configuration in a two-dimensional space at particular time t. If we assume that the evolution of the configuration is determined by the (elastic) energy according to the functional E[X], the functional is expressed in the following form: where ℘ is a variation symbol that distinguishes the symbol from the Dirac-delta function δ. Taking the Fréchet derivative of E at X(s, t), we can define the force density that is generated by the material F(s, t) as follows: The elastic energy functional is given by where σ is the stiffness coefficient and | · | is the L 2 -norm. Subsequently, we obtain Thus, the elastic force density can be rewritten as where is the tension of the configuration and τ(s, t) = ∂X/∂s |∂X/∂s| , is the unit tangent vector to the configuration. This is the only formulation for the elastic force, and it can be modified by adding other energy functional to Equation (3) in order to consider the corresponding force.
Let u = (u(x, t), v(x, t)) be the velocity at the position x and time t. Subsequently, we can define the velocity of the configuration U(s, t), as follows: where δ(x) denotes the two-dimensional Dirac-delta function δ(x)δ(y). Note that the material derivative of u is defined as and it implies that Conversely, the system would naturally evolve to minimize the time integration of the difference between the kinetic and potential energy functionals: where v(x, t) is the variational derivative of the configuration at position x and time t, which implies v(X(s), t) = ℘X(s, t). Because the expression contains both Lagrangian and Eulerian variables, we should eliminate the Lagrangian variables by defining the fluid mass density ρ(x, t) and fluid force density f(x, t), as follows: f(x, t) = F(s, t)δ(x − X(s, t))ds.
For simplicity, we assume a constant mass density, which is, M(s) ≡ M and ρ(x, t) ≡ ρ. The equations of fluid motion are governed by the Navier-Stokes (NS) equation, as follows: where p(x, t) is the pressure, µ is the viscosity, and Ω = (0, L) 2 is the square domain. Let Re = ρU * L * /µ and We = ρ(U * ) 2 L * /σ be the Reynolds and Weber numbers with the characteristic velocity U * and length L * , respectively. Subsequently, the non-dimensionalized version of Equations (14) and (15) can be expressed, as follows:

Discretization
Various numerical schemes have been developed to solve the NS equation. In this study, we employ the implicit scheme, as follows [20]: Here, ∇ h is the discrete gradient operator using the central difference and ∆ h is the discrete Laplacian operator defined as For the convection term u n · ∇ d u n = (uD ± x u + vD ± y u, uD ± x v + vD ± y v), we use the skew symmetry, as follows: Other terms are computed in a similar manner.
Applying the projection method [21], we split Equation (18) into two equations with the temporal velocityũ, as follows:ũ Taking the divergence operator into Equation (24), the Poisson equation for the pressure is derived by Equation (19), as follows: Equation (23) can be solved directly, because it has an explicit formulation, and the fast Fourier transformation is used to solve Equations (24) and (25). The details can be found in [20]. Here, the fluid forcing term f n is given by where G h is the set of indices for the Lagrangian points, δ h is the discrete Dirac-delta function, which was discussed in the later section, and ∆s is the step size of the label of the configuration. Here, the discrete force density for the structure F n (s) is defined by where D s is the discrete differential operator with respect to s. Finally, the discretization of the transport equation for IB is given, as follows: where g h is the set of Eulerian points.

Code and Simulation
In this section, we discuss the MATLAB code of the IB method for solving a moving boundary problem with a simple and basic example. The code updates the fluid velocity u n+1 and boundary position X n+1 from the given values u n and X n .
To explain how to obtain the updated boundary position, we first set a domain of fixed Cartesian coordinates for a fluid and let the IB structure on the fixed domain. Here, the following ellipse is chosen for the initial condition of the boundary: where L is the length of the domain. Next, the elastic boundary force density is evaluated from the IB structure. It can be obtained from a minus gradient of the elastic energy functional with respect to X. An elastic energy function can be classified based on three types of contributions: stretching, bending, and tether energies. In the example, the force is only generated from the stretching energy with a zero rest length. 1 %% generate the boundary force dencity 2 F=(X(kp,:)+X(km,:)-2 * X)/(ds * ds); Using the evaluated boundary force density in the previous step, we spread out the boundary force density to a fluid force density while using the Dirac-delta function. This is one of the key steps for the fluid-structure interaction in the IB method. There are some choices of the Dirac-delta function, and its details will be discussed in the next section.  1)). * phi2(r (2)); 10 f(i1,i2,1)=f(i1,i2,1)+(c * F(k,1)) * w; 11 f(i1,i2,2)=f(i1,i2,2)+(c * F(k,2)) * w;  1)). * phi2(r(2)); 9 U(k,1)=sum(sum(w. * u(i1,i2,1))); 10 U(k,2)=sum(sum(w. * u(i1,i2,2)));  Figure 1 shows the vorticity and pressure contours around the damped vibration elastic body after code execution. Note that the vorticity is a physical quantity that represents a local circulation in the fluid and it is defined as a flow rate rotation, and the pressure difference is derived from the boundary force density. If the reader uses the code to simulate other problems, then it is important to modify the initial boundary and force density generation. We will introduce the simulation of the rigid and flexible boundary problems in Sections 4.3 and 4.4, respectively.

Discrete Delta Function
The Dirac-delta function, δ(x), plays an important role in the interaction of fluids and boundaries in the IB method. In Equations (8) There are some common choices for the discrete delta functions or regular functions φ i , such as 2-point, 4-point, 4-point cosine, and 6-point functions. Figure 2 shows the profile of the various regular functions, φ i , i ∈ {2, 4c, 4, 6} that form the delta functions, δ i . Among the choices, the only regular function φ 4 , which satisfies all of the above five conditions, is defined, as follows: Other regular functions can be obtained by changing the support size or adding/deleting conditions. For example, the four-point cosine approximation does not satisfy the second moment condition (5th condition in the list). The detailed formulas mentioned above can be found in [23]. We constructed a dumbbell-shaped boundary, as shown in the upper left panel of Figure 3, to investigate how the selection of the delta function affects the simulation. The distance between boundaries is equal to h/4 at X = 2. The initial boundary comprises a single closed curve and it is given a stretching force. Thus, with time, the boundary must change to a circular shape. Figure 3 shows the simulation results of the selection of the delta function in the IB method. When using δ 4c , it is the most affected by close distance, and when using δ 2 , it is the least affected. The delta function δ 6 has wider support than δ 4c . However, the actual impact thickness shows that function δ 6 is narrower than function δ 4c from the simulation results. Figure 4 shows the area within the boundary over time in each case that is covered in Figure 3. In general, the traditional IB method has a disadvantage of causing volume loss within boundaries. The use of the δ 2 function in this simulation result appears to be the most vulnerable to volume retention within boundaries.

Effect of the Reynolds Number
Next, we apply the IB method for modeling flows around fixed rigid bodies. A simple approach for implementing a fixed rigid body using the traditional IB method is to place the target point(X T ) at the desired location and attach the boundary(X) whlie using a rigid spring. Here, the tethered force density is defined by where k T is the stiffness constant. The flow around obstacles, such as circular or square cylinders, has been extensively studied for Newtonian fluids [2,5,[24][25][26]. Flows around these obstacles exhibit a variety of phenomena, such as boundary layer separation, vortex shading, lift, and drag forces. In general, the flow around a circular cylinder can be categorized by the range of the Reynolds number, as follows: • Re ≤ 5, regime of unseparated flow, • 10 ≤ Re ≤ 40, a fixed pair of Föppl vortices in the wake of the cylinder, • 40 ≤ Re ≤ 150, periodic the von Karman vortex street, and • 150 ≤ Re ≤ 300, transition to turbulence in vortex street. Figure 5 shows vorticity and streamline distribution for fluid flow across 2D cylinders, depending on the various Reynolds numbers. The flow around the circular cylinder according to the Reynolds number is consistent with the literature. We can see that vortex shading is produced behind the obstacle when Re = 100 and 200. Moreover, the length of the vortex shedding and distance between Karman vortices in the square cylinder appears longer than in the circular cylinder. Figure 6 shows the drag coefficients of circular and square cylinders as a function of the Reynolds numbers. The drag coefficient is a dimensionless number used to quantify the drag of an object in a fluid environment. The higher the Reynolds number, the smaller the drag coefficient. The drag coefficient is higher in the square cylinder than in the circular cylinder. Table 1 lists the drag coefficients and lift coefficients calculated for Re = 100 and 200. The value of the present result is higher than that of the references. This is attributed to constraints caused by treating rigid bodies while using the IB method. In the traditional IB method, the stiffness constant was set high to indicate a fixed rigid body in response to the flow. This makes numerical calculations unstable and makes it difficult to calculate the correct coefficients. Despite these limitations, simulating the rigid body using traditional IB methods has the advantage of dealing with complex boundaries. Figure 7 shows the flow around the obstacle, which includes not only circular and rectangular cylinders, but also star-shaped boundaries of complex boundaries. Looking at the flow around the star-shaped boundary, we can see that the pointed part of the shape affects the flow. In other words, it is also sufficient to use traditional IB methods to identify the behavioral characteristics of flows around complex boundaries.

Biological Problem
Simulations of eel, sperm, bacteria, and other swimmers have been developed whlie using the IB method. In this section, we discuss the simulation results of various motions of the C. elegans and present a simulation of the free swimming using the randomwalk. C. elegans is a nematode of approximately 1 mm length that lives in temperate soil environments. To simulate the movement of the C. elegans, we used one line representing the C. elegans body as the center skeleton. Two forces that produce movement are applied to the center line; the stretching force and bending force. The stretching force is used for representing the unextended body; the bending force is used for representing the body shape. The stretching energy is shown in Equation (3), and the bending energy is shown, as follows: where σ b is the bending coefficient, N is the unit normal vector to the body, and c 0 is the reference wave curvature. Various movements may be indicated, depending on the construction of c 0 . For example, to simulate a forward nematode, c 0 is organized as a traveling wave, as follows: where k = 2π/λ is the wave number with wavelength λ, A is the oscillation amplitude, and ω is the oscillation frequency. Figure 8 shows the simulation results for various movements (forwarding motion, backward motion, rest, and coil turn) of C. elegans according to the c 0 configuration. The left side of each panel represents the wave shape, and the right side represents the snapshot of C. elegans over time. Unlike in the top panels, the wave shape in the bottom panels is consistently applied over time. In-silico model of C. elegans combines each of these motions to move forward or rotate to the desired point. Note that changes in wave curvature alone lead to different movements.
Randomwalk is a process that comprises a sequence of steps determined by chance. As shown in Figure 8, we simulated the randomwalk with the combination of each motion of C. elegans. Figure 9 shows the results of the randomwalk. Each motion consists of forward, rest, and coil turns, assuming a motion decision cycle of 0.25 s. In each panel, 'S' represents the starting point and 'E' represents the ending point. In the top panel, it can be seen that the progress of C. elegans changes direction immediately after the coil turn. In the lower panels, the trajectory of motion was compared according to the ratio of determination of moving (forward) and stop (rest and coil turn) motion. Each simulation was performed for 20 s and the ratio of rest and coil turns in stop motion set equal, which is 1 : 1. Table 2 shows the moving distance and turn angle according to the ratio of the moving and stop motions. We can notice that a higher rate of moving motion yields larger displacement, and a higher ratio of stop motion causes more rotation.

Conclusions
In this paper, we briefly reviewed the IB method in order to investigate the fluid-structure interaction problems that are governed by the NS equation and provided the MATLAB codes for numerical implementation. The essence of the IB method is its simplicity to represent the (not limited to simple) interaction between fluid and structure by coupling Eulerian and Lagrangian representations. The equations of motion and discretization were described in the IB framework. The computational results are shown as an example to facilitate the understanding of the IB method, including the choice of discrete delta functions, the effect of the Reynolds number, and the application for the biological problem. In a typical example or the volume conservation, it was shown that the choice of the discrete delta function affects the dynamics. The error estimation of regular functions is a topic that has been actively discussed. We also showed that the use of the conventional IB method is sufficient to identify the behavioral characteristics of flows around complex boundaries by considering the effect of the Reynolds number. Nevertheless, it remains an open problem to simulate flows around complex boundaries in the case of high Reynolds numbers.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

IB
immersed boundary NS Navier-Stokes

Appendix A. Code Structure
The MATLAB codes for our simulations are available from the corresponding author's webpage: https://sites.google.com/view/sglee/research.