AN AUGMENTED IMMERSED INTERFACE METHOD FOR MOVING STRUCTURES WITH MASS.

We present an augmented immersed interface method for simulating the dynamics of a deformable structure with mass in an incompressible fluid. The fluid is modeled by the Navier-Stokes equations in two dimensions. The acceleration of the structure due to mass is coupled with the flow velocity and the pressure. The surface tension of the structure is assumed to be a constant for simplicity. In our method, we treat the unknown acceleration as the only augmented variable so that the augmented immersed interface method can be applied. We use a modified projection method that can enforce the pressure jump conditions corresponding to the unknown acceleration. The acceleration must match the flow acceleration along the interface. The proposed augmented method is tested against an exact solution with a stationary interface. It shows that the augmented method has a second order of convergence in space. The dynamics of a deformable circular structure with mass is also investigated. It shows that the fluid-structure system has bi-stability: a stationary state for a smaller Reynolds number and an oscillatory state for a larger Reynolds number. The observation agrees with those in the literature.


Introduction
The immersed boundary (IB) method is an effective mathematical method for solving problems involving fluid-structure interactions. Such interactions between a flexible body and a fluid are quite common in nature, such as a flapping flag in the air and deformation of red blood cells in the blood. The IB method was originated by Peskin in 1970s to model the blood flow in the human heart [18,19]. For the past decades it has been successfully applied to many other problems. For an extensive list of applications. see the references [20,25,16,4].
The IB method is known to be only first-order accurate and smear out the solution around an interface. The immersed interface (II) method was developed by LeVeque and Li to address the first order accuracy in the IB method for the problems with sharp interfaces [11,12]. The II method also shares some common feature with the Fictitious Domain (FD) method for fluid-particle interactions by which the flow computation is done on a fixed mesh and the rigid body motion of the particle is enforced through Lagrange multipliers in the equations [1]. We refer the reader to [1,17,2,3] for the FD method and its recent applications.
Most of the time-stepping schemes used in the IB/II method in the literature are explicit where the elastic force is calculated on the known configuration of the structure at each time step. Because the fluid-structure interaction is in nature a stiff problem, an explicit IB/II method suffers a severe restriction on time step size. To overcome the stiffness in the explicit IB/II methods, much effort has been made in the recent years to develop implicit or semi-implicit IB/II methods [4,5,22,6,16,7,15]. The augmented immersed interface method is essentially an implicit method by which the the unknown augmented variables are coupled with the flow velocity and the pressure, and computed implicitly during each time step.
In the recent work in the paper [15], Li and Lai presented an augmented immersed interface method for inextensible massless interfaces by which the surface tension ρ(s, t) is treated as an augmented variable and a controlled tangential force is also introduced. Their augmented method preserves both the area enclosed by the interface and the length of the interface. In the work in [7], the authors proposed an augmented method for the incompressible fluid in irregular domains by which the jump in the normal derivative of the velocity is set as the augmented variable. See the references [7,15,13] for the detailed discussion on the augmented method.
In this paper, we present an augmented immersed interface method for simulating the dynamics of a deformable structure with mass in an incompressible fluid. We should point out that the structure was given mass in the context of the immersed boundary method in the previous work in [16,9,26,10]. The fluid is modeled by the Navier-Stokes equations in two dimensions. The surface tension of the structure is assumed to be a constant for simplicity. However, The acceleration of the structure due to mass is coupled with the flow velocity and the pressure. As an implicit method, we treat the unknown acceleration as the only augmented variable so that the augmented immersed interface method can be applied. Unlike the previous work in [15], a controlled tangential force is not needed in our method. We use a modified projection method used in [15] that can enforce the pressure jump conditions corresponding to the unknown acceleration. The augmented equation requires that the acceleration must match the flow acceleration along the interface. The proposed augmented method will be tested against an exact solution with a stationary interface. The dynamics of a deformable circular structure with mass will also be investigated by the augmented method. This paper is organized as follows. In Section 2, we give mathematical formulation for the augmented II method. In section 3, we propose an augmented scheme and briefly describe the numerical solvers. In Section 4, we present the numerical results for the test problems. In Section 5, conclusions are drawn.

Mathematical formulation for the augmented II method
We consider the incompressible fluid in a bounded domain Ω ⊂ R 2 that contains a moving interface Γ(t). See Figure 1 for a diagram of the geometry for the interface problem. The unit normal vector and tangential vector are denoted by n and τ locally on the interface Γ; the angle between the normal vector and the the x-axis is denoted by θ.
The fluid is modeled by the Navier-Stokes equations (1) with the fluid incompressibility constraint (2) where u is the fluid velocity, p is the pressure, ρ f is the fluid density, μ is the fluid viscosity, and F(x, t) is the force term. The interface problem is complete with additional boundary conditions and initial conditions. The singular force F(x, t) on the interface Γ has the following form (3) where σ(s, t) is the surface tension, σ s is the structure density, a(s, t) is the acceleration of the structure, δ(.) is the Dirac function, and X(s, t) is a parametric representation of the moving interface Γ with s being the arc-length. The force term −ρ s a(s, t) is the d'Alembert force due to mass of the structure, as used in the references [10,5]. The acceleration a(s, t) satisfies the following equation (4) Both the surface tension σ(s, t) and the acceleration a(s, t) can be unknown, coupled with the fluid velocity u and the pressure p. As we know, the force term due to σ(s, t) can be rewritten as (5) where κ is the interface curvature. The tension force generally acts in both normal and tangential directions. When σ(s, t) is a constant, the tension force is only in the normal direction, and proportional to the curvature κ.
Using the similar idea of the augmented method in [15], we present an augmented immersed interface method for moving structures with mass by which the acceleration a(s, t) is treated as the augmented variable. In this paper, for simplicity we assume that the surface tension σ(s, t) is a constant. Unlike the previous work in [15], a controlled tangential force is not needed in our augmented method.
With the acceleration a(s, t) and the surface tension σ(s, t), we have the following jump conditions across the interface Γ: In the above, we us the bracket [.] to denote a jump; a n and a τ are the components of the acceleration a in the normal and tangential directions respectively; p n and u n are two directional derivatives in the normal directions.

A modified projection method
In our augmented immersed interface method, the augmented variable is the acceleration a(s, t). The surface tension σ(s, t) is assumed to be a constant for simplicity. Given the jump conditions (6) and (7), the flow can be solved by modifying the numerical method in the paper [15]. The augmented equation is the acceleration constraint (4) on the interface. The novelty of this paper is that mass of the interface is taken into account, and the acceleration is coupled with the velocity and the pressure, and set as the only augmented variable. In the previous work, the interface is massless, and a controlled tangential force is introduced as an augmented variable. Here we use a modified projection method to enforce the pressure jump conditions. Similar to the numerical method in [15], we apply the divergence operator to the momentum equation (1) and reformulate the singular force term by the jump conditions to obtain (8) Therefore, our modified projection method from time t k to t k+1 can be written as (10) Additionally, the acceleration constraint below needs to be satisfied (15) Note that proper boundary conditions for u* are needed on the rectangular domain, depending on a given flow. After the velocity u k+1 has been calculated, we can update the control points of the interface to the new positions and form a new interface by a cubic spline representation.

An augmented approach
We use an augmented approach to solve the equations (10) - (15). Choose the acceleration a as the augmented variable. Let Q be the discrete form for a, and Ũ for u and P. During each time step, given Q the discrete solution Ũ satisfies a linear system (16) where the matrix A corresponds to the prediction step (10)-(11) and the matrix B corresponds to the correction terms in the immersed interface method due to the jump conditions. The residual R(Q) is defined as (17) where Q̂ is interpolated from Ũ by equation (15). The solution Ũ and Q satisfies the linear system (18) The Schur-complement system for Q is (19) Practically, we do not need to explicitly form the matrices A, B, C and D. Assume that the mesh size is N × N, and the total number of the control points on the interface is N 1 .
Compared to the dimension of Ũ (O(3N 2 )), the dimension of Q (O(2N 1 )) is much smaller. We can use a direct method or an iterative method to solve the Schur-complement system for Q. When GMRES is used, the matrix-vector product can be calculated by We can form the coefficient matrix (D −CA −1 B) of the Schur-complement system by setting Q= e l , the lth unity vector, = 1, 2, …, 2N 1 . See the reference [7] for details.

Validation and convergence analysis
We first validate our augmented method by checking the accuracy and convergence of the method against an exact solution in the domain [−1, 1] × [−1, 1]. The exact solution is taken from the reference [14] as follows: where . The interface is the circle . It is easy to check that the velocity satisfies the incompressibility, and it is continuous but has a finite jump in the normal direction across the interface Thus, the normal and tangential forces are Outside the interface, the nonzero and bounded source term is derived from the exact solution.
In our numerical calculations, we take h(t) = 0.20. Note that the exact solution is a steady one. We present the error and convergence analysis by grid refinement at t = 0.5 in Table 1.
The error for velocity is defined as the infinity norm of the difference between the calculated and exact values at t = 0.5. The numerical order of convergence is calculated from the two consecutive errors for a coarse mesh and a twice finer mesh. The mesh sizes are N = 32, 64, 128, and 256 respectively, and the corresponding time steps are 0.05, 0.01, 0.002, and 0.0004. The three calculated orders are 1.95, 2.07 and 1.89 respectively. The calculated average order is 1.97. It indicates that the augmented immersed interface method converges to the exact solution and the order of convergence with repect to space is 2. The plot of the x-component velocity −u(x, y) at t = 0.5 is given in Figure 2.  and near the end of the simulation (below). The interface has the larger deformations during the early stage while it has smaller deformations near the end of the simulation. By the end of the simulation at t = 20, the relative error for area preserving is only 0.37%, which is very small. Both cases in the figure indicate the tapping motion (a stretch and compression) of the interface, as stated in the reference [21]. Figure 4 shows the deviation of a diameter of the circular structure changed over time with different Reynolds numbers. In this test case, the initial interface is a circle of the radius r = 0.3. We chose a diameter (called AB) with two fixed points A and B on the interface. Due to the deformation of the circular interface, the length of the line segment AB changes over time. The deviation of a diameter of the circular structure refers to the difference between the length of the line segment AB and the diameter (2r). In the case Re = 1 (left), the curve shows a stationary state (almost motionless) with no oscillation. In the case Re = 100 (right), the curve clearly shows a stable oscillation after t = 6. The period of the oscillation is approximately equal to 2.46, corresponding to the frequency of 0.4. The fluid-structure system exhibits bi-stability: a stationary state with a smaller Reynolds number and an oscillatory state with a larger Reynolds number. The simulation results are in agreement with the findings in the literature [23,8,21].

Conclusions
In this paper, we develop an augmented immersed interface method for simulating the dynamics of a deformable structure with mass in an incompressible fluid. The acceleration of the structure due to mass is coupled with the flow velocity and the pressure. The surface tension of the structure is assumed to be a constant for simplicity. In our method, we treat the unknown acceleration as the only augmented variable so that the augmented immersed interface method can be applied. We use a modified projection method that can enforce the pressure jump conditions corresponding to the unknown acceleration. The acceleration matches the flow acceleration along the interface. The proposed augmented method has been tested against an exact solution with a stationary interface. The numerical method has a second order of convergence in space. The dynamics of a deformable circular structure with mass has also been investigated. It shows that the fluid-structure system has bi-stability: a stationary state for a smaller Reynolds number and an oscillatory state for a larger Reynolds number. The numerical observation is in good agreements with those findings in the literature. The proposed augmented immersed interface method is also more stable with larger time steps, as it is an implicit method. A diagram of the geometry for the interface problem.  Snapshots of deformations of the interface superposed upon the initial position during the early stage of the simulation (top) and near the end of the simulation (below); the relative error for area preserving is 0.37%. Re = 100. Deviations of a diameter of the circular structure changed over time in a fluid with different Reynolds numbers. It shows two different stable states: a stationary state with a smaller Reynolds number (Re = 1)(left) and an oscillatory state with a larger Reynolds number (Re = 100)(right).