An immersed transitional interface finite element method for fluid interacting with rigid/deformable solid

An immersed transitional interface finite element method (ITI-FEM) is proposed to simulate fluid-structure interaction (FSI). In the framework of finite element method (FEM), the Navier-Stokes equations and the dynamic equation for solid are integrated using the Galerkin method, and the velocity and traction of the fluid are interpolated with those of the solid using the finite element interpolation function. Since the immersed finite element method (IFEM) generates an unphysical velocity/pressure field within the overlapping fluid domain that leads to possible accumulative errors and difficulties in convergence, a ghost fluid domain is introduced to replace the unphysical domain so that the unphysical fluid velocity and pressure are not involved in the equations of the FSI system. A transitional interface is then established to smooth the oscillating solution, along with a momentum-forcing term incorporated into its N-S equations to compensate for the induced errors. Without the inference from the unphysical velocity/pressure, the ITI-FEM has good robustness and accuracy. To validate the proposed ITI-IFEM, a flow over a stationary solid at two different Reynolds numbers and a flow over a moving rigid solid are simulated, examples of the flow-induced interaction with small deformation and finite deformation are also given. The calculated results generally agree with the published results. The proposed method exhibits good capabilities in bio-mechanical engineering application.


Introduction
Fluid-structure interaction (FSI) problems span multiple engineering fields and have been researched for several decades (Akbarian et al., 2018;Faizollahzadeh Ardabili et al., 2018;Zienkiewicz & Newton, 1969). In problems particularly related to aerospace engineering and marine engineering, the computational domain boundaries that usually have complex and time-dependent geometric features may suffer from difficulties in mesh generation and implementation of adaptive solution algorithms. In particular, when a body-fitted method is applied and the solid is moving fast or deforming significantly, the mesh near the solid boundaries may become highly distorted. Although the degradation in mesh quality can be mitigated by using adaptive remeshing algorithms, those algorithms are complex and time-consuming, especially for large scale computations. To circumvent the mesh generation and continuous remeshing, non-body-fitted approaches were subsequently proposed, including the immersed boundary method (IBM) (Peskin, 1977), fictitious domain method (FDM) (Glowinski, Pan, & Hesla, CONTACT Lisheng Liu liulish@whut.edu.cn 1999), and ghost fluid method (GFM). In these methods, the physical domain is embedded into a larger Eulerian fluid domain. The Navier-Stokes equations are solved on the Cartesian fluid grid. The presence of the structure is enforced by adding a continuous/discrete forcing to the N-S equations. The velocities and tractions of the fluid and solid in these non-body-fitted methods can be transferred efficiently with a fast solver, whereas they also have several major disadvantages: the spurious pressure oscillation in IBM, modeling problem of the structure in FDM, and limited capability of simulating the FSI with a moving solid in GFM. Based on this understanding, many novel immersed approaches were explored. One of them is the immersed interface method (IIM) (Lee & Leveque, 2003), which incorporates a part of the continuous/discrete force in the IBM into the jump conditions so that the accuracy of the pressure computation can be improved. Another immersed approach is the sharp-interface IBM, which employs a ghost-cell methodology to satisfy the boundary conditions in the IB for handling a moving/deformable body (Mittal et al., 2008). A current immersed approach, which is improved and extended in the finite element (FE) discretization, is known as the IFEM (Zhang, Gerstenberger, Wang, & Liu, 2004). In the IFEM, the N-S equations are solved using the continuous fluid mesh and the solid stress is calculated with a continuum constitution. Using the reproducing kernel particle function (RKPF) (Hu, 2009), the velocity and coupling force of the solid are discretely interpolated with those of the fictitious fluid within the overlapping fluid domain (Zhang & Gay, 2007) so as to realize the interaction. Nevertheless, the fictitious fluid domain in the IFEM does not physically exist and hence, errors in the velocity/pressure field of this unphysical fluid are inevitably introduced to those of the real fluid and the solid displacement field. When the accumulation of such errors reach a certain threshold value, mismatch will result. To reduce such errors, Wang and Zhang (Wang & Zhang, 2013) have proposed a volume correction algorithm. Nevertheless, the correction still requires iterations and the solution depends on the precision of the interpolation. Considering the interpolation is adversely affected by the unphysical velocity and pressure, thus it is difficult for the corrected solution to converge. Furthermore, the computational cost on the correction scheme becomes uncertain. Therefore, the accuracy and efficiency of the IFEM still present significant challenges.
Based on these considerations, an ITI-FEM was developed in this study, by eliminating the unphysical velocity/pressure field in the overlapping fluid domain and replacing the fictitious fluid with the ghost fluid (Zhang et al., 2004). A transitional interface is then established (Ghias, Mittal, & Dong, 2007) to smooth the physical jump between the ghost fluid domain and the real fluid domain. In addition to introducing a momentum forcing term, the errors induced in the transitional interface are compensated. The proposed method comprises three key component techniques: a signed distance to define the fluid in the transitional interface, a predictorcorrector algorithm to conduct the momentum forcing, and a finite element interpolation function to transmit the FSI messages.
There are several noteworthy features of the proposed ITI-FEM for the FSI problems involving rigid or deformable solids. Firstly, the overlapping fluid domain almost never participates in the computations and the ITI-FEM only couples the interface embedded in the Cartesian fluid mesh with the solid boundary. The computational cost is reduced due to the elimination of the DOFs within the unphysical domain in the numerical computations. Secondly, the transitional interface can stabilize the liquid solution, and smooth the physical jumps between the ghost fluid and the real fluid in space and time. Thirdly, the ITI-FEM avoids the interference of the unphysical fluid in the N-S equations, and the accuracy of the equation system is only related to the numerical computations in the physical domain. The errors induced by the assumed transitional interface can be reduced with the momentum forcing. With these techniques, the ITI-FEM method shows good capability in the FSI problems with both rigid and deformable solids.
This paper is organized into the following sections. In Section 2, the governing equations and the technique ingredients are presented. The FE formulation and the implementation of the techniques are given in Section 3. The efficiency of the ITI-FEM is analyzed in Section 4. In Section 5, several two-dimensional cases are presented to validate the proposed method. The simulation results are compared with the related experimental and numerical results reported in the literature. The paper ends with a brief conclusion in Section 6.

Governing Equation
We consider the model problem of an incompressible Newtonian fluid spanning in a background domain , upon which there is an immersed solid in S . The physical field of the fluid is represented by the Eulerian coordinates x ∈ : ⊂ R d , ∀ d ∈ [2, 3]. The position of the solid in the reference configuration is represented by the Lagrangian coordinates 0 X S ∈ 0 S ; it can also be described by the Eulerian coordinates x( t X S ) at arbitrary time t ∈ [0, T] ( Figure 1).
The Navier-Stokes equations are written as where ρ is the fluid density, v is the velocity vector, σ is the Cauchy stress tensor, g is the body force vector.
Symbol ∇ x represents the gradient operator based on the Eulerian coordinate x, while the operator (·) ,t represents the time derivative of an arbitrary variable (·).
Considering the turbulent effect in high Reynolds number flow, the turbulent stress is included in the Cauchy stress σ as (2) where p is the pressure, I is the second-order identity tensor, μ is the dynamic viscosity, μ t is the Smagorinsky sub-grid scale (SGS) eddy viscosity (Smagorinsky, 1963)  in large eddy simulation (LES), and μ t is expressed as where l s is the mixing length for SGSs, C s is the Smagorinsky constant which is 0.1 in this study (Akin, Tezduyar, Ungor, & Mittal, 2003), and s is the filter width equating to the mesh size (Lilly & Lilly, 1992). The general equation that represents the dynamic behavior of the solid including the structural vibration is written as where σ s is the Cauchy stress tensor of the solid, ρ s is the density, v s is the velocity vector, f s is the body force vector, s is the constant in the Rayleigh damping model (Chowdhury & Dasgupta, 2003) that is used to describe the damping effect of the solid vibration in the FSI system (Yang, Wang, Krane, & Zhang, 2016). In small deformation, the Cauchy stress of the solid is expressed as where D is the fourth-order tensor of material moduli expressed in the initial configuration, ε is the Cauchy strain tensor and ε = ∂d s /∂ 0 X S with d s denoting the displacement vector of the solid. For finite deformation, the solid stress is associated with the initial configuration and where F is the deformation gradient tensor, F = x( t X S )/ 0 X S = I + ∂d s /∂ 0 X S and J is its determinate; S is the second Piola-Kirchhoff stress tensor, relating to the strain energy function W and the right Cauchy-Green deformation tensor C (Bonet & Wood, 1997) as The hyper-elastic material is considered for the case of finite deformation, and the related mathematical derivation are specified in Appendix A.

Technique ingredients
The subdomains and the liquid materials within are different due to the employment of a ghost fluid (Fedkiw, Aslam, Merriman, & Osher, 1999) and a transitional interface (Wang, Wang, & Zhang, 2012). As shown in Figure 2(a), the subdomain F is assumed to be filled with the real fluid, and the part of the volume that overlaps with S is the ghost fluid domain VF . Since VF does not physically exist, the fluid velocity and pressure fields within it almost never participate in the computations of the N-S equations and the interpolations between the fluid and solid. Between VF and F , there is an immersed sharp interface ISI , where a physical jump occurs. Owing to the no-slip condition between F and S , the condition ISI = S should always be satisfied. Let us represent the neighboring region of ISI as MIX , then MIX+ is a part of VF and MIX− is a part of F . To smooth the physical jump, the so-called transitional interface TF is established with a definite width to replace MIX (Figure 2(b)), and this region is filled with a transitional fluid. Here MIX is uncertain and its dimension depends on TF .
Based on the GFM (Fedkiw et al., 1999) principle, nodes in the ghost fluid domain VF are classified into 'ghost nodes' and 'inactive nodes', as shown in Figure 3(a). The ghost nodes are the nodes near S and are assumed to have ghost masses and moments, which are infinitesimal. Their physical parameters are specified as infinitesimal values of the order 10 −19 . The inactive nodes are the interior nodes within S , whose velocity  and pressure are never calculated in the N-S equation calculations and interpolations.
As illustrated in Figure 3(b), the position of the physical jump also changes with time. At any time t, some fluid nodes near S would become ghost nodes at t + t, and vice versa. Such a sharp interface usually leads to the geometrical non-conservation and spurious pressure oscillation (Seo & Mittal, 2011). Thus, the transitional interface aims to alleviate this effect.
The transitional interface is approximately symmetrical to S , and has linearly distributed physical properties. As illustrated in Figure 4, the density ρ and the dynamic viscosity μ within the transitional interface TF are defined as where ρ f and μ f denote the physical parameters of the real fluid while ρ vf and μ vf are the specified infinitesimal values, λ is the transitional indicator, d is the half width of the transitional interface; φ(x, t) is a signed distance, which is negative if x is within TF− , and positive if x is within TF+ ; x s,FSI corresponds to the coordinate of the solid node on S , which is closest to the chosen fluid node.
As the mass and momentum of the ghost nodes are infinitesimally small, the motion of the ghost fluid within the ghost fluid domain can satisfy the N-S equations (Fedkiw et al., 1999). Thus, the velocity and pressure fields of the nodes within MIX can be solved by Equation (1). The mass balance within the transitional interface TF has been proved in Appendix B, and so, the mass equation for MIX is equivalent to that for TF . Due to the property differences between the fluids in TF and MIX , there are discrepancies between their mathematical formulations of the equations of motion. Mathematical compensation should be done for TF .
For MIX , the equation of motion is where the density ρ mix , the velocity v mix and the stress σ mix are different within MIX− and MIX+ and are specified as After the introduction of the transitional interface, with MIX ≡ TF , Equation (10) can be rewritten as Considering that MIX ≡ TF and v mix ≡ v tf based on the no-slip condition imposed at the interface ISI , the equation of motion for TF can be written as where ρ tf is the density of the transitional fluid, v tf is the velocity, σ tf is the stress, and F is a momentum forcing term. According to the Equations (1) and (13), the motion of the transitional fluid can also satisfy the N-S equations. Meantime, the momentum forcing term F has to be additionally incorporated into Equation (1) to compensate for the momentum discrepancy induced by the transitional interface. Then, the velocity and pressure fields within TF become equivalent to the case with the physical jump, which is analogous to the real situation.

General interpolations
The two-way coupling between fluid and solid are achieved through interpolations of physical variables upon ISI and S . The no-slip and no-penetrating conditions are satisfied via the kinematic and dynamic matching between the fluid and solid: (1) Kinematic matching. The velocity message is transmitted from ISI to S as (2) Dynamic matching. The traction message is transmitted from S to ISI as where (σ s · N)| S is the traction imposed upon S and (σ f · n)| ISI is the liquid traction obtained from the immersed sharp interface ISI ; N and n are the unit outward normal vectors of S and ISI , respectively, and their directions are opposite.

Discretized formulation
Let us first assume a suitably defined finite-dimensional trial solution and test function space for the velocity and pressure, denoted by S v , S p , V v and V p = S p (Elias, Coutinho, & Martins, 2006). The finite element formulation of the N-S equations, using the Petrov-Galerkin stabilized technique to meet the inf-sup condition (Hughes, Franca, & Balestra, 1986) can be written as follows: find v ∈ S v and p ∈ S p such that δv ∈ V v and δp ∈ V p : where τ SUPG , τ PSPG and τ LSIC are the stabilization parameters (Takizawa & Tezduyar, 2012); ie is a counter of elements, while ne is the total element number.
In the discretization process, the equal-orderinterpolation velocity-pressure elements are used. The element considered for the purpose is Q1Q1, which has continuous bilinear velocity and pressure (Tezduyar, Mittal, Ray, & Shih, 1992). The velocity v f and pressure p are interpolated along with the test functions δv and δp as where N I denotes the shape function at node I ∈ {1, . . . , nen} and nen is node number per element, c is the test function vector for the velocity vector and q is the test function for the pressure.
Combining Equations (17) and (18), the FE formulation can be obtained. When δv I and δp I are arbitrary, the residual formulation can be written implicitly as with where |J| is the determinant of the Jacobian matrix, W ip is the weight of the Gauss integration at the point ip ∈ [1, nip] and nip is the number of integration points.
In this paper, the Jacobian-free Newton-Krylov (JFNK) method (Knoll & Keyes, 2004) is employed for the nonlinear Equations (19), following the recent studies of the computational fluid dynamics (Chisholm & Zingg, 2009) and IFEMs (Wang et al., 2012;Zhang et al., 2004). In each Newton-Raphson iteration, the residual of the linearized systems of equations are solved using the generalized minimum residual (GMRES) method (Schultz, 1986). Then, the matrix-vector product is computed per GMRES iteration, using a simpler FDM-based calculation. Instead of computing this Jacobian matrix, using the matrix-vector multiplication avoids the need for matrix storage, thereby easing computational requirements, while retaining the Newton-like properties of JFNK method that manifests in a quadratic convergence rate. Additionally, the JFNK is a fully-coupled approach, which does not rely upon fixed point iteration, operator splitting, or loose coupling. For a multiphysics analysis, such an approach could naturally accommodate multiple partial differential equations without requiring specialized elements that couple several unknowns. Based on these premises, the JFNK method with a preconditioning technique was deemed more popular and hence, employed for this study. The preconditioner involved a simple approximation to the sub-blocks of the Jacobian J along the diagonal. For the detailed derivations of the JFNK/GMRES, one can refer to the literature (Sheldon Wang, 2007).
Assuming that d s is the displacement vector of the solid, then v s = d s ,t and the variation is δd s . The weak form of Equation (4) is written as The damping term in Equation (21) can be transformed into the linear combination of the mass and stiffness terms. Using the standard Galerkin method (Smith, Griffiths, & Margetts, 2004), the FE formulation is expressed as with where M s is the mass matrix, C s is the damping matrix, K s is the internal stiffness matrix, and F s is the external force vector; B I is the strain-displacement matrix at node I; f m and f k are scalars, designated as the 'Rayleigh' damping coefficients; D is the matrix of material constants in the current configuration given for the linear elastic problem. Further derivation regarding the nonlinearity of geometric and material can be seen in Appendix C. The time discretization of the equations system (22) is written as where θ is a constant having the value 0.5. This method is known as Newmark's 'β = 1/4' method and is equivalent to the Crank-Nicolson method used in first order problems (Smith et al., 2004). The solution steps for Equation (24) are given as follows: (1) Compute the left-hand side matrix K s = f m + 1 θ t M s + (f k + θ t)K s , and store the matrices using the 'skyline' technique; (2) Factorise K s using the Cholesky factorization to facilitate step (4); (3) Solve the matrix-by-vector multiplications on the right hand side: (4) Solve the linear equations K s d s,t = F by performing the Cholesky forward and back-substitutions on the global matrix; (5) Compute: (6) Increment t and go to (3) for small deformation or (1) for finite deformation.

Implementation
The flowchart of the ITI-FEM is shown in Figure 5. The procedure is similar to that in the literature (Wang & Zhang, 2013), including a solid module, a fluid module, and an FSI module. In the solid module, an independent solid solver is adopted. In the fluid module, the solution procedure consists of a prediction stage and a correction stage, in which F is solved in the first stage and then implemented in the transitional fluid interface in the second stage. In the FSI module, the kinematic and dynamic variables of the fluid and the solid are transmitted via the interface embedded in the fluid and the solid boundary. For all the above modules, the codes are written in FORTRAN 90.

Transitional interface
The quadrilateral elements are used in the computational domain , hence the transitional interface is discretized as stair-step shaped mesh. To obtain the nodes of these interfacial elements, a level-set technique is adopted, and several iso-lines of the level-set function are established, as shown in Figure 6(a).
With the signed distance φ, the transitional indicator values λ of the nodes are calculated and the coordinates  x within each iso-line are determined. The discrete interface FSI,e embedded into the Cartesian mesh is composed of the nodes closest to S , hence its position is then captured, as shown in Figure 6(b). The elements of the real fluid, ghost fluid and transitional-fluid can be further identified (Figure 6(c)). This technique is simple for convex-shaped S,e , and also suitable for the concaveshaped case if an additional operation applied for holes or gaps on S .
When the solid moves, the width of the transitional interface is correlated with t and v s . To smooth the physical jump, the mesh width n of the transitional interface should follow where h e is the mesh size and nh e = d.
It is mentioned in Appendix B that the mass balance of the transitional interface is based on the assumption that the perimeter is far larger than the width 2d. Since ISI = S is always true in order to satisfy the no-slip condition, the perimeter of ISI should equal to that of ISI . Thus, it is evident that L s 2d where L s represents the solid characteristic length. The further stipulation of L s /2d ≥ 10 is reasonable, and the mesh width n should satisfy the relation In Equations (25) and (26), h e and the time step size t are controllable parameters when constructing the FSI model, and never influence the kinematic state of the solid.

Momentum forcing
The momentum forcing is imposed using a predictorcorrector algorithm (Rüberg & Cirak, 2011). The algorithm divides the fluid solution process into the two stages: (1) Stage 1. In the prediction stage, the N-S equations are solved using the second-order fractional-step FEM (Codina, 2001) to obtain F.
Assuming that t+ tv and t+ tp are the solutions at arbitrary time t + t, the corresponding intermediate variables can be expressed as where t+ξ t τ is the deviatoric stress; t+ t v * is the prediction value of velocity. Notation ξ is the intermediate Combining Equations (27) and (1), and the momentum equation is written as Then, the prediction value t+ t v * is obtained. The pressure Poisson equation is written as where η is another intermediate coefficient and t+ t p * is the current predictor value. Consequently, F is calculated by taking t+ t v * and t+ t p * to Equation (14).
(2) Stage 2. In the correction stage, F is brought into the N-S equations to balance the difference induced by the transitional fluid. In this stage, the FSFEM is not adopted. Solving the N-S equations, the values of t+ t v and t+ t p are the final solution.
In the two stages, the prescribed boundary conditions, including the velocity and traction conditions, are constant. The nodes on FSI,e and S,e contribute few to the unknown DOFs (Kim & Moin, 1985). It is worthy to note that reasonable iterations can improve convergence.

Interpolation
The interpolation is described through a representative FSI model, including a fluid element F,e i and a solid element S,e j , as shown in Figure 7. The FE interpolation function (Wang & Zhang, 2010) is used for the interpolation near the sharp fluid-solid interface. The FE basis function satisfies the reproducing conditions and involves no correction function relative to the RKPF. This function is ready-made as the shape function in the Galerkin method. The interpolations are implemented as below: (2) Traction distribution is conducted as

Efficiency analysis
In this section, the efficiency of the proposed ITI-FEM is analysed using the classical complexity theory (Biirgisser, Clausen, & Shokrollahi, 1997) since there have been no open-source codes of the IFEMs. We choose the latest version of the IFEM or the m-IFEM to assess the efficiency. The sums of the computational complexities in the mentioned three modules are compared.
In the solid module, both methods employ an independent solid solver. DOFs in the dynamic equations are the same for the same solid model. Therefore, the computational complexity of the two methods differs little regardless of using Newmark-β algorithm (Smith et al., 2004) or α-method (Hughes, 2008) (Strang, 2006).
In the other two modules, their complexities are given in Table 1. In the fluid module, the ITI-FEM also solves the N-S equations with the combination of the GMRES and JFNK, similarly to the IFEMs (Sheldon Wang, 2007). Comparing the ITI-FEM with the m-IFEM, the evident differences in the fluid module are the DOFs in the N-S equations and the imposition of the momentum forcing. For a 2-D case, a single fluid element has 3 DOF per node. Assuming that N F ,N VF and N TF represent the total node numbers within F ,¯ VF and TF , respectively, and the DOFs are 3(N F −N VF ) for the ITI-FEM and 3N F for the m-IFEM. So, the IFEM has 3N VF of DOFs more than the ITI-FEM. Again, we let n nr denote the Newton-Raphson iteration, and n in and n out represent the inner and outer iterations in the GMRES method. Namely, the efficiency of the ITI-FEM is validated.

Flow over a rigid circular cylinder (Re = 45 and Re = 361)
In this section, the benchmark problems of flow over a rigid solid are solved. Flow characteristics of different Reynolds numbers (Re = ρ f U inlet D/μ f ) are quantitatively compared with the published results to validate the ITI-FEM (Singha & Sinhamahapatra, 2010;Thom, 1933;Wang & Zhang, 2013). The first case is a steady analysis at Re = 45, for which both experimental and numerical results have been obtained. As specifications of the numerical model are more lucid (Singha & Sinhamahapatra, 2010), the geometry of this case is chosen to match it. As shown in Figure 8(a), there is a channel of length L = 35 cm and height H = 8 cm, and the H/D ratio is chosen considering its effect on the pressure (Mou, He, Zhao, & Chau, 2017), which is observed to be the least in the work of Singha and Sinhamahapatra. A circular    Table 2. Firstly, we study the number of GMRES iterations and Newton iterations in the flow problem at Re = 45, as shown in Figure 9(b)−9(d). For the GMRES iterations, both the inner and outer iterations are studied. It shows that the norms of the residuals finally converge to the infinitesimal value (10 −5 ) with the increase in either the inner or outer iterations, which indicates the stability of the solutions. The variations of the Newton-Raphson iterations follow the same trend. Considering the computational cost and accuracy, we examine 30 inner and outer iterations in each Newton-Raphson iteration, and the Newton iterations is set as 10. The computation consumes 0.14 s for only one inner iteration, one outer iteration, and one Newton iteration. It requires a memory storage of 217 MB. And it is implemented using an Intel Core i7-4790HQ processor with a main frequency of 3.60 GHz.
The pressure results for different cases at Re = 45 are presented in Figure 9(a). It shows that our result agrees well with the experimental result (Thom, 1933), whereas the resulting pressure is smaller than that in the given numerical result (Singha & Sinhamahapatra, 2010). The resulting drag coefficient C d = F drag / 1 2 (ρ f v 2 ) of 1.35 in the ITI-FEM is also smaller than the numerical value of 1.43 reported by Singha and Sinhamahapatra. However, both are within the range of [1.33, 1.52] at Re ∈ [40, 60] (Apelt, 1958). Different element types and discretization schemes could explain the discrepancy between the results calculated by the ITI-FEM and the published unstructured collocated grid finite volume method (FVM). In summary, the agreement between our result and previously reported results indicate the sufficient accuracy of the ITI-FEM fluid solver in dealing with flow problems at low Reynolds numbers.
To compare with an immersed-finite-element type method, the second case considered here is the transient analysis of the flow around a circular cylinder at Re = 361 (Wang & Zhang, 2013). Since a 2D model is typically used to represent the 3D case with the planar assumption (Pianet & Arquis, 2008), our model is built to match the planar geometry of the model available in the literature. As shown in Figure 10(a), the aspect ratio of the channel is 4 and the diameter of the cylinder is 0.5 cm. The cylinder is placed 1.5 cm away from the inlet. The mesh is denser at the subdomain of [0, 3] × [0, 2] with the minimum element size of 0.02 cm. The properties and discretisation parameters are given in Table 2. We set the GMRES (outer) iterations as 30 and the Newton-Raphson iterations as 10, for which the solution is convergent (the default value for inner iteration is 30 for all the presented cases). It costs 0.875 s for each GMRES iteration and each Newton iteration, with a memory usage of 103 MB, when using an Intel Core i7-4790HQ processor with a main frequency of 3.60 GHz.
Comparing the resulting C d and C l with those published in the work of (Wang & Zhang, 2013), our results match well with that calculated using m-IFEM. The mean value of C d is 1.066, which is close to the experimental data of [0.994, 1.14] in the vicinity of Re = 361. The resulting C d is also close to 1.15 which is available in the previously published literature (Wang & Zhang, 2013). The mean value of C l is −0.00172, which is approaching zero. Regarding their vortex shedding frequencies, the Strouhal number St = f 0 D/U inlet is calculated as 0.227, where f 0 denotes the vortex shedding frequency equal to 45.4 s −1 . The resulting St is also close to the previously published value 0.205 in the work of Wang and Zhang, with the deviation −10.3%. Figure 11 illustrates the Karman vortex street phenomenon. It shows that the ITI-FEM can capture the process of vortex formation and shedding at the back of the cylinder. This phenomena matches the results obtained using the m-IFEM, especially with regard to the orientation of the tail of the separation vortex. Again, the accuracy of the ITI-FEM in simulating the transient flow at a higher Reynolds number is verified.

Flow over a rigid moving circular disk
Considering that the transitional interface influences the solid motion, the effects of the influence factors on the solid dynamics are analysed in the second example. As shown in Figure 12, this example is a 2-D disk dropping in a fluid channel due to gravity. The diameter D of the disk is 1cm, and the thickness is 0.05 cm. The size of the channel is 15 × 4 cm. The disk is positioned 0.7D away from the top boundary of the channel. The no-slip condition is set upon the sidewalls of the channel. The disk is initially motionless, and its rotation is constrained. The properties and discretization parameters of the computational model can be seen in Table 3. The step size for time is 0.001 s and each time step costs 240 s when  using an Intel Core i7-4790HQ processor with a main frequency of 3.60 GHz. The corresponding memory usage is 413 MB. The numbers of inner and outer iterations are 30 and 10, respectively; it is found that the final residual norm converges to 5 × 10 −5 ∼ 6 × 10 −5 after 10 Newton iterations.
The time-dependent variation of velocity and acceleration are plotted in Figure 13. It is seen that the width of the transitional interface has a more significant effect on the solid acceleration (Figure 13(b)) while the momentum forcing having a more significant effect on the velocity (Figure 13(c)). The width can smooth the oscillating pressure imposed upon the solid. The resulting acceleration, which is the smoothest when d = h e , also suggests that Equations (25) and (26) should be satisfied. When d = h e , the momentum forcing is shown to help the velocity converge to a lower value and has a more apparent effect than the iterations do.
The solid dynamics agree well with the prediction reported in the literature (Lapple & Shepherd, 1940).
The velocity of the disk increases continually until the decreasing acceleration ultimately fluctuates near zero. To be specific, the drag coefficients of the ITI-FEM and the reported are compared quantitatively. When d = h e and the momentum forcing is imposed, the drag coefficient C d is 1.18 at Re = 4620 for more iterations and 1.23 at Re = 4550 for few iterations. Both drag coefficients approximate to the experimental results (Concha & Barrientos, 1986), where C d ∈ [1.18, 1.20] at Re ∈ [2250,9020]. This agreement proves that the incorporation of the momentum forcing into the N-S equations can improve the solution. Moreover, the proposed method can obtain satisfactory results for fluid flow over a moving rigid solid.

Flow over a linear-elastic unilateral leaflet with small/definite deformation
The third example is a study of the aero-elastic response of a unilateral leaflet in a 2-D channel, as shown in Figure 14. This example is presented to validate the ITI-FEM applicable to the FSI problem with deformable solid. The resulting displacement of the unilateral leaflet are compared with those of the IFEM, m-FEM (Wang & Zhang, 2013) and ALE method. To facilitate the comparison, the geometry and material properties of the presented model are the same as that in the work of Wang and Zhang. As shown in Figure 14, the size of the channel is 8 × 2 cm; the leaflet has the size of 0.5 × 1ï£¡ï£¡ï£¡cm, positioned 3 cm away from the inlet.   Table 4.
The discretization details are given in Table 5. The ALE method is implemented by ANSYS Workbench and ABAQUS using the re-meshing technology. Due to the usage limitation of the software, a three-dimensional (3-D) body-fitted model is built for the ALE method. Only two elements are set in the z-direction to reduce the 3-D effect. The mesh size of the solid is different from that in the non-body-fitted model. The time step is 0.001 s for all these computations. Finally, it takes 3.25 h with a memory usage of 295 MB for the ANSYS Workbench, 2.25 h with a memory usage of 1.72 GB for ABAQUS, and 2.7 h with a memory usage of 74.5 MB for the ITI-FEM. Of the three, simulations ABAQUS consumes the least   CPU time whereas it requires the largest memory; compared to other methods, the ITI-FEM shows good efficiency if the specifications of the hardware available are limited. Thus, the efficiency of the method is verified. When using the ITI-FEM, the numbers of inner and outer iterations are 30 and 10, respectively; the number of Newton iteration is 10. Here, all the computations using the software and the ITI-FEM program created in this study are implemented on an Intel Core i7-4790HQ processor with a main frequency of 3.60 GHz.
For case 1, the quantitative comparison of different methods can be seen by measuring the x-coordinate at the top right corner, as shown in Figure 15(a). All timedependent variations of the x-coordinate follow the same trend and reach steady state after 1.4 s. In details, our result is precisely between the results in the IFEM and the m-IFEM, as well as close to those calculated by the ALE method. The no-slip conditions are also satisfied (Figures 15(b) and 15(c)). Therefore, the accuracy of the ITI-FEM can be verified in the FSI problem with small deformation.
In case 2, the rotation of the solid is more apparent than that in case 1 due to the softer material, which also indicates the limitations of applying the Workbench and ABAQUS software to such cases. The resulting configuration of the solid at the steady state is only compared with that of the m-IFEM, as shown in Figure 16(a). It shows that our result agrees with the published result, except for a slight difference in the volumes of the deformed solid. As shown in Figure 16(b), the volume change is 3.2%   Table 6. Physical properties and discretizations for the biomechanical model.

Material
Discretization Properties Blood 44,589 elements ρ f = 1.0 g/cm 3 44,000 nodes μ f = 0.035 g/(cm·s) Quadrilateral Single leaflet 120 elements ρ s = 2.7 g/cm 3 155 nodes µ 01 = 2.0 g/(cm·s 2 ) Quadrilateral µ 10 = 1.0 g/(cm·s 2 ) λ s = 0.0 g/(cm·s 2 ) when the displacement of point P converges to a stable solution. The compressibility is so small that the width of the transitional interface varies little, and the mass balance can be still satisfied. The computational cost and the memory usage in case 2 are the same as that in case 1, since the settings of the GMRES and Newton iterations as well as the hardware and computational scale are the same.
According to the obtained results of the two cases, the displacement in the m-IFEM is more conservative than the values calculated in this study for small deformations. It agrees with the previously mentioned prediction that, coupling of the nodes of the unphysical domain may lead to an additional and excessive constraint in the solid response, namely, the possible accumulative error. In the case of finite formation, the agreement between the m-IFEM and ITI-FEM shows such constraint upon the solid response is not obvious, which may attribute to the iterations. However, it can still be seen that a ghost mesh treatment on the overlapping domain obtains a satisfactory result in the FSI problem with deformable solid. In particular, the ITI-FEM can obtain a better solution for the small deformation.

Flow over hyper-elastic bilateral leaflets with material and geometric nonlinearities
In this section, the ITI-FEM is employed to simulate a bio-mechanical problem of blood-induced valve motion in a vein (Hossler & West, 1988). The bio-mechanical model is simplified as a straight tube with bilateral leaflets in geometry, as shown in Figure 17. The size of the channel is 4.0 × 0.7 cm. The size of the leaflets is 0.04 × 0.3 cm. The leaflet is positioned symmetrically 0.5 cm from the inlet. There is an opening of size 0.1 cm between the leaflet tips. We assume a muscular pressure of 2.4 e + 04 dyne/cm 2 is prescribed at the inlet, along with the uniform velocity of 4.5 cm/s. The no-slip condition is set upon the wall. This simulation lasts for 0.2 s and only considers the opening phase of the valve cycle (Lurie, Kistner, & Eklof, 2002).
In the computation, the material and geometric nonlinearities of the leaflet are considered (Appendix C). The leaflet is defined as an incompressible hyper-elastic material, using the mentioned Mooney-Rivlin model. The physical properties (Martinez, Fierro, Shireman, & Han, 2010) and discretization details of the bio-mechanical model are illustrated in Table 6. The time step size is 0.0001 s and each time step consumes 180 s, with the memory usage of 200 MB. Ten iterations of the GMRES as well as Newton-Raphson iterations are considered. This case is simulated using an Intel Core i7-4790HQ processor with a main frequency of 3.60 GHz. Figure 18 shows the blood flow characteristics and the valve motion at different times. Subjected to a strong pressure gradient between the leading and trailing edges, the bilateral leaflets deform symmetrically like an opening door. The blood flow flux increases at the opening gap and the flow is fast across the orifice. Pairs of symmetrical, shedding Karman vortices are generated alternately from the bi-leaflet tips (Figure 18(a) ∼ 18(b)) since the Reynolds number is increasing, as shown in Figure 19(a). With time marching, the orifice becomes greater (Figure 18(c) ∼ 18(e)) and the deformation is nonlinear, as shown in Figure 19(b). In addition, the hyper-elastic leaflets are always incompressible during the presented period, see Figure 19(c). These calculated results match with the existed description of the opening phase.
The pressure imposed upon one leaflet at the time of 0.2 s is presented in Figure 20(a). The strong pressure gradient near the leaflet is 1.9 kPa. The maximum value of wall shear stress appears at the base region of the leaflets, which is consistent with the position reported in the work of Soifer and Weiss et al. (Soifer, Weiss, Marom, & Einav, 2016). Again, the first principal stress is also calculated 1.45 kPa at the base region (see Figure 20).
In summary, the above simulations are, in general, challenging due to the nonlinearities present. The proposed method provides an acceptable approach to predict the details of the valve dynamics in the opening phase. The simulation illustrates that the proposed method can accommodate flexible-structure-fluid problems with geometric and material nonlinearity.

Conclusions
In this study, a novel method, the ITI-FEM, is proposed based on the FE discretization of the immersed approach. By using the transitional interface, the immersed approach can be used for the resolution of several problems or classes of problems pertinent to physics and mathematics. Firstly, compared to other IFEMs, the ITI-FEM does not generate an unphysical fluid in relation to the N-S equations, and is qualitatively more efficient with a smaller number of DOFs that need to participate in the calculations related to the FSI system. Secondly, as the unphysical fluid velocity and pressure are eliminated, the solution of the equations in the FSI system can be physically accurate, especially with a momentum forcing term into the N-S equations of the transitional interface. Thirdly, the quantitative and qualitative comparisons in the examples indicate that the proposed method has good robustness and accuracy. When the width of the interface is reasonable, the oscillating pressure in the N-S equations can be smoothed and the FSI solution converges easily to the steady state. Using the ITI-FEM, the resulting response of the deformable solid is close to that calculated by the ALE method.
It is notable that this method is still in its preliminary stages, and currently, it can address only 2-D problems. There is still ample room for improvement with respect to several aspects, such as its applicability for compressible fluids, 3-D problems, advanced/novel solvers, and the enrichment of the solid constitutions. Keeping those limitations in mind, our future work is to improve the ITI-FEM theory and broaden its applications to civil engineering problems. A particularly interesting potential application is the turbulent analysis at high Reynolds numbers which involves a 3-D analysis of the flow characteristics, commonly required in aerospace engineering. With respect to the design of the bio-valves, the fatigue failure and fracture may also be considered by utilizing new types of material models in the solid module. Figure A1. A transitional interface with arbitrary shape. and tangent quantities as σ ij ≡ (2µ 01 + 4µ 10 )δ ij + 2µ 01 ij , c ijkl = λ s δ ij δ kl + 8µ 01 I ijkl δ ij δ kl (A6) which agrees with the formulation of the linear elasticity in (Taylor, 2012). The above derivations show that parameters µ 01 and υ can be transformed to the Lame moduli and give the same response in small deformations as a linear elastic material. For finite deformation or J =1, the modified Lame constants can still describe the material nonlinearity when the material moduli is time-varying with the right Cauchy-Green deformation tensor C or the left Cauchy-Green deformation tensor b.