A ghost-cell immersed boundary method for unified simulations of flow over finite- and zero-thickness moving bodies at large CFL numbers

A ghost-cell immersed boundary method for unified simulations of flow over finite- and zero-thickness moving bodies at large Courant-Friedrichs-Lewy (CFL) numbers is presented. In order to handle such bodies in a unified manner, algorithms for interface construction and cell demarcation are proposed. The main challenge in treating zero-thickness bodies is to maintain sharpness and accuracy even at large CFL numbers with diminished spurious force oscillations. Thus, the effect of large CFL numbers on the solution accuracy of fluid-structure interaction (FSI) problems involving zero-thickness bodies is investigated and necessary treatments to preserve solution accuracy even at large CFL numbers are suggested. The present study suggests two treatments which are important in preserving the accuracy and stability of the solution: backward time integration for computational cells called ‘swept-cells’ and pressure boundary condition with mass conservation. Composite implicit time integration for the dynamic equation of a thin elastic structure is employed for a stable simulation of FSI at large CFL numbers. By using large time step sizes, the present method not only enhances computational efficiency, but also suppresses spurious force oscillations while maintaining the sharpness of an infinitesimally thin body. The efficacy and accuracy of the present method are examined through numerical examples.


Introduction
In many fluid-structure interaction (FSI) problems found in nature, moving bodies interacting with the fluid are combinations of objects with both finite and infinitesimally small thicknesses. Several examples include fins of a swimming fish (Cui et al., 2020;Zhang et al., 2018), wings of a flying insect (Oh et al., 2020;Yuan et al., 2015;Zhang et al., 2015), valves of a cardiovascular system (Liu & Liu, 2019;Tang et al., 2012), or a wing of a falling maple seed . Such objects are also observed in various engineering applications such as the trailing edge of a heaving and pitching foil (Fairuz et al., 2013;Ghalandari et al., 2019;You et al., 2019), blades of a rotating turbine (Ezhilsabareesh et al., 2018), or suspenders of an oscillating suspension bridge (Kim & Choi, 2019a). Therefore, a unified modeling of both finite-and zero-thickness bodies is important for accurately predicting their dynamics. The main difficulty is treating zero-thickness bodies undergoing large motions, which requires an accurate and computationally efficient numerical method (Salih et al., 2019).
Flow over zero-thickness bodies is usually modeled using continuous forcing immersed boundary (IB) CONTACT Donghyun You dhyou@postech.ac.kr methods (Fauci & Peskin, 1988;Jung & Peskin, 2001;Kim & Peskin, 2006). These methods initially assumed zero-thickness bodies to be massless, but methods that spread mass of the structure to the ambient fluid have been suggested to take the inertia of the structure into account (Kim & Lai, 2010;Kim & Peskin, 2007Zhu & Peskin, 2002). To better handle mass and smearing of force, a feedback forcing approach has been proposed (Goldstein et al., 1993) and has been employed in various simulations of zero-thickness structures (Huang et al., 2007). However, this method requires tuning of two problem-dependent parameters, which affect the accuracy and stability of the solution. In addition to such arbitrariness associated with the parameter tuning, severe restriction in the computational time step size has been a major drawback for practical applications (Goldstein et al., 1993).
To overcome such limitations, discrete forcing IB methods have been also introduced to simulate flow over zero-thickness bodies. These methods are attractive due to relieved restriction in the time step size (Fadlun et al., 2000;Kim et al., 2001;Mittal & Iaccarino, 2005) and sharp representation of IBs (Huang & Tian, 2019;Kim & Choi, 2019b;Mittal & Iaccarino, 2005). However, the methods often suffer from spurious force oscillations on the moving body. In order to reduce the spurious oscillations, Luo and his colleagues Tian et al., 2014;) exploited a flow reconstruction remedy which can sacrifice the sharpness of the interface as noted by Lee and Choi (2015); Seo and Mittal (2011). Lee and Choi (2015) avoided any treatment on the fluid node in the vicinity of the IBs to maintain sharpness, and obtained the hydrodynamic force directly from the Navier-Stokes equations to attenuate spurious force oscillations. This method was able to predict flow over an elastic slender body at moderate CFL numbers (0.75−1.0), but it is only applicable to a slender body having a finite thickness, which requires a number of cell layers within the body. On the other hand, a zerothickness structure can be treated using IB methods by Tullio and his colleagues (de Tullio & Pascazio, 2016;Spandan et al., 2018Spandan et al., , 2017. However, smearing of the interface as well as restrictions in the time step size are not prevented due to the spread of momentum forcing terms to Eulerian fluid cells (Kim & Choi, 2019b). A ghost-cell IB method has been proposed by Mittal and his colleagues (Dong et al., 2010;Mittal et al., 2008;Rips & Mittal, 2019;Zhang et al., 2015;Zheng et al., 2013), which is capable of simulating a zero-thickness structure interacting with the fluid, while maintaining sharpness at the interface.
Although the sharp interface ghost-cell IB method can simulate flow over zero-thickness bodies at moderate time step sizes, increasing the time step size to larger values is necessary for reducing spurious force oscillations (Lee et al., 2011) and improving computational efficiency. Studies on flow over zero-thickness bodies have mostly adopted moderate CFL numbers within 1.0 (de Tullio & Pascazio, 2016;Zheng et al., 2013;Zhou & Mittal, 2020). However, as will be presented in the following sections, we have found that the flow field around zerothickness bodies is significantly influenced by the CFL number, and may even deteriorate at large CFL numbers. To sum up, the main challenge in treating zerothickness bodies is to maintain sharpness and accuracy even at large CFL numbers with diminished spurious force oscillations.
Thus, the aim of the present work is to investigate the effect of large CFL numbers on the solution accuracy of FSI problems involving zero-thickness bodies, and suggest necessary treatments to preserve solution accuracy even at large CFL numbers. The present study suggests two treatments which are important in preserving the accuracy and stability of the solution: backward time integration for computational cells called 'swept-cells' and pressure boundary condition with mass conservation. It is shown that such treatments are essential for simulating zero-thickness bodies at large CFL numbers, which have not been previously reported. Composite implicit time integration for the dynamic equation of a thin elastic structure is employed for a stable simulation of FSI at large CFL numbers. In addition, algorithms for interface construction and cell demarcation are proposed to handle both finite-and zero-thickness bodies in a unified manner.
The paper is organized as follows. In Section 2, numerical methods for simulating FSI problems with finite-and zero-thickness structures at large CFL numbers are presented in detail. In Section 3, the efficacy and accuracy of the present IB method for simulations of flow over finite-and zero-thickness bodies are examined, followed by concluding remarks in Section 4.

Navier-Stokes equations for fluid flow
The non-dimensionalized Navier-Stokes equations are given as follows: where u i is the velocity component in the i-direction.
Variables are non-dimensionalized by a length scale L and a reference velocity U ∞ . Pressure is normalized by ρ f U 2 ∞ where ρ f is the density of the fluid. The Reynolds number is defined as Re = U ∞ L/ν where ν is the kinematic viscosity of the fluid. q in the mass conservation equation (Equation (2)) is the mass source and sink term introduced to better conserve the fluid mass in the vicinity of IBs (Kim et al., 2001). Pressure and velocity are arranged and discretized on a colocated grid. The Crank-Nicolson scheme is used for both convective and diffusion terms. A second-order central difference scheme is used for spatial discretization.

Dynamic equation for the motion of a thin body
The non-dimensionalized governing equation for the motion of a thin body is as follows: where force), and V * = h s1 s2/L 3 (nondimensional volume of each element). Variables are nondimensionalized by the characteristic length L and velocity U. Dimensional form of Equation (3) and meaning of variables are described in the Appendix 1. Note that the non-dimensional tension or shearing coefficient K T ab should be chosen to be large enough to satisfy the inextensibility condition, but also needs to be selected sufficiently small for numerical stability (Huang et al., 2007;Huang & Tian, 2019;Tian et al., 2014).
For computation of Equation (3), the time integration scheme should be carefully selected in order to guarantee numerical stability with large time step sizes. Huang and co-researchers (Huang et al., 2007;Huang & Sung, 2009 integrated the structure equation in time with a backward finite-difference (BFD) scheme, and reported that the balance between K T ab and the time step size is important for the stability. Lee and Choi (2015) employed an implicit second-order finite-difference (ISFD).
However, the effect of numerical integration methods for the structure equation on the stability of coupled FSI simulations is not clearly identified in the previous research where the time step size used was sufficiently small to guarantee the stability. In addition to above-mentioned BFD and ISFD methods, in the present study, other well-known numerical methods for the structure equation such as the Newmark (NWMRK) method (Newmark, 1959), the generalized-α (GENALP) method (Chung & Hulbert, 1993), and the composite implicit time integration (CITI) method (Bathe & Baig, 2005) are compared for the stability characteristics for simulations with large time step sizes. It is found that the CITI method produces the most stable solution, and therefore, the CITI method is used for the forthcoming FSI simulations. Details of the integration methods are summarized in the Appendix 2.

Fluid-structure coupling
The coupling process can be summarized as follows supposing that the time step advances from n to n + 1: (1) At time step n, the hydrodynamic force F n+1 is obtained on structure surfaces X n by interpolating the pressure and calculating the shear stress.
(2) The dynamic equation for the motion of a thin body (Equation (3)) is solved by one of the numerical integration methods described in Section 2.1.2. (3) After moving the position of the body, the cell demarcation algorithm which will be described in Section 2.2 is conducted to determine types of Eulerian cells. (4) X n+1 ,Ẋ n+1 , andẌ n+1 are provided to determine u n+1 and p n+1 .
This procedure is repeated until the end of the simulation.
Note that the procedure is weakly-coupled, since the hydrodynamic force F n+1 is calculated using the solution of u n and p n . In general, a weakly-coupled procedure itself has stability issues, which may cause numerical instabilities when dealing with a problem involving strong added mass effects (low solid-to-fluid density ratios). In common with other ghost-cell IB methods, the present method is also vulnerable to instabilities caused by strong added mass effects. However, the numerical stability discussed in the present study is a separate issue arising from thin bodies moving at large time step sizes rather than the coupling procedure. Thus, the present study proposes two special treatments that will enhance solution accuracy and stability for thin bodies moving at large time step sizes. For this reason, examples having weak added mass effects are selected so that the effect of the coupling procedure could be excluded.

Interface construction and cell demarcation algorithms for unified simulations of finite-and zero-thickness bodies
The first prerequisite for FSI problems using IB methods is to construct the interface and determine cell types of the computational cells so that appropriate boundary conditions are exchanged at the interface. In order to impose boundary conditions from solid to fluid and vise versa, immersed boundaries are often represented with a number of triangular elements called markers (Kim et al., 2018;Lee & You, 2013;Mittal et al., 2008). The approach to treat both finite-and zero-thickness bodies in a unified way is illustrated in Figure 1. Figure 1(a) describes a blunt circle which is continuously modified to a zero-thickness plate. Note that arrows represent outward normal vectors to the fluid. From the illustration, a zero-thickness body is inherently identical to a finitethickness body meaning that they can be treated in the same manner. The key idea is that, in the pre-processing stage, the surface elements for a zero-thickness body are divided into markers which are overlapped but with outward normal vectors in the opposite directions as shown in Figure 1(b). The idea of using the overlapped markers for representing a zero-thickness body is different from the previously reported approaches which employed an artificial thickness (Luo et al., 2012), a thin but finite-thickness (Lee & Choi, 2015), or a single layer of non-overlapped markers (de Tullio & Pascazio, 2016;Gilmanov et al., 2015). The present approach has an advantageous feature for cell demarcation which is important in the ghost-cell IB methods.
In the previous ghost-cell IB methods (Lee & You, 2013;Mittal et al., 2008), once the source and sink cells that are cut by the markers are found, solid-fluid cell demarcation is conducted using dot-product, and then ghost-cells are determined as solid cells that contact at least one fluid cell. In the present unified ghost-cell IB method, demarcation of solid, fluid, and ghost cells is conducted using the ray-tracing. The detailed algorithm for ray-tracing cell demarcation is as follows: first, the source and sink cells and the neighbor cells which are the cells a layer apart from the source and sink cells are selected. Those cells are, in our parallel computing framework, divided into several subdomains which need to be searched by the ray-tracing. The ray is the line between two cell-centers as depicted in the lower part of Figure 2. Based on the information of the current cell, the type of the next cell is determined. If the ray crosses an even number of markers including zero, the type of the next cell is the same with that of the current cell. If the ray crosses an odd number of markers, the type of the next cell is different to that of the current cell. If the current cell is a fluid cell and the number of markers penetrated by the ray is non-zero, the next cell becomes a ghost-cell.
This searching process continues along the snake-like pattern until the pattern reaches the end of the subdomain (see the upper part of Figure 2). In Figure 2, S and E denote the start and the end of the pattern, respectively. The snake-like pattern between each plane is connected by A − A and B − B so that the whole subdomain is connected with a single line. This algorithm is more robust than dot-product and valid for a body with highly curved surface or zero-thickness surface without exception. Figure 3 shows a few degenerate cases produced by dot-product and how those cases are avoided by adopting the present ray-tracing algorithm. Note that the two parallel lines represent two layers of markers as described before. As shown in the figure, a vector connecting the cell-center and the nearest marker-center and the outward normal vector are used in dot-product. If the dot-product produces a positive (negative) value, the cell is considered to be inside (outside) the solid body. However, as can be seen in Figure 3(a), the cell A and the cell B, which are positioned near the highly curved surface, are determined as a solid cell and a fluid cell which in fact are a fluid cell and a solid cell, respectively. Moreover, cell C cannot be determined for the type for a zero-thickness body because there are two normal vectors with opposite signs. On the other hand, as can be seen in Figure 3(b), these degenerate cases are avoided through the ray-tracing algorithm used in the present study. Unlike in the previous ghost-cell IB methods (Lee & You, 2013;Mittal et al., 2008) where fresh-cells are determined as the cells that are previously solid cells, but currently become fluid cells due to a motion of a body, in the present method, a swept-cell concept is newly introduced to define cells which were swept by markers during a time step advancement to consider motions of both finite-and zero-thickness bodies.

Backward time integration for swept-cells
The first requirement for preserving solution accuracy at large CFL numbers is to obtain the solution at the swept cells by taking account of the trajectory of markers. One of the important topics for simulation of flow over a moving body using a ghost-cell IB method is how to assign appropriate values to fresh-cells. A fresh-cell is defined as the cell which used to be a solid cell in the previous time step but has become a fluid cell in the current time step due to the motion of a body. Fresh-cells have no values in the previous time step because the cells were inside the IBs. For this reason, values at fresh-cells have been commonly determined by interpolating velocities at neighbor cells. However, when the fluid cells are swept by a thin body, they contain values of the previous time step, which is why no special attention has been paid to such swept cells. Swept cells were either considered as regular fluid cells (de Tullio & Pascazio, 2016) or treated as fresh cells for a thick body by assuming an artificial (Tian et al., 2014) or finite thickness (Lee & Choi, 2015).
In the present study, the fluid cells swept by both zeroand finite-thickness bodies are treated in a consistent manner. Lee and You (2013) reported and overcame the instability when a large time step size or CFL number is used in the simulation of moving body due to the occurrence of multiple layers of fresh-cells by developing a backward time integration method for fresh-cells. The concept of the fresh-cells is extended to handle not only the cells that change their states from the solid to the fluid but also the cells that swept by the markers during the motion of a body as shown in Figure 4.
Once the swept-cells are defined, the backward time integration developed by Lee and You (2013) is employed. The method is known to be stable even though multiple layers of fresh-cells are generated by the use of large CFL numbers and this advantageous feature is retained even for a zero-thickness body.
Once a swept-cell and the corresponding sweeping marker are selected, assuming a linear velocity profile on the surface of the marker, the pressure gradient at time step t n+α can be calculated as follows: where 0 ≤ α ≤ 1. Note that the right hand side of the Equation (4) is the material derivative which includes the convection term. On the other hand, the velocity vector at time step t n+α can be obtained by using the velocity vector of the marker. Thus, the discretized momentum equations are integrated using the backward time integration scheme to assign velocity vectors to the swept-cells at t n+1 as follows: where represent the convection and the diffusion terms, respectively, and t α = t n+1 − t n+α . U j in the convection term refers to the face velocity. The merit of introducing the swept-cell is that the solution retains its sharpness even with the use of a large time step size while the solution is deteriorated without such treatment as demonstrated in Figures 15 and 20 in Sections 3.2 and 3.3, respectively.

Pressure boundary condition with mass conservation
Another treatment required for preserving solution accuracy at large CFL numbers is to strictly conserve mass via pressure boundary condition with mass source and sink. While simulating flow over a moving body, effects of explicit or implicit imposition of the pressure boundary condition at the interface on the result are rarely reported in the literature. However, it is found that, in the present study, as the thickness of the body becomes thinner, explicit imposition of the pressure boundary condition becomes more important, especially to maintain the sharpness of the interface when a large time step size is employed.
whereû k is the intermediate velocity. As shown in the Figure 5(a), when the pressure boundary condition is not explicitly imposed, the Poisson equation is solved by using the same stencils as in Equation (6). On the other hand, as shown in the Figure 5(b), when the pressure boundary condition is explicitly imposed, the pressure gradient terms (δp/δx)| i−1/2,j and (δp/δy)| i,j−1/2 are substituted by −du x,IB /dt and −du y,IB /dt, respectively, by assuming a linear velocity profile on the surface of the IBs as in Section 2.3.1. The significance of explicit imposition of the pressure boundary condition to retain sharpness of the interface, especially for a large time step size and for a reduced thickness is shown in Figure 12 of Section 3.2.
In addition to the left-hand side of the Poisson equation (Equation (6)) where the pressure boundary condition is concerned, the right-hand side of the equation should be treated properly for simulating zero-thickness bodies using large CFL numbers. The conservation of fluid mass in the vicinity of IBs is strongly enforced by employing the mass source and sink algorithm. The present mass source and sink algorithm is based on the method proposed by Kim et al. (2001), in addition, due to the existence of a zero-thickness body, the virtual merging used in Seo and Mittal (2011) is also exploited. In the previous ghost-cell IB methods (Lee & You, 2013;Mittal et al., 2008;Seo & Mittal, 2011), computational cells containing both fluid and solid fractions are denoted as the source and sink cells. However, in the presence of zero-thickness bodies, the source and sink cell can contain only fluid cut by the markers.
As shown in Figure 6(a), in a two-dimensional configuration for instance, the mass source and sink term per unit span can be expressed in a discrete form as follows: where x s , y s , and l IB are solid fractions of cell faces in each direction and the surface area of the immersed body cut by the cell faces per unit span. U IB is the velocity of the surface of a moving body which is zero for a stationary body. The original mass flux of a control volume shown in Figure 6(b) is modified by excluding the mass source and sink term of the solid fraction ( Figure 6(a)) so as to solely calculate the mass flux of the fluid fraction shown in Figure 6(c). As shown in Figure 6(c), the meaning of Equation (2) is nothing but the mass flux of the fluid fraction normalized by the control volume and it corresponds to finite volume formulation of the mass conservation equation. Therefore, it is identical whether calculating When zero-thickness bodies are included, the latter which directly calculates the fluid fraction of the mass flux is more convenient since there are two fluid domains in   the same cell as shown in Figure 6(d-e). Superscripts f 1 and f 2 in Figure 6(d-e) represent fluid fractions 1 and 2 discriminated by the zero-thickness boundary.
When bodies are immersed in the domain, computational cells can be divided into non cut-cells and cut-cells (or mass source and sink cells) which can be classified further into small cut-cells and regular cutcells depending on the solid volume ratio. In order to avoid ill-conditioned matrices and inconsistent number of cells to be solved in each time step occurred in the cut-cell method, a 'virtual' cell-merging technique which transfers source and sink terms of small cut-cells to adjacent non cut-cells or regular cut-cells was proposed (Seo & Mittal, 2011).
For finite-thickness bodies, the virtual cell-merging technique may not be necessary depending on the numerical method for the pressure Poisson equation. However, the virtual cell-merging technique is inevitable to embrace infinitesimally thin bodies and to conserve mass at the same time. In the present study, the technique is generalized to be consistently applicable to both zerothickness and finite-thickness bodies. Figure 7(a-b) show virtual transfer of the source and sink terms of small cutcells to adjacent non cut-cells or regular cut-cells for a blunt body and a zero-thickness body, respectively. Note that, for a zero-thickness body, the small cut-cell and the regular cut-cell exist in the same cell.
The present unified algorithm can be used even though both zero-and finite-thickness bodies are present simultaneously in a single computational domain. For each mass source and sink cell, directions of source and sink terms to be transferred is pre-determined as outward normal from the regular cut-cells while executing cell demarcation algorithm, as shown in Figure 7(c). Therefore, in the virtual merging step, the source and sink terms are easily transferred to the pre-defined directions.
By combining the mass source and sink algorithm and the virtual merging technique for a zero-thickness body, the present IB method can retain the mass conservation property for a moving body with both finite and zero thicknesses (demonstrated in Figures 14 and 19 in Sections 3.2 and 3.3).

Results and discussion
The advantage and accuracy of the present numerical method is verified through various examples. All of the cases considered in the following sections show that the present IB method is capable of accurately simulating FSI phenomena for moving finite-and zero-thickness bodies with large CFL numbers, thereby enhancing computational efficiency and diminishing spurious force oscillations while retaining the sharpness of the IB.

A hanging filament without fluid flow
A hanging filament without fluid flow under a gravitational force is considered to investigate the effect of temporal integration methods on the time step size for simulation of the dynamic equation of a thin body (Equation (3)). In this example, only the Lagrangian coordinate s1 is used. Note that, in Equation (3), the term F H becomes 0 without fluid flow and dividing all the remaining terms with ρ * yields the non-dimensional dynamic equation of a thin body without any fluid variables as follows: where K T 11 = κ T 11 /ρ s U 2 h and K B 11 = κ B 11 /ρ s U 2 L 2 h. In the following simulations, ρ s h = L = U = 1 and N e = 100 are used. θ is the initial angle of the filament with respect to the vertical axis.  (9)); the present solution with the numerical inextensible condition (Equation (8)): In order to validate the accuracy of the solver, K B 11 = 0 is adopted as in the analytical solution derived by Huang et al. (2007). Fr = 10 and θ = 0.01π are considered. The analytical solution of the non-dimensional transverse displacement component of the filament Y(s * 1 , t * ) is as follows: where J i is the Bessel function of the first kind of order i, z i is the i-th positive root of J 0 (z). As reported by Huang et al. (2007), modeling the inextensibility condition requires a sufficiently high tension coefficient, K T 11 . Thus, two different tension coefficients K T 11 = 1000 and K T 11 = 10000 are used to investigate the effect of a high tension coefficient on accuracy. This is shown by a comparison between the analytical solution (Equation (9)) and the present numerical result (Equation (8)) in Figure 8. Although some oscillations are observed at K T 11 = 1000 (blue line of the figure), both cases agree well with the analytical solution.
The effect of the five temporal integration methods described in Section 2.1.2 on the numerical stability is examined. Huang et al. (2007) investigated, for the IFD method with K B 11 = 0.01 and θ = 0.1π, the relationship among the time step size, the elastic coefficients, and the stability. Table 1 shows that the CITI method is capable of maintaining the stability with the largest time step size among the methods which is an advantageous feature when coupled with the fluid solver with a large CFL number. Figure 9(a) shows time traces of the freeend displacement predicted by different methods. Note that solutions plotted are obtained using the maximum time step size allowing a stable solution for each method. Figure 9(b) shows the motion of a hanging filament from 14.1L/U to 14.9L/U. In the following sections, the CITI method is employed to solve dynamic motions of thin bodies.

Flapping filament in a free stream
A flapping filament in a free stream is chosen to evaluate the capability of the present method for simulation of a zero-thickness body and to identify the effect of the explicit imposition of the pressure boundary condition, the mass source and sink with the virtual merging and the swept-cell treatment on the solution especially at a large time step size. The size of the computational domain is [−2L, 6L] × [−4L, 4L] as shown in Figure 10 where L is the length of the filament. A filament is pinned at the leading-edge (X| s 1 =0 = (x l , y l ) = (0, 0), ∂ 2 X/∂ 2 s 1 | s 1 =0 = 0) with the initial inclination angle of 0.1π radians. Free-end boundary conditions (∂X/∂s 1 · ∂X/∂s 1 | s 1 =L = 1, ∂ 2 X/∂ 2 s 1 | s 1 =L = 0, ∂ 3 X/∂ 3 s 1 | s 1 =L = 0) are imposed at the trailing edge. Grid spacings of x = y = 0.005L are uniformly distributed in the vicinity of the flapping filament while it is gradually stretched to the boundaries of the domain. The numbers of grid points are 496 and 349 in the streamwise (x) and transverse (y) directions, respectively. Dimensionless parameters for the filament are chosen as ρ s /ρ f = 100, K T 11 = 2500, K B 11 = 0.1, and Fr = 0.5. Note that, unlike other parameters, K T 11 is an arbitrary value for satisfying the inextensibility condition. In the present simulation, K T 11 of 2500 is chosen to guarantee extension of the filament length less than 1% throughout the simulation. The number of elements N e of the filament is fixed to 100 which makes the length of each element as 0.01L so as to have sufficient resolution for interaction with the fluid. Figures 11(a-b) shows time traces of drag and lift coefficients of the filament compared with those reported in de Tullio and Pascazio (2016); Lee and Choi (2015) where drag and lift coefficients are defined as C D = 2F D /ρ f U 2 A f and C L = 2F L /ρ f U 2 A f , respectively. F D , F L , and A f are the drag force, the lift force, and the frontal area. Figure 11(c) shows contours of the instantaneous vorticity around the flapping filament at four different time steps in a cycle. Vortices are shed alternatively in accordance with the flapping motion of the filament. Note that a favorable agreement of the results among researchers is observed while CFL numbers employed for those simulations are significantly different. Figure 12(a) shows time histories of lift coefficients as a function of the filament thickness and the time step size without explicit imposition of the pressure boundary condition. The black dashed line refers to the result obtained with a zero-thickness filament by imposing the  Although the method of Huang et al. (2007) is able to deal with a zero-thickness body by using the feedback forcing method for the fluid and the BFD method for the structure, severe restriction in time step sizes for both fluid and structure solvers was reported. de Tullio and Pascazio (2016) also deal with a zero-thickness body; however, the CFL number was limited to 0.1 due to interface smearing originated from the distribution of the momentum forcing to nearby fluid cells (Kim & Choi, 2019b). The present method is capable of handling a zerothickness body even with large time step sizes while maintaining the sharpness of the IB by using the unified fully implicit ghost-cell method for a fluid and the CITI method for a structure. The merit in dealing with a zero-thickness body independently of the grid spacing is that one only needs to consider physical parameters such as the Reynolds number. As shown in Figure 13, in order to simulate the flapping filament with the thickness of h/L = 0.01, Lee and Choi (2015) used the grid spacing of x/L = 0.005 to have two grid cells in the thickness direction whereas the equivalent result can be obtained with the grid spacing of x/L = 0.04 in the present study. The advantage of using a large time step size for a moving body with a finite-thickness in reducing spurious force oscillations as well as enhancing computational efficiency has been reported in the previous literature (Lee et al., 2011;Lee & You, 2013). In the present method, the advantageous feature is maintained even for a zero-thickness body without losing sharpness by the swept-cell treatment, modification of the mass source and sink and virtual merging algorithms, and explicit imposition of the pressure boundary condition as discussed in Sections 2.3.1 and 2.3.2. Figure 14 shows the effect of the mass source/sink algorithm combined with virtual merging by plotting time traces of the lift coefficient for three different time step sizes. Note that the time step size increases from Figure 14(a-c). Three cases are considered. Firstly, the red line refers to the case without the mass source/sink algorithm, which shows the largest spurious force oscillation and eventual divergence of the solution. Secondly, the blue line refers to the case with the mass source/sink algorithm, but without virtual merging. This case shows a reduced spurious force oscillation at t = 0.0005L/U (Figure 14(a)) and maintains stability at a larger time step size of t = 0.0015L/U (Figure 14(b)) compared to the first case. Nevertheless, the solution blows up for a much larger time step size of t = 0.0025L/U (Figure 14(c)). Lastly, the black line refers to the case with the mass source/sink algorithm combined with virtual merging. As a result of this algorithm, spurious force oscillations are further improved, and the solution does not blow up even at t = 0.0025L/U. The mass source/sink algorithm combined with virtual merging not only relieve spurious force oscillations but also improve stability at a large time step size. Moreover, the spurious force oscillations are found to be further reduced by increasing the time step size. The amount of reduction of spurious force oscillations is C L rms ∼ (1/ t) 0.6 as shown in Figure 14(d), which is similar to the result reported for an oscillating cylinder by Seo and Mittal (2011).
The effect of the swept-cell treatment on the time history of the lift coefficient is also examined in Figure 15. Without the swept-cell treatment, the lift coefficient is heavily affected by the time step size; however, with the swept-cell treatment, the lift coefficient shows little change as the time step size is increased. Note that the peak values of the lift coefficients are not influenced by the change in time step size. This is a feature different from other IB methods with momentum forcing, which generally show a dependance of peak values on the time step size. This means that the swept-cell treatment plays a significant role in the solution accuracy especially when the number of swept-cells is large due to the use of a large time step size.
In the present study, the simulation with tU/L = 0.0075 which corresponds to the maximum CFL number of 4.0 produces an accurate result whereas CFL numbers of 1.0 and 0.1 are used in de Tullio and Pascazio (2016); Lee and Choi (2015), respectively. For time integration of the structure equation, tU/L = 0.0025 was the maximum time step size for the SECIMP method whereas tU/L = 0.0075 was the limit for the CITI method. This result again stresses that the numerical integration method for a structure must be carefully selected for an FSI simulation of a thin body with a large time step size.

Hovering and flapping flexible wing
A hovering and flapping flexible wing is considered to examine the capability of the present method for simulating a case with active motion as well as passive deformation of a body. The present case was also considered by Lee and Choi (2015); Yin and Luo (2010). The computational configuration is shown in Figure 16  quiescent and all boundaries are treated with the Neumann boundary condition.
The flexible wing with the length of L is clamped at the leading-edge and has prescribed translational and rotational motions. The hydrodynamic, elastic, and inertia forces govern the passive deformation of the flexible wing. The prescribed active motions at the leading-edge consist of a translational motion (x l (t) = (A 0 /2)cos(2πft), y l (t) = 0) which governs the position (X| s 1 =0 = (x l , y l )) and a rotational motion (α(t) = α 0 + βsin(2π ft)) which determines the angle (∂X/∂s 1 | s 1 =0 = (cosα(t), sinα(t))), where A 0 = 2.5L, α 0 = π/2, and β = −π/8. The position and the angle define the clamped boundary condition at the leading-edge. The freeend boundary conditions are used at the trailing-edge (∂X/∂s 1 · ∂X/∂s 1 | s 1 =L = 1, ∂ 2 X/∂ 2 s 1 | s 1 =L = 0, ∂ 3 X/∂ 3 s 1 | s 1 =L = 0). The Reynolds number Re = U m L/ν is 150 based on the maximum velocity of the leading-edge U m = π A 0 f . The dimensionless parameters are chosen as ρ s /ρ f = 100, K B 11 = 32.298, K T 11 = 4998, h/L = 0.01, and N e = 100. Note that these values are different from those of Lee and Choi (2015), where h/L of 0.06 which allows approximately six fluid cells in the thickness direction is used. In the present study, the thickness h/L of 0.01 is used for the structure equation while the thickness is treated as zero for the fluid simulation.
As shown in Figure 17(a-b), the present method predicts, with the maximum CFL number of 3.0, the temporal variation of the trailing-edge angle and and the lift coefficient comparable to those of Lee and Choi (2015) which were obtained with the maximum CFL number of 0.75. Figure 17(c) shows superposition of the shape of the hovering and flapping flexible wing from 13T to 14T. Figure 18 shows the instantaneous vorticity contours around the hovering and flapping flexible wing at eight different time steps during a cycle. Vortices are shed downstream alternatively which make positive lift forces to the wing.
As in Section 3.2, in order to examine the effects of the mass source and sink algorithm combined with a virtual merging and the swept-cell treatment with the increasing time step size on the numerical result, a rigid wing with the same active motion is also simulated. Figure 19(a-c) shows results obtained with three different time step sizes. As in the previous section, only the case with the mass source and sink and the virtual merging treatments produces an accurate result for large time step sizes. The spurious force oscillations are further diminished by increasing the time step size. Figure 19(d) shows that the amount of decrease in the spurious force oscillations is C L rms ∼ (1/ t) 0.5 in the present case.
The effect of the swept-cell treatment on the time history of the lift coefficient is also examined. Figure 20 shows that, as the time step size increases, deviations among lift coefficients become larger when the swept-cell treatment is not applied while the cases with the sweptcell treatment show almost the same results. Note that the time step size of 0.01L/U m corresponds to the maximum CFL number of 3.0 which is much larger than the CFL number used in other previous studies (Lee & Choi, 2015;. The present result again assures the effectiveness of the swept-cell treatment for an FSI simulation of a zero-thickness body with a large time step size.

Three-dimensional flapping flag in a free stream
To examine the capability of the present method for three-dimensional problems, simulations of a zerothickness flapping flag in a free stream, which was also considered by many researchers (de Tullio & Pascazio, 2016;Huang & Sung, 2010;Lee & Choi, 2015;Li et al., 2019;Tian et al., 2014), are conducted. For the configuration shown in Figure 21, the smallest grid spacings of x = y = 0.01L and z = 0.02L are uniformly distributed in the vicinity of the flapping flag and gradually stretched to the outer boundaries of the domain. The total number of cells is 234 × 209 × 100. The flag has the initial inclination angle of 0.1π radians and is pinned at the leading-edge (X| s 1 =0 = X 0 , ∂ 2 X/∂ 2 s 1 | s 1 =0 = 0) while the other edges are treated with the free-end boundary condition. The inextensible condition (∂X/∂s 1 · ∂X/∂s 1 | s 1 =L = 1, ∂X/∂s 2 · ∂X/∂s 2 | s 2 =0,W = 1), no bending moment (∂ 2 X/∂ 2 s 1 | s 1 =L = 0, ∂ 2 X/∂ 2 s 2 | s 2 =0,W = 0), and no shearing force (∂ 3 X/ ∂ 3 s 1 | s 1 =L = 0, ∂ 3 X/∂ 3 s 2 | s 2 =0,W = 0) in the thickness  The dimensionless parameters of the flag are chosen as ρ s /ρ f = 100, K T ab = 2500 (for a = b), K T ab = 100 (for a = b), K B ab = 0.01 (regardless of a and b). The number of elements for the flag is N e = 50(s 1 ) × 50(s 2 ) which makes the length of each element be 0.02L. Note that the minimum grid spacings are chosen to be twice larger than those of Lee and Choi (2015) since no grid cells are necessary inside the flag. Figure 22(a) shows the time trace of the transverse displacement of the trailing-edge, which agrees well with those of de Tullio Pascazio (2016) ;Lee Choi (2015). Figure 22(b-c) show time histories of the drag and lift coefficients. Agreement to those of de Tullio and Pascazio (2016); Lee and Choi (2015) is satisfactory except some visible discrepancies in the drag coefficient. Similar discrepancies have also been observed by de Tullio and Pascazio (2016). Figure Table 2 summarizes the maximum amplitude of the trailing-edge displacement at the middle of the flag, the Strouhal number, and the time step size reported by several researchers. The maximum displacement and the Strouhal number are in favorable agreement with those of other researchers. However, it is worth noting that the discrete IB methods (de Tullio & Pascazio, 2016;Lee & Choi, 2015;Tian et al., 2014) generally predict higher Strouhal numbers than the continuous IB method (Huang & Sung, 2010). The advantage of the present method is clearly identified by the size of the time step, which can be orders of magnitude larger than those employed in the past research.
The flapping flag under gravity has also been simulated to see the difference in the behavior of the flag and to validate the present method once more. The Froude number taken in the present simulation is 0.05. Figure 24    shows time traces of displacements at three different locations of the trailing-edge, which show good agreements to those of Lee Choi (2015). Differently from the flag without gravity, the flag with gravity loses symmetry in the xy-plane and the lower part of the trailing-edge moves more than the upper part. Figure 25 shows the instantaneous vortical structures and the corresponding flag configurations. The flapping flag under gravity shows more complex behaviors than the flag without gravity which is shown in Figure 23.

Flapping flag with a pole in a free stream
The capability of the present method for simulation of both finite-and zero-thickness bodies can be demonstrated by considering a flapping flag attached to a pole in a free stream. The computational configuration is shown in Figure 26. The size of the computational domain is [−L, 7L] × [−4L, 4L] × [−2.5L, 1.5L]. Since the same domain sizes in the streamwise x and transverse y directions are used and only the size in the spanwise z direction is doubled compared to the case in the previous section, the same grid spacings of x = y = 0.01L and z = 0.02L yield 234 × 209 × 200 cells. The time step size tU/L is 0.01 which corresponds to the maximum CFL number of 1.9. Note that all the initial and boundary conditions of the flag are the same with those for the flapping flag with gravity in the previous section. A pole with the length and the diameter of 3L and 0.1L, respectively, is located at the origin of the xy-plane. The   pole is fixed at the bottom of the domain while the top end of the pole is modeled by a hemisphere with the same diameter.
Dirichlet boundary conditions of u = U and v = w = 0 are applied at the inlet, while the no-slip condition (u = v = w = 0) is applied at the bottom. Neumann boundary conditions are applied for all the other boundaries. The Reynolds number is 200 based on the free stream velocity and the length of the flag.
The flapping flag with the pole includes various FSI phenomena such as flow over a flexible splitter plate, flow over a surface mounted cylinder, and interaction of wing-tip, leading-edge, and trailing-edge vortices. In the present research, the problem is employed to examine the capability of the present method for handling both finite-thickness and zero-thickness bodies simultaneously. The flapping flag without a pole is simulated in the same computational domain to identify the effect of the pole. Figure 27 shows time traces of the transverse displacement of the trailing-edge of the flapping flag with and without a pole. T is the oscillation period of the trailingedge at the middle. The Strouhal numbers are found to be 0.216 and 0.232 for the flag with and without a pole,  respectively. The figure shows that amplitudes of displacements with a pole are generally smaller than those of the flag without a pole. In addition, the upper part of the trailing-edge without a pole smoothly changes the direction at the peak while that with a pole abruptly change the direction due to tip vortices shed from the hemisphere top end of the pole. Figure 28(a-b) show time traces of the drag and lift coefficients, respectively. The flag without a pole shows higher drag and lift coefficients than those of the flag with a pole. It is probably because of the direct interaction between the free stream and the flag without a pole while such interaction is obstructed by the pole for the flag with a pole. These phenomena are observed in the instantaneous vorticity contours in Figure 28(c). Figure 29 shows the instantaneous vortical structures around the flag with a pole at four time instants during a period. The combinatory effect of the asymmetric configuration of the flag with a pole as well as gravity makes fluid flow over the flag highly complex.

Concluding remarks
In the present study, a ghost-cell immersed boundary method, which is capable of simulating fluid flow interacting with a moving structure having a finite-or zerothickness, has been presented. To allow large CFL numbers in simulations of flow over zero-thickness bodies, backward time integration for swept-cells and pressure boundary condition with mass conservation are applied, which results in effective reduction of spurious force oscillations without losing sharpness and accuracy.
Composite implicit time integration method is used to obtain a stable solution for the dynamic equation of a thin elastic structure when a large time step size is employed. Interface construction and cell demarcation algorithms to handle both finite-and zero-thickness bodies in a unified manner are devised.
The efficacy of the present method has been throughly examined in several canonical examples. In the simulations of a flapping filament in a free stream and a hovering and flapping flexible wing, it has been found that backward time integration for swept-cells and pressure boundary condition with mass conservation have a small influence on the solution at small CFL numbers. However, numerical experiments show that these treatments become essential for obtaining an accurate solution at increased CFL numbers; without these treatments, the lift coefficient deviates significantly from validated values at large CFL numbers. Furthermore, three dimensional examples including a flapping flag with a pole in a free stream have been simulated to present a unified simulation involving both finite-and zero-thickness bodies at large CFL numbers. The present immersed boundary method is capable of suppressing spurious force oscillations while maintaining the sharpness of an infinitesimally thin body at large CFL numbers. This feature is expected to allow efficient and accurate computation of many FSI problems involving thin structures.
Despite these advantages, the present study has limitations that must be improved in the future research. The present method is vulnerable to instabilities caused by low solid-to-fluid density ratios (strong added mass effects) which are encountered in many biological systems. Many methodologies have since been proposed to overcome such instabilities, but the behavior of these methodologies for thin bodies moving at large CFL numbers has not been investigated. In the future research, the effect of the present method on the methodologies should be analyzed to conduct simulations of zero-thickness bodies accompanying strong added mass effects at large CFL numbers.