3D Common-Refinement Method for Non-Matching Meshes in Partitioned Variational Fluid-Structure Analysis

# 3D Common-Refinement Method for Non-Matching Meshes in Partitioned Variational Fluid-Structure Analysis

Y. L. Li, Y. Z. Law V. Joshi R. K. Jaiman Department of Mechanical Engineering, National University of Singapore, Singapore 119077
###### Abstract

We present a three-dimensional (3D) common-refinement method for non-matching meshes between discrete non-overlapping subdomains of incompressible fluid and nonlinear hyperelastic structure. The fluid flow is discretized using a stabilized Petrov-Galerkin method, and the large deformation structural formulation relies on a continuous Galerkin finite element method. An arbitrary Lagrangian-Eulerian formulation with a nonlinear iterative force correction (NIFC) coupling is achieved in a staggered partitioned manner by means of fully decoupled implicit procedures for the fluid and solid discretizations. To begin, we first investigate the accuracy of common-refinement method (CRM) to satisfy traction equilibrium condition along the fluid-elastic interface with non-matching meshes. We systematically assess the accuracy of CRM against the matching grid solution by varying grid mismatch between the fluid and solid meshes over a cylindrical tubular elastic body. We demonstrate second-order accuracy of CRM through uniform refinements of fluid and solid meshes along the interface. We then extend the error analysis to transient data transfer across non-matching meshes between fluid and solid solvers. We show that the common-refinement discretization across non-matching fluid-structure grids yields accurate transfer of the physical quantities across the fluid-solid interface. We next solve a 3D benchmark problem of a cantilevered hyperelastic plate behind a circular bluff body and verify the accuracy of coupled solutions with respect to the available solution in the literature. By varying the solid interface resolution, we generate various non-matching grid ratios and quantify the accuracy of CRM for the nonlinear structure interacting with a laminar flow. We illustrate that the CRM with the partitioned NIFC treatment is stable for low solid-to-fluid density ratio and non-matching meshes. Finally, we demonstrate the 3D parallel implementation of common-refinement with NIFC scheme for a realistic engineering problem of drilling riser undergoing complex vortex-induced vibration with strong added mass effects.

###### keywords:
3D common-refinement; Non-matching meshes; FSI; Partitioned staggered; Nonlinear iterative force correction.

, , Corresponding author

## 1 Introduction

Many scientific and engineering simulations that involve interaction of multiple physical fields often require an accurate and conservative scheme to transfer the physical data across non-matching discrete meshes. These problems include electromagnetics [1, 2], contact dynamics [3, 4, 5], conjugate heat transfer [6], and fluid-structure interaction (FSI) [7]. In particular, FSI applications generally rely on different mesh requirements for fluid and solid subdomains to capture the interaction physics accurately, which involves multiple scales and complex multi-modal coupled dynamics. Such requirements of non-matching meshes are common in FSI applications spanning from aircraft wings, deep-water drilling riser, mooring lines, tendons and subsea pipelines [8, 9, 10] to blood blow in arteries and various biomechanical problems. A non-conservative data transfer and locally inaccurate interpolation and projection may lead to poor estimation of flow-elastic response and instability prediction, especially when the frequency is close to the natural frequency of the structure. Therefore, high-fidelity coupled fluid-structure simulations require accurate and conservative treatment of interface boundary conditions across non-matching surface meshes.

For body-fitted Eulerian-Lagrangian coupling, two main approaches exist for the numerical modeling of FSI problems, namely monolithic [11, 12, 13, 14, 15] and partitioned [16, 17, 18, 19, 20]. In the monolithic approach, the flow and structure equations are solved together in a fully coupled manner by assembling the coupled equations into a single block [11, 13, 12]. While the monolithic formulations offer good numerical stability for problems involving very strong added mass effects, the schemes lack the advantage of flexibility and modularity of using existing stable fluid or structural solvers [13, 14, 21]. On the other hand, the partitioned schemes solve the fluid and structure equations in a sequential manner over two decomposed subdomains [22, 19], facilitating the coupling of the existing fluid and structural codes with suitable choices of spatial and temporal discretizations. The partitioned schemes can be classified as either strongly-coupled [23, 24, 25] or loosely-coupled [22, 18] and the surface boundary data must be exchanged or transferred through the interface meshes between the fluid and structure fields.

In a typical partitioned-based FSI simulation, surface meshes at the fluid-structure interface are generally non-matching [17, 26]. This means that their connectivity arrangements are different, and their geometric coordinates may not be coincident due to discretization requirements. Such non-matching meshes and associated data transfer problems also exist in other situations [27], such as adaptive meshing and multi-grid considerations. There can be numerous conservative and non-conservative ways to interpolate and project data across non-matching meshes. The data transfer must be numerically accurate and physically conservative in FSI simulations, especially for those that are time-dependent. This is because errors may accumulate over iterations and long time integration, and a scheme that is both accurate and conservative tends to introduce smaller errors and deliver an improved convergence than non-conservative or locally inaccurate approaches. In our work, we focus on the 3D implementation of common-refinement scheme based on the weighted residual or minimization process [28]. There are two questions one needs to answer when dealing with the spatial coupling methods across non-matching meshes [7, 29]: (a) How to interpolate and project tractions in conservative manner across fluid-solid interface with non-matching meshes? (b) How to integrate the traction vector defined on fluid mesh over the interface elements of solid mesh with their respective shape functions? In [7, 29], a detailed survey of point-to-element and common-refinement based scheme was provided for one-dimensional (1D) and two-dimensional (2D) problems with both flat and curved boundaries. It was demonstrated that the point-to-element schemes can lead to significant errors and sensitivity to grid-mismatch due to a violation of regularity of quadrature rule. These errors can be non-convergent during simultaneous grid refinement of fluid and solid input meshes and can impact the local accuracy along the fluid-solid boundary.

A common-refinement overlay mesh is a surface mesh composed of elements that subdivide the elements of both fluid and solid input meshes simultaneously [30, 31], or simply the intersections of the elements of the input meshes. The common-refinement discretization enables accurate integration of functions that depend on the shape functions of the two meshes [7, 29]. Defined as the topological intersection of the source and target meshes, consistency of the integrations is obtained by performing the numerical quadrature over the common-refinement overlay surface. The common-refinement overlay mesh is constructed such that both the source functions and the target functions are continuous in each of its elements, yielding both accuracy and conservation via accurate integration [29, 32]. For these reasons, the common-refinement scheme is important for stable, accurate and conservative computations of fluid-structure interaction. This scheme is deemed desirable, but can be somewhat complex to implement in three dimensions as compared to other simple nearest neighbor or point-to-element projection methods. Therefore, it has been avoided by application scientists and engineers in their coupled partitioned multiphysics analysis. Indeed it is challenging to compute a common-refinement surface, since the geometrical realizations of the meshes are defined by distinct surfaces with arbitrary mesh intersections. Detailed computational geometry issues related to the construction of common-refinement of two three-dimensional discrete surfaces with curvature and sharp features are provided in [31, 33].

The objectives of this paper are two folds. The first is to quantify the error introduced by 3D common-refinement and compare against the matching reference counterpart. The scope of the present work is to remove the necessity of matching meshes and to redesign a more general projection scheme for non-matching fluid and solid nodes along curved three-dimensional surface. Earlier investigations in [29, 7] were performed in one- and two-dimensional configurations and the fluid flow was considered to be compressible. As shown in [34, 35] for a model elastic plate, the added mass of a compressible flow system has a dependency on the length of time interval, whereas the added mass of an incompressible system asymptotically approaches a constant value as the length of the time interval goes to zero. This fundamental difference in the behavior of compressible and incompressible flows has an implication on the design of partitioned staggered algorithms [25, 36, 37] and therefore it is worth investigating the stability of common-refinement scheme for an incompressible flow interacting with an elastic structure. Furthermore, in the earlier works [29, 7], the common-refinement method was implemented for finite volume fluid and finite element solid solvers. In the present contribution, we employ the common-refinement interface between two consistent Galerkin-based variational formulations for fluid and solid subdomains. We assess the accuracy and convergence of 3D common-refinement for a circular cylinder tube problem in both static and transient situations.

The second objective of this work is to extend our NIFC implementation [25] with the common-refinement scheme for large-scale problems in FSI simulation. The proposed computational framework integrates an ALE-based filtered Navier-Stokes solver, an implicit nonlinear hyperelastic structure solver, and the common-refinement scheme with nonlinear iterative force correction [25, 38]. While the nonlinear structure model is discretized using a continuous Galerkin (CG) finite element discretization, a fluid solver using Petrov-Galerkin finite element spatial discretization and semi-discrete time stepping has been considered for the incompressible fluid flow. The temporal discretization of both the fluid and the structural equations is embedded in the generalized- framework by employing the classical Newmark approximations in time [25]. Owing to domain decomposition strategy in the partitioned iterative procedure, we independently construct the three-dimensional meshes for the fluid and the solid subdomains. The forces from the fluid are applied to the structural boundary as surface tractions, and the structure displacements give a deformation of the fluid subdomain. The fields are advanced explicitly and the interface force correction is constructed at the end of each fluid sub-iteration. During the nonlinear sequence transformation, approximate interface force corrections are dynamically formed through sub-iterations to satisfy the force equilibrium while maintaining the velocity continuity condition along the fluid-solid interface. This iterative sequence coupling relies on the generalized Aitken’s iterated process and the dynamic sequence parameter, which provides a fluid-structure stability at low structure-to-fluid mass ratio [25]. We demonstrate the applicability of the new variational formulation based on the hybrid CRM-NIFC technique for laminar and turbulent flows and compare against the reference solutions. Finally, we demonstrate the 3D non-matching FSI computational framework for the vortex-induced vibration (VIV) prediction of long flexible cylinder in a viscous incompressible flow with strong added mass effects.

The organization of this manuscript is as follows. Section 2 summarizes the flow and structural governing equations with interface coupling conditions. Section 3 gives the spatial discretization of the governing equations and the interface coupling conditions. Discretization details of common-refinement scheme and NIFC-based coupling are described in Section 4. A systematic study on the spatial accuracy of the common-refinement scheme is presented in Section 5, which is followed by the demonstration of the accuracy of common-refinement scheme for the FSI benchmark of cylinder-foil problem in Section 6. Section 7 gives a realistic engineering application of long flexible riser using the coupled framework based on the common-refinement and NIFC schemes. Finally, the work is concluded with some key findings in Section 8.

## 2 Governing fluid-structure equations

Before the presentation of 3D common-refinement scheme, we provide for completeness a brief review of the fluid-structure system. The governing equations for the fluid are applied in an arbitrary Lagrangian Eulerian form while the dynamical structural equation is formulated in a Lagrangian way, and the interface conditions are enforced between the two physical fields.

### 2.1 Incompressible Navier-Stokes with ALE formulation

To simulate the interaction of incompressible fluid flow with a flexible structure, the body-fitted moving boundary based approach is considered in this study. Let be a fluid subdomain at time , where is the space dimension. The motion of an incompressible viscous fluid in is governed by the Navier-Stokes equations given by

 =∇⋅¯¯¯σf+∇⋅σsgs+bf, on Ωf(t), (1) ∇⋅¯\boldmathuf =0, on Ωf(t), (2)

where and represent the fluid and mesh velocities respectively, defined for each spatial fluid point , is the density of the fluid and is the body force applied on the fluid and represents the extra stress term due to the subgrid filtering procedure for large eddy simulation. Here, is the Cauchy stress tensor for a Newtonian fluid, written as , where represents the filtered fluid pressure, is the dynamic viscosity of the fluid. The spatial and temporal coordinates are represented by and , respectively. The first term in Eq. (1) represents the partial derivative of with respect to time while the ALE referential coordinate is kept fixed. Based on the formulation in [38], the filtered Navier-Stokes equations (1-2) in the weak form can be written as

 ∫Ωf(t)ρf(∂t¯\boldmathuf+(¯\boldmathuf−\boldmathwf)⋅∇¯% \boldmathuf) \vspace−2cm=∫Ωf(t) \boldmathbf⋅\boldmathϕf(\boldmathx)d\boldmathΩ+∫Γfh(t)\boldmathhf⋅\boldmathϕf(\boldmathx)d% \boldmathΓ, (3) ∫Ωf(t)\boldmath∇¯% \boldmathufq(\boldmathx)d\boldmathΩ=0. (4)

Here denotes the partial time derivative operator , and are the test functions for the fluid velocity and pressure, respectively. represents the non-interface Neumann boundary along which , where is the normal to the fluid boundary. The update of the deformable fluid subdomain is performed by means of the ALE formulation [39, 40]. The movement of the internal finite element nodes is achieved by solving a continuum hyperelastic model for the fluid mesh such that the mesh quality does not deteriorate as the displacement of the body increases.

### 2.2 Nonlinear hyperelastic structure

We present the principle of virtual work to express the equations of motion and equilibrium of stresses acting on the structure. The principle of virtual work forms the basis for the finite element method for the dynamics of solids, which will be discussed later in the next section. For a dynamically deforming structure with large strains, we use a nonlinear hyperelastic formulation [41] in the coupled fluid-structure system. Consider a solid with mass density that undergoes deformation under external load by fluid flow. Each point on the solid is specified by its position vector. Let denote the initial reference position of a point in an undeformed solid, while denote the displacement of the point in the deformed solid after some time . The function is thus a mapping from initial position to position at time , which completely specifies the change in shape of the solid. The velocity field , which is defined as , describes the motion of the solid under the deformation. The external load, is applied on part of the boundary of the solid. We use to denote the solid boundary that is subjected to external force (Neumann boundary condition) at the reference configuration , to denote the rest of the boundary (Dirichlet boundary condition). Both and form the boundary , such that . By using the principle of virtual work [41], the following weak form is obtained:

 ∫Ωs0τsijδLsijdΩ−∫Ωs0ρsbsiδvsidΩ+∫Ωs0ρs∂usi∂tδvsidΩ−∫Γs2tsiδvsiηsdΓ=0 (5)

Here, is the Kirchhoff stress, is the Jacobian of deformation gradient tensor ; is the Cauchy stress; is the virtual velocity gradients, which satisfies along boundary ; is the body force per unit mass, is virtual velocity field; is an inverse surface Jacobian which relates the boundary surface of and as , where is the isoparametric coordinate. Note that the usual summation on are considered in Eq. (5). The Cauchy stress is related to the left Cauchy-Green stress via the neo-Hookean constitutive law as

 σsij=μs(Js)5/3(Bsij−13Bskkδij)+Ks(Js−1)δij,withBsij=FsikFsjk (6)

where is the Cauchy-Green tensor, and are the shear modulus and the bulk modulus of the solid respectively, and the deformation gradient tensor corresponding to a given displacement field is given as

 Fsij=δij+∂dsi∂xj. (7)

This completes the description of the hyperelastic structure used in our fluid-structure formulation.

### 2.3 Coupling interface boundary conditions

Here we present a short description of the coupling interface conditions for the 3D FSI problem which consists of a fluid domain , a solid domain , and a common interface boundary . For simplicity we only consider the external load through fluid flow as the solid Neumann boundary, i.e. . Two interface boundary conditions corresponding to the continuity of tractions and velocities must be satisfied along . Let be the fluid-solid interface at and be the interface at time . The required conditions to be satisfied are as follows:

 ¯uf(\boldmathφs(xs,t),t) =us(xs,t), (8) ∫\boldmathφs(γ,t)σf(xf,t)⋅ndΓ(xf)+∫γ\boldmathtsdΓ =0 (9)

where denotes the position vector mapping the initial position of the flexible body to its position at time , is the fluid traction acting on the body, and is the structural velocity at time given by . Here, is the outer normal to the fluid-body interface, is any part of the interface in the reference configuration, denotes the differential surface area and is the corresponding fluid part at time . The above conditions are satisfied such that the fluid velocity is exactly equal to the velocity of the body along the interface. The motion of the flexible body is governed by the fluid forces which includes the integration of pressure and shear stress effects on the body surface.

## 3 Partitioned variational fluid-structure system

For the sake of completeness, we next present the discretization using a stabilized variational procedure with equal order interpolations for velocity and pressure. The coupled equations are presented in the semi-discrete variational form for the turbulent fluid flow interacting with a large deformation hyperelastic solid. For a partitioned treatment of the coupled fluid-structure interaction problems, the coupled system is independently discretized with the aid of suitable and desired types of formulations for fluid and structural subdomains, the interface conditions associated with the force equilibrium and no-slip conditions.

### 3.1 Petrov-Galerkin finite element for fluid flow

By means of finite element method, the fluid spatial domain is discretized into several non-overlapping finite elements , , where is the total number of elements. In this paper we adopt a generalized- method to integrate in time , which can be unconditionally stable as well as second-order accurate for linear problems simultaneously. The scheme enables user-controlled high frequency damping, which is desirable and useful for a coarser discretization in space and time. This scheme is implemented by means of a single parameter called the spectral radius which is able to dampen the spurious high-frequency responses and retain the second-order accuracy. With the aid of the generalized- parameters , the expressions employed in the variational form for the flow equation are given as [42]:

 ¯uf,n+1h =¯uf,nh+Δt∂t¯uf,nh+γfΔt(∂t¯uf,n+1h−∂t¯uf,nh), (10) ¯uf,n+αfh =¯uf,nh+αf(¯uf,n+1h−¯uf,nh), (11) ∂t¯uf,n+αfmh =∂t¯uf,nh+αfm(∂t¯uf,n+1h−∂t¯uf,nh) (12)

where

 αfm=12(3−ρf∞1+ρf∞),αf=11+ρf∞,γf =12+αfm−αf. (13)

Let the space of the trial solutions be denoted by and the space of test functions be . The variational form of the flow equations can be written as: find such that :

 ∫Ωfh(tn+1)ρf(∂t¯\boldmathuf,n+αfmh+(¯\boldmathuf,n+αfh−wf,n+αfh)⋅\boldmath∇¯\boldmathu% f,n+αfh)⋅\boldmathϕfd\boldmathΩ + ∫Ωfh(tn+1)(¯\boldmathσf,n+αfh+\boldmathσsgs,n+αfh):\boldmath∇\boldmathϕfd% \boldmathΩ−∫Ωfh(tn+1)\boldmath∇⋅¯\boldmathuf,n+αfhqd\boldmathΩ + nel∑e=1∫Ωeτfm(ρf(¯\boldmathuf,n+αfh−wf,n+αfh)⋅\boldmath∇\boldmath% ϕf+\boldmath∇q)⋅(ρf∂t¯\boldmathuf,n+αfmh +ρf(¯\boldmathuf,n+αfh−wf,n+αfh)⋅\boldmath∇¯\boldmathu% f,n+αfh−\boldmath∇⋅¯\boldmathσf,n+αfh−\boldmath∇⋅\boldmathσsgs,n+αfh−\boldmathbf(tn+αf))dΩe + nel∑e=1∫Ωe\boldmath∇⋅\boldmathϕfτfcρf\boldmath∇⋅¯\boldmathuf,n+αfhdΩe = ∫Ωfh(tn+1)bf(tn+αf)⋅\boldmathϕfd\boldmathΩ+∫Γfh(tn+1)hf⋅\boldmathϕfd\boldmathΓ, (14)

where the lines 3, 4 and 5 represent the stabilization terms applied on each element locally. The remaining terms and the right-hand side constitute the Galerkin terms. The stabilization parameters and appearing in the element level integrals are the least-squares metrics, which are added to the fully discretized formulation [43]. The least-squares metric for the momentum equation is defined as:

 τfm=⎡⎣(2ρfΔt)2+(ρf)2(¯\boldmathufh−wfh)⋅G(¯\boldmathufh−wfh)+CI(μf+μt)2G:G⎤⎦−12, (15)

where is the constant coming from the element-wise inverse estimate, is the size of element contravariant metric tensor and is the turbulence viscosity. The contravariant metric and the least-squares metric are defined as:

 G=∂ξT∂\boldmathxf∂ξ∂\boldmathxf,τfc=1(trG)τfm, (16)

where and are the space coordinate and its parametric counterpart respectively, denotes the trace of the contravariant metric tensor. The element metric tensor intrinsically deals with different element topology for different mesh discretizations.

Linearization of the variational form is carried out by Newton-Raphson technique. Let and denote the increment in the velocity and pressure variables. The linearized matrix form of Eq. (14) is written as:

 \boldmathMfΔ¯\boldmathuf+θf\boldmathKfΔ¯\boldmathuf+θf\boldmathNfΔ¯\boldmathuf+θ\boldmathGΔ¯pf =\boldmathRm (17) −(\boldmathGfM)TΔ¯% \boldmathuf−(\boldmathGfK)TΔ¯\boldmathuf+\boldmathCfΔ¯pf =Rc (18)

where is the mass matrix, is the diffusion matrix, is the convection matrix, is the pressure gradient operator. , , and are the contribution of mass matrix, stiffness matrix and pressure matrix for the continuity equation respectively. is a scalar, in which is the spectral radius that acts as a parameter to control the high frequency damping [38]. and are the right hand residual vectors in the linearized form for the momentum and continuity equations respectively.

### 3.2 Galerkin finite element for hyperelastic structure

Using the principle of virtual work in Eq. (5), we present the finite element approximation for large structural deformations. Without the loss of generality, we consider a hyperelastic material for the structure [44]. The system of structural equations is solved employing the standard Gakerkin finite element technique by means of isoparametric elements for curved boundaries. For an element with nodes, we denote the coordinates of each node by , where the superscript is an integer ranging from to , the subscript is an integer ranging from 1 to 3. The displacement vector at each nodal point will be denoted by . The displacement field and virtual velocity field at an arbitrary point within the solid is specified by interpolating between the nodal values as

 dsi(\boldmathxs)=n∑a=1Na(%\boldmath$x$s)ds,ai,δvsi(\boldmathxs)=n∑a=1Na(\boldmathxs)δvs,ai. (19)

where denotes the coordinates of an arbitrary solid point in the reference configuration and denotes the shape function. By substituting the appropriate deformation measure, the Kirchhoff stress can be calculated. Note that the Kirchhoff stress depends on the displacements through the deformation gradient. Substituting Eq. (19) into the virtual work equation, we obtain the following variational equation for nonlinear structure:

The volume and surface integrals in the above virtual work equation are taken over the reference configuration. Using the inverse surface Jacobian the deformed configuration is mapped back to the reference configuration. The virtual work Eq. (20) gives a set of nonlinear equations with unknowns due to the geometric terms associated with the finite deformations. Notably is a functional relationship, which relates the Kirchhoff stress through the deformation gradient. This nonlinear virtual work equation is solved by means of Newton-Raphson iteration within each time step. We next briefly summarize the linearization process to construct a system of equations.

Let the correction of the displacement vector be . After linearizing the virtual work equation Eq. (20) with respect to , we have the following system of linear equations:

 \boldmathMsΔ¨\boldmathds+\boldmathKsΔ\boldmathds=% \boldmathRs (21)

where is the finite element mass matrix, is the finite element stiffness matrix, and the net force vector , which consists of the external force and the residual terms. The expressions for the mass matrix and stiffness matrix are given by

 Msab =∫Ωs0ρsNbNadΩ (22) Ksaibk =∫Ωs0∂τsij∂Fskl∂Nb∂xl∂Na∂xmFs,−1mjdΩ−∫Ωs0τsij∂Na∂xmFs,−1mk∂Nb∂xpFs,−1pjdΩ −∫ΓfstsiNa∂ηs∂ds,bkdΓ, (23)

whereas the net force vectors is as follows

 {\boldmathRs}ai =∫Ωs0ρsbsiNadΩ−∫Ωs0ρsNbNa∂2ds,bi∂t2dΩ −∫Ωs0τsij∂Na∂xmFs,−1mjdΩ +∫ΓfstsiNaηsdΓExternal fluid force, \boldmathfs. (24)

Through the inverse surface Jacobian , we can relate the nominal and true traction as on the reference boundary . After the current time step solution is obtained, we integrate these equations with respect to time to update the solution in time. For our fluid-structure problems, we consider implicit Newmark time integration method to handle this dynamical nonlinear system. The stiffness matrix is rebuilt at each iteration within time-stepping. In the spirit of partitioned treatment and domain decomposition, the fluid and solid subdomains are solved iteratively. We next present the coupled fluid-structure matrix formulation and the iterative force correction procedure.

### 3.3 Interface force correction scheme

In this section, we present the coupled matrix form of the variational finite element equations defined in the previous subsection at the semi-discrete level for non-overlapping decomposition of two subdomains of fluid and structure. A variational problem of fluid-structure system discretized by a finite element method gives a coupled linear system of equations with the unknowns of fluid and structure in the form , where is a given right-hand side and is the vector of unknowns for the fluid-structure system. Corresponding to the domain decomposition, the set of degrees of freedom (DOF) is decomposed into interior DOFs for the fluid-structure system and interface DOFs for the Dirichlet-to-Neumann (DtN) map. Using the coupled fluid-structure Eqs. (17), (18), (21) and the Dirichlet-to-Neumann map along the interface, the resultant block decomposition of the linear system can be expressed in the following abstract form:

 ⎡⎢ ⎢ ⎢ ⎢ ⎢⎣\boldmathAss00% \boldmathAfs\boldmathAdsI000\boldmathAdq\boldmathAqq000\boldmathAfq\boldmathAff⎤⎥ ⎥ ⎥ ⎥ ⎥⎦⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩Δ\boldmathdsΔ\boldmathdIΔ\boldmathqfΔ\boldmathfI⎫⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪⎭=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩[]c\boldmathRs\boldmathRDI\boldmathRq\boldmathRNI⎫⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪⎭ (25)

where is structural displacement, denotes the fluid unknown variables, and are the displacement and force along the coupling interface. On the other side, and are the right hand side of the corresponding solid and fluid equations; and are the residual errors representing the imbalances during the enforcement of the Dirichlet (kinematic) condition and Neumann (dynamic) condition between the non-overlapping decomposed fluid and solid subdomains. The left-hand side matrix represents the derivatives of the fluid, solid, interface equations with respect to their state variables. The subscripts denote the interior DOFs for the solid and fluid and represent the interface DOFs for the displacement and force. While the block matrix corresponds to the mass and stiffness matrix of the structure in Eq. (21), corresponds to the coupled fluid velocity and pressure linear system in Eq. (17) and (18). is an extraction matrix which maps the solid displacement to the interface, is a matrix which relates the displacement of the interface to the fluid through ALE. is the computation of the force and its mapping to the interface. is a matrix that gets the solid load vector from the fluid-solid interface force. During the partitioned Dirichlet-to-Neumann coupling, there is no explicit availability of the Jacobian matrix entered in the coupled fluid-structure system Eq.  (25).

By eliminating the off-diagonal term in Eq. (25) via static condensation, we can obtain the following reduced linear system:

 ⎡⎢ ⎢ ⎢ ⎢ ⎢⎣\boldmathAss000\boldmathAdsI000\boldmathAdq\boldmathAqq0000˜\boldmathAff⎤⎥ ⎥ ⎥ ⎥ ⎥⎦⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩Δ\boldmathdsΔ\boldmathdIΔ\boldmathqfΔ\boldmathfI⎫⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪⎭=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩[]c\boldmathRs\boldmathRDI\boldmathRq˜\boldmathRNI⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (26)

where the derivation and the terms and are described in detail in [38]. The idea is to compute the iterative correction for the interface fluid force over the deformed ALE configuration as the structure moves. In the nonlinear iterative force correction, we construct a force correction vector in the following manner to correct the previous force at sub-iteration:

 Δ\boldmathfI=˜\boldmathA−1ff˜\boldmathRNI (27)

Here, the force correction vector at the sub-iteration can be constructed by successive approximation, which essentially estimates the coupled effects along the fluid-solid interface. The off-diagonal terms are not explicitly formed and the scheme instead proceeds in a predictor-corrector format by constructing an iterative interface force correction at each sub-iteration. This iterative force correction relies on the input-output relationship between the displacement and the force transfer at each sub-iteration. When the brute-force sub-iterations lead to severe numerical instabilities during strong added-mass effects, the NIFC-based correction provides a stability to the partitioned FSI coupling [38]. The present force correction scheme can be interpreted as a generalization of Aitken’s extrapolation scheme [45, 46] to provide convergent behavior to the interface force sequence generated through the nonlinear iterations between the fluid and the structure. The geometric extrapolations with the aid of dynamic weighting parameter allow to transform a divergent fixed-point iteration to a stable and convergent iteration [38, 25]. We next present the common-refinement scheme to transfer the fluid traction over the solid boundary across non-matching meshes.

## 4 3D common-refinement scheme

In this section, we will address the central topic of this paper focusing on the spatial coupling between the fluid and the structure for non-matching meshes via common-refinement. The fluid and the structural equations are coupled by the continuity of velocity and traction along the fluid-solid interface. We consider that the fluid-structure boundary is discretized independently by the same polynomial order for both fluid and solid subdomains. During the ALE updates, the deformed fluid subdomain using isoparametric elements follows the deformed structure and the discrete fluid and solid meshes remain coincident to the physical fluid-structure boundary.

### 4.1 Variational interface condition

To satisfy the traction equilibrium condition, the momentum flux from the fluid flow must be transferred to the structural surface through the surface traction. To formulate the load transfer operator for the common-refinement method, the weighted residual based on minimization is considered. Let and denote the standard finite element shape functions associated with node of the fluid and node of the solid interface meshes respectively, while and denote the approximate nominal tractions at the corresponding nodes of the discrete fluid interface and solid interface respectively. The continuum traction fields and over and are interpolated as follows:

 tf(\boldmathxf)≈mf∑i=1Nfi~tfi,ts(\boldmathxs)≈ms∑j=1Nsj~tsj. (28)

where and are the number of fluid and solid nodes on the fluid and solid interface meshes respectively. Once we have , , and , we can obtain the transferred distributed loads by solving for . We can measure the residual , by minimizing the norm employing the Galerkin weighted residual method. Multiplying both sides with a set of weighting functions , and integrating over the interface boundary , we obtain:

 ∫ΓfsNsitsdΓ=∫ΓfsNsitfdΓ, (29)

By using finite element approximations, the traction equilibrium condition is:

 ∫ΓfsNsiNsj~tsjdΓ=∫ΓfsNsiNfj~tfjdΓ, (30)

which gives solid-side tractions as

 (31)

where represents the mass matrix for the solid interface elements and is defined by using:

 [Msij]=∫ΓfsNsiNsjdΓ, (32)

and is the nodal force vector given as:

 {fsi}=mf∑j=1~tfj∫ΓfsNfjNsidΓ, (33)

This completes the general formulation of the weighted residual method for extracting the load vector on the solid side interface. While the construction of mass matrix requires only the solid side shape functions, the load vector integral consists of shape functions from both the fluid and the solid sides. For matching grids, this would not cause any problem. However for non-matching grids, the inconsistency of shape functions will lead to integrations across discontinuities. To resolve this issue, there is a need for the common-refinement surface which allows to perform integrations consistently. Further details of common-refinement construction for three dimensional surface meshes can be found in [31].

### 4.2 Algorithmic details

The common-refinement scheme is an important and special data structure, for transferring data between meshes that have different mesh ratios. As shown in Fig. 1, the common-refinement surface of two meshes consists of polygons that subdivide the input boundary meshes of fluid and solid subdomains simultaneously. Every sub-element of a common-refinement mesh has two geometrical realizations, in general, which are different but must be close to each other, to obtain a physically consistent data transfer. In the finite element form, the spatial configuration of the fluid and solid interface meshes can be written as:

 \boldmathxf≈mf∑i=1Nfi(\boldmathx)~\boldmathxfionΓfh,\boldmathxs≈ms∑j=1Nsj(\boldmathx)~% \boldmathxsjonΓsh. (34)

Within this paper, the topology of the common-refinement sub-elements are defined by the intersection of the surface elements of input meshes. These 3D sub-elements are illustrated in Fig. 1. We notice that the intersection of two arbitrary triangles or two hybrid surface elements can be quite complex.

During FSI simulation, within the common-refinement scheme, the load vector over the common-refinement mesh nodes is computed as follows:

 fsj=ec∑i=1∫σciNsj~tfdΓ, (35)

where represents the total number of sub-elements of the common-refinement overlay surface, and represents its th sub-element. We use the Gaussian integration to determine the integration point locations and their weight functions. The basic steps of CRM are summarized as follows:

### 4.3 Common-refinement with NIFC scheme

In this section, we summarize the partitioned iterative coupling of the ALE fluid-turbulence solver with the hyperelastic structure solver, as illustrated in Fig. 2. The solution to the hyperelastic structure equations provides a predictor displacement. The ALE fluid equations with turbulence are then solved to evaluate and correct the forces at the fluid-solid interface. Let the structural displacement be denoted by due to the turbulent fluid forces at time . The first step of the iterative procedure at iteration involves the prediction of the displacement of the hyperelastic structure due to the fluid forces. The computed structural displacements are then transferred to the fluid side by satisfying the ALE compatibility and the velocity continuity conditions at the interface in the second step. This is elaborated as follows: suppose is the mesh displacement at time . The mesh displacements are equated to the structural displacements at the fluid-solid interface to avoid any gaps and overlaps between the non-matching fluid and solid mesh configurations

 \boldmathdm,n+1=\boldmathds   on  Γfs. (36)

The fluid velocity is then equated with the mesh velocity to satisfy the velocity continuity on as

 ¯\boldmathuf,n+αf=\boldmathwf,n+1=\boldmathdm,n+1−\boldmathdm,nΔt    on  Γfs. (37)

In the third step of the iteration , the ALE Navier-Stokes equations with subgrid LES filtering are solved at the mid-point moving mesh configuration to evaluate the fluid forces. The computed forces are finally corrected using the NIFC filter and transferred to the hyperelastic structural solver in the fourth step of the nonlinear iteration. When the solver has achieved the convergence criteria, the fluid-structure solver is advanced in time after updating the variable values at .

In this numerical study, we employ Newton-Raphson technique to minimize the linearization error at each time step and the flow and ALE mesh fields are updated in time by the generalized- method [42]. The resulting incremental velocity, pressure and mesh displacement coming from the finite element discretization are evaluated by solving the linear system of equations via the Generalized Minimal RESidual (GMRES) algorithm proposed in [47], which relies on the Krylov subspace iteration and the modified Gram-Schmidt orthogonalization. Note that to solve the linear matrix system, we do not explicitly form the left hand-side matrix, rather we perform the needed matrix-vector product of each block matrix in pieces for the GMRES algorithm. The solver relies on a hybrid parallelism for the solution of partitioned NIFC-based FSI solver for parallel computing. The parallelization employs a standard master-slave strategy for distributed memory clusters via message passing interface (MPI) based on a domain decomposition strategy. The master process extracts the mesh and generates the partition of the mesh into subgrids via an automatic graph partitioner. Each master process performs the computation for the root subgrid and the remaining subgrids behave as the slaves [48, 49]. While the local element matrices and the local right-hand vectors are evaluated by the slave processes, the resulting system is solved in parallel across different compute nodes [50]. The hybrid or mixed approach provides the benefit of thread-level parallelism of multicore architecture and allows MPI task accessing the full memory of a compute node. After solving the Navier-Stokes equations, we assemble the tractions acting on the hyperelastic structure by means of MPI, and then transfer them onto the surface of the hyperelastic structure through common-refinement overlay surface.

## 5 Error analysis and convergence study

In a typical FSI simulation, the surface of fluid and solid is always intact and coincident without any gaps or overlaps. The fluid load is projected onto the solid surface, while the displacement of the solid is projected onto the fluid surface. Such data transfer is repeated over each time step. There are two types of errors that we are interested during the transferring of data. First, the error in a single transfer of data from one surface to another surface with non-matching mesh. Second, the error in a repeated transfer of data from one surface to another. The analysis for these two types of errors is conducted and described in the following sections.

### 5.1 Static data transfer

The error introduced in a single transfer of data from one surface to another surface with non-matching mesh is investigated here. To quantify this error, we set up two intact surfaces with different mesh size. For the consistency of the notation for FSI, we define one of the surface as fluid, while the other surface as solid. Both surfaces resemble the geometry of a cylinder surface with diameter and height , as shown in Fig. (a)a. Both surfaces are discretized uniformly into elements in -direction. Without the loss of generality, we choose in our analysis. Then, it is discretized into and elements along the circumference of the cylinder. The resultant number of elements on the fluid and the solid surface is thus and respectively. The distributions of element sizes are uniform within each surface. Fig. (b)b shows a typical patch of the non-matching fluid and solid surface meshes.

We define the size of the element on each surface through the area of each element as follows:

 As=πDHwzNs,Af=πDHwzNf. (38)

The degree of non-matching is quantified by using the mesh ratio between these two surfaces, . Several meshes with different mesh sizes are generated by setting