{'fraction_non_alphanumeric': 0.032478285849471086, 'fraction_numerical': 0.01492108553854473, 'mean_word_length': 4.729402652725336, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A strategy to interface isogeometric analysis with Lagrangian finite elements - application to fluid-structure interaction problems

Isogeometrically enriched finite elements offer efficient localized isogeometric analysis (IGA) enrichment for numerical simulations involving large computational domains. This is achieved by employing surface enriched elements to interface isogeometric elements with classical Langrangian finite elements. In this paper, we explore their applicability and merits for fluid-structure interaction (FSI) analysis. The implemented approach not only offers an enrichment of the finite element space, but also offers a framework for discretizing and analyzing fluid and structure with different finite element approaches, namely, classical Lagrange finite elements and IGA. In this context, a monolithic solution approach with an explicit grid update mechanism is implemented for FSI. The applicability and the impact of the isogeometric enrichment approach on the accuracy of the numerical solution is assessed by comparing the obtained results with existing reference solutions of FSI benchmark examples involving two- and three-dimensional incompressible fluid flow past hyper-elastic solids.


Introduction
Fluid-structure interaction (FSI) belongs to a special class of multi-physics where a deformable solid structure interacts with a body of fluid. The interactions between the constituent media occur at the interface, where the two media meet, and typically involve transfer of interfacial forces from one media to the other. These forces then determine the kinematical response of the constituent media. These interactions are typically strong for many physical and engineering applications, and the absence of a coupled analysis will lead to misleading, and possibly catastrophic conclusions. Blood flow in arteries and veins, design of mechanical cardiovascular valves and pumps, flutter and fatigue of aircraft wing and turbine blades under aerodynamic loadings, and safety calculations of suspension bridges and tall buildings under strong winds, are only a few of the many examples that demand a coupled analysis approach that takes into account the cross-interactions between solid and fluid media.
Mathematical models of FSI problems are often too complex to be analyzed analytically. Therefore such problems are largely investigated through experiments and numerical simulations. In the context of numerical analysis, the available modeling approaches for FSI analysis can broadly be categorized into three groups, namely field elimination, partitioned approaches and monolitic approaches. Field elimination methods (e.g. see Belvins [1990], Dütsch et al. [1998], Ohayon [2004]) constitute the simplest of these approaches, where one or more of the associated field variables are eliminated from the governing model through reduction techniques such as substitution, integral transforms or model reductions. Such approaches, however, are restricted to the modeling of linear problems where a decoupling between the field variables can be formulated. On the other hand, available FSI solvers are generally classified as a partitioned (Felippa et al. [2001]) or a monolithic (Bendiksen [1991]) approach, depending upon the mechanism by which the kinematic and the kinetic description is transfered between solid and fluid media at the fluid-solid interface.
Partitioned approaches make use of existing stand-alone fluid and structure solvers by invoking them in a sequential manner. The coupling interactions at the interface are enforced as forcing effects communicated through prediction, substitution and synchronization techniques. Whether these communications are performed iteratively or in a staggered manner, partitioned approaches are further categorized as strongly (Tezduyar et al. [2006], Dettmer and Perić [2006], Matthies et al. [2006], Küttler and Wall [2008]) or weakly coupled schemes (Piperno et al. [1995], Farhat [2004], Wang et al. [2004], van Zuijlen et al. [2007]), respectively. In contrast, monolithic solvers cast the governing equations of the fluid and the solid model into a singular unit and attempt a solution in a single iteration with consistent time integration schemes (Hübner et al. [2004], Walhorn et al. [2005], , Klöppel et al. [2011], Mayr et al. [2015]).
Partitioned and monolithic approaches have been compared by many researchers (e.g. see Rugonyi and Bathe [2001], Heil et al. [2008], Küttler et al. [2010]), however the choice of a preferred approach for a particular application predominantly remains problem and resource specific. The strongest merit of partition approaches is the flexibility to incorporate validated and optimized stand-alone fluid and solid solvers in the coupled analysis framework, irrespective of the spatial discretization methodology. Finite volume method (FVM) is a popular choice for fluid flow problems, while finite element method (FEM) is used extensively in the structural mechanics community. However, the sequential invocation of the solvers leads to time-inconsistent field description, which if left untreated, can lead to numerical instabilities (e.g. the added mass effect Causin et al. [2005], incompressibility dilemma at the interface Küttler et al. [2006], etc.). Although monolithic methods lead to time-consistent solutions and subsequently to a more robust solution approach, they demand larger computational resources due to storage and solution of a larger system of equations. Moreover, due to a combined cast for the governing equations of both models, monolithic schemes tend to lead to ill-conditioned global matrices for cases where material parameters of the constituents are orders of magnitude different from one-another. For such cases, preconditioning strategies (such as Heil [2004], Gee et al. [2011]) should be employed or developed before the monolithic system is solved.
In this paper, we present a monolithic FSI solution procedure in the framework of isogeometric finite elements. Initially proposed as a subset of the FEM in Hughes et al. [2005], isogeometric analysis (IGA) offers an improvement to the classical FEM by incorporating CAD discretizations based on NURBS (Piegel and Tiller [1997]) and T-splines -rather than interpolatory Lagrange polynomials -to the trial and test function spaces of the discrete finite element formulation. In this manner, the original CAD description (and consequently the original geometry) is retained in the discrete finite element setting and a direct correspondence between CAD and finite element analysis is maintained. Although initially formulated to speed-up the design-to-analysis transition, the use of spline-based functions provides additional merits for numerical analysis over classical interpolatory polynomials due to many desirable properties of the spline basis (such as accurate geometry representation, higher continuity across elements and the variation diminishing property). As a results, IGA has lately been the focus of intense technical research with wide-range engineering and biomedical applications. For an elaborate overview, interested readers are invited to consult Cottrell et al. [2009].
From biomedical applications (Zhang et al. [2007], , Chivukula et al. [2014]) to the design of wind turbines (Bazilevs et al. [2011(Bazilevs et al. [ , 2012(Bazilevs et al. [ , 2016), IGA has been successively applied to various FSI problems. The superiority of IGA over interpolatory FEM for uncoupled problems has been investigated by many researchers, e.g. , Akkerman et al. [2008], Motlagh et al. [2013], Rasool et al. [2016] for fluid dynamics and Cottrell et al. [2006], , , Morganti et al. [2015] for structural mechanics applications. For many of these studies, IGA yields higher accuracy per degree-offreedom (dof) than the classical FEM. However, modeling the entire computational domain with spline-based isogeometric elements may not be desirable for every problem. Due to recursive evaluation of the basis functions (which for the case of NURBS and T-splines can be distinct over each element), IGA demands additional computational resources for discretizations with comparable dofs. Moreover, it has also been observed that the extended continuity of the basis across element boundaries incurs additional burden on linear solvers (see Collier et al. [2012Collier et al. [ , 2013). For many problems, especially true for simulation of physical flow in unbounded regions, it is sufficient to have an enriched analysis only in certain localized regions of the discretized domain, while the remaining bulk can be represented with a low-order computationally efficient discretization, e.g. Rasool et al. [2016] and Harmel et al. [2017]. Such applications motivates a framework where an interpolatory finite element discretization is enriched with isogeometric elements only at regions where an enhanced representation and analysis is beneficial.
The motivation to blend isogeometric elements to an interpolatory finite element discretization is certainly not novel. Perhaps, the earliest reference in this regard can be found in Bazilevs et al. [2012], where T-spline shells for solids were interfaced with low-order Lagrangian finite elements for fluid at the fluid-solid interface. Such a treatment was achieved by weakly enforcing the coupling conditions at the interface, which also relieved the necessity of a matching discretization at the interface. However, the interfacing between isogeometric and Lagrange finite elements is only available at the fluid-solid interface. Interpolating B-spline curves with Lagrangian functions, transformation maps were developed in Lu et al. [2013] and Ge et al. [2016] to construct "Blended elements" which can interface two-and three-dimensional NURBS-based elements with Lagrangian elements.
Another strategy developed with a similar motivation is the isogeometric zone enrichment approach. It was proposed in Rasool et al. [2016] for fluid flow problems and in Harmel et al. [2017] for thermal analysis. In this paper, we explore its applicability and benefits for FSI applications. Isogeometrically enriched surface elements (see Sauer [2014, 2015]) are used to interface IGA elements with Lagrange finite elements, yielding the possibility of a locally enriched analysis. The strategy additionally offers a monolithic mechanism, where one medium can be analyzed in the framework of IGA, while the other in the classical finite element setting. The potential of such an analysis is also investigated.
The remainder of this paper is organized in the following manner: The framework of the governing monolithic FSI model is discussed in Section 2. It includes the continuous differential model, its discretized weak form and a grid motion mechanism for incorporating moving surfaces. The isogeometric zone enrichment strategy is explained in Section 3, while results from several numerical experiments are presented in Section 4. Conclusions from the analysis are discussed in Section 5.

2 Monolithic FSI model
The theoretical framework for the developed monolithic FSI solver is discussed in this section. It begins with the description of the governing equations that describe momentum and mass balance for the entire system. The equations are first expressed in strong form. Their variational (weak) description follows consequently. The implemented finite element formulation, in conjunction with the stabilization technique used in this study, is then presented. A note on the explicit grid motion mechanism concludes this section.

Governing equations in monolithic setting
Let us consider disjointed spatial domains B f ⊂ R d and B s ⊂ R d at any given time t ∈ [0, T ]. Here d denote the number of spatial dimensions involved, while B f and B s represent the solid and the fluid domain respectively.
where ρ f is the density of the fluid at time t, σ f denotes the Cauchy stress tensor and b f represents external volumetric forces acting on B f . The variables a f and v f represent the acceleration and the velocity of the fluid at a given point x. In order to accommodate moving surfaces, typically at the fluid-solid interface, we express the conservation laws for the fluid media in an arbitrary Lagrangian-Eulerian (ALE) framework (for details, refer to Donea and Huerta [2003]). Within the ALE description of motion, the acceleration of a particle is given as where c := v f −v, withv being the velocity of the fluid mesh. The notation (·) χ denotes an observation made in the referential domain described by χ, relative to which both the material X and the mesh x can deform.

Conservation law for the solid
Let us express the solid body B s with a purely Lagrangian representation, which can also be inferred as a special case of ALE with coincident referential and material coordinates for all t.
The conservation of momentum for B s then follows, where ρ s is the solid density, σ s is the Cauchy stress tensor and b s is the cumulative external body force acting on B s . The acceleration a s , in terms of the displacement u s is given as It is often convenient to express Eq. (4) over a reference configuration, for which we will particularly consider the initial configuration (distinguished as B so ) at t = 0 such that where ρ so is the solid density at t = 0, F = ∂x/∂X is the deformation gradient and S s is the second Piola-Kirchhoff stress tensor.

Constitutive laws
For the fluid model, let us consider an incompressible Newtonian fluid, for which the shear stress is modeled as a linear function of the strain-rate tensor (symmetric gradient of the fluid velocity). Air and water are typical examples of Newtonian fluids. For such fluids, the relation between stress and velocity is expressed as, where p f is the pressure, µ f is the dynamic viscosity, v f is the velocity of the fluid and I is the identity tensor. For B s , we consider a rubber-like hyper-elastic solid capable of undergoing large deformations. In keeping consistency with the benchmark problems considered in this study, we employ the Saint Venant-Kirchhoff material model for the constitution of the solid, i.e., Here λ s and µ s represents the material parameters known as the Lamé constants, while E s is the Green-Lagrange strain tensor. Having formulated S s , the Cauchy stress tensor can be readily obtained by employing the push-forward operation with J = det(F ). The above discussed material models are chosen specifically to suit the benchmark configurations and represent only two of the many possible constitutive models for the continuum. It should be noted that the enrichment strategy discussed in this paper is a discretization strategy and can be easily used together with alternative material constitutions.

Initial, boundary and interface conditions
In order to close the continuum model, we need to specify a suitable combination of initial and boundary conditions for the continuum. For the fluid, we assume a velocity description v fo that essentially satisfies Eq.
(2) at t = 0. For the solid, we assume a stress-free initial state and refer to it as the reference configuration, i.e., x = X at t = 0. Hence the initial conditions are formulated as while the boundary conditions, comprising the Dirichlet and the Neumann boundaries, are Here n f and n s denote the outward normal vectors at the Neumann boundary surface for the fluid and solid media respectively, while the entities with an over-bar represent the imposed fields at the boundary. The interface boundary should ensure a continuous transfer of kinematic and kinetic field, which is enforced through the following conditions: The kinematic condition on the position of the grid points x demands a pure Lagrangian description at the interface. It is automatically satisfied across the entire B s , since conservative laws for the solid continuum are expressed in the Lagrangian frame of reference. For B f , this condition at the interface is accommodated through the ALE description. It should also be noted that the interfacial condition on the velocity field, together with Eq. (2), necessitates a divergence-free velocity of the solid interface. For a partitioned FSI solution procedure such a description is not always guaranteed and may lead to numerical artifacts (e.g. see Küttler et al. [2006]).

Weak form and the finite element model
The continuum model presented in the previous section is analyzed in the framework of the FEM. To obtain the discrete finite element equations, let us first express Eqs.
(1-2) and Eq. (6) in its variational, i.e., weak form. We define the following functional spaces for the unknown field variables (v f , p f , u s ) and their respective test functions (w, q, δu): where L 2 (·) represents a collection of functions that are square-integrable over the given domain, while H 1 (·) is a subclass of L 2 (·) comprising functions that additionally posses finite squareintegrable first-order derivatives. The weak form can now be obtained by multiplying Eqs.
(1-2) with the test function w and Eq. (6) with δu, and integrating over the entire domain. After certain mathematical manipulations, we arrive at the following weak form: Find v f ∈ S v , p f ∈ S p and u s ∈ S u , such that ∀ w ∈ V v , q ∈ V p and δu ∈ V u .

Spatially discrete stabilized formulation
Let us now discretize the spatial domains B f and B s with n f e and n se non-overlapping finite elements such that constitute a reasonable approximation of the respective continuous domains.
Here Ω e denotes a unique finite element identified with the index e, while the superscript h denotes a discretized entity. We further define a discrete subset of our trial and test function spaces, denoted as such that the chosen functions now only span in the vicinity of element e and vanish elsewhere. The discrete formulation, also known as the Bubonov-Galerkin finite element formulation, is readily obtained by formulating the weak form of Eqs. (22-24) over the discretized space.
The Bubonov-Galerkin formulations leads to a discretization scheme that is similar in character as the central-difference scheme. Central-difference schemes are known to be numerically unstable for problems where the advection operator dominates the diffusion term, such as convective fluid flows. These instabilities manifest themselves in the form of spurious node-to-node oscillations, primarily in the velocity field. Moreover for incompressible fluid flow, the peculiar form of Eq. (2), which is devoid of the pressure, leads to a saddle-point problem where the pressure merely acts as a Lagrange multiplier to satisfy the incompressibility constraint. Such problems are constrained by a solvability condition, known as the LBB (Ladyzhenskay-Babuška-Brezzi) condition, the violation of which leads to spurious oscillations in the pressure field.
Several solution strategies were developed to circumvent the numerical instabilities associated with the Bubonov-Galerkin form for problems involving incompressible fluid flow. For this study, we employ the consistently stabilized approach, the streamline upwind Petrov-Galerkin (SUPG) (Brooks and Hughes [1982]) to overcome the difficulties associated with the convective operator, while the pressure stabilizing Petrov-Galerkin (PSPG) (Hughes et al. [1986]) approach is used to circumvent the LBB constraint. Additionally, the least-square incompressibility stabilization (LSIC) (Hansbo and Szepessy [1990]) approach is utilized for ensuring mass conservation at the discrete level. The discrete stabilized formulation of the coupled problem is given as: while τ m and τ c are stabilization parameters which are also known as the intrinsic time scales. The optimal definition of the stabilization parameters has been an area of intense research and over the years several alternate definitions have been proposed in the context of stabilized finite elements for fluid dynamics (e.g. see Brooks and Hughes [1982], Hughes et al. [1986], Codina [2000], Franca and Valentin [2000], Tezduyar and Osawa [2000], ). Here, we employ the definition presented in  in the context of the variational multiscale approach. Thus, for τ m , we have with the tensor G defined by Here, ∂ξ/∂x can be construed as the inverse Jacobian matrix of the mapping between the parametric and physical element, ∆t is the discrete time step size and ν f = µ f /ρ f is the kinematic viscosity of the fluid. Additionally, τ m also depends on the constant m e , where m e = min 1 3 , C inv . 8 The value of C inv is chosen such that the stabilized weak form is coercive for a given discretization. Estimates for C inv for triangular and quadrilateral Lagrangian elements were formulated in Harari and Hughes [1992], which were later shown to be reasonable estimates for NURBSbased finite elements as well in Gamnitzer [2010]. For τ c , we have where Within a unique finite element e, the discrete value of the involved variables are obtained through a linear combination of basis functions chosen from the discrete functional space. We retain the isoparameteric concept and represent all variables with similar basis functions. Hence, the value of the test functions and the trial solutions within the element e are estimated as respectively, and a similar representation of the geometry is employed, i.e., The basis functions for the element e are collected in the vectors N andÑ, such that where n n denotes the total number of nodes and control points of element e. The definition of the individual basis functions N A is discussed in Section 3. The scalar and the vector variables are collected as and similar arrangements for δu e , v e f , p e f , u e s , x e and X e . The approximations of Eqs. (36-38) are substituted into Eqs. (26-28) and a finite element assembly operation is performed to obtain the discrete balance equations at the global level. Considering that the weak form is valid for all admissible choices of (w, q, δu), we obtain the nodal force balances at the non-Dirichlet nodes of the FSI model: The terms f con and f sint (and subsequently f SUPG and f PSPG ) are inherently nonlinear in nature due to the solid material model and the Eulerian description of the fluid, respectively. A Newton-Raphson method based predictor-multicorrector algorithm is used to linearize and obtain the solution of this system. A detailed derivation and the definition of the individual terms involved in Eqs. (43-45) is presented in Appendix A.

Temporal discretization and linearization
Eqs. (43-45) represent a set of first-and second-order ODEs in time, respectively. In order to integrate these equations in time, we employ the generalized-α method, which is an extension of the Newmark method. It was first presented in Chung and Hulbert [1993] in the context of second-order systems (solid mechanics) and was later extended to first-order Navier-Stokes in Jansen et al. [1999]. It has since been successively applied to many FSI problems (e.g., see Kuhl et al. [2003], Bazilevs et al. [2010], Eken and Mehmet [2016]).
Let us proceed by first expressing the nodal force balance of Eqs. (43-45) in a unified residual vector such that, where an over-dot denotes a time-derivative. The generalized-α method requires the evaluation of the residual at intermediate time-levels n + α m and n + α f to advance the system from time-level n to n + 1, i.e., where the intermediate time-levels are given as and consequently the intermediate description of the field variables is obtained aṡ Additionally, the discrete Newmark formulae are used to estimate the field description at time level n + 1, i.e., v n+1 Here, α m , α f , γ and β are real-valued parameters that determine the character of the timeintegration scheme. A subsequent Taylor expansion of a two-step time-integration scheme obtained through the generalized-α method reveals a second-order accurate method if and a subsequent evaluation of the eigenvalues associated with the amplification matrix of the scheme guarantees unconditional stability if 10 Using the spectral radius ρ ∞ of the amplification matrix of the two-step method, a family of methods were proposed in Chung and Hulbert [1993] to obtained second-order, unconditionally stable schemes for second-order system. It was suggested to have for second-order systems, while a similar analysis in Jansen et al. [1999] recommends for the first-order Navier-Stokes system. A value of ρ ∞ = 1 results in a non-dissipative scheme, where all frequencies are retained in the analysis, essentially requiring a very fine time discretization. While ρ ∞ = 0 corresponds to the asymptotic annihilation case, where high-frequency components are annihilated in a single solution step. It is evident from Eqs. (58-59), that in the context of FSI analysis, there is a mismatch between α m 2 and α m 1 . Such an arrangement has the tendency to lead to asynchronous evaluation of time-derivatives, which can subsequently lead to the loss of accuracy. We employ the definition given by Eq. (59) for both the fluid and the solid system, as suggested in . Such an arrangement guarantees consistent evaluation of all the time-derivatives and provides optimal numerical dissipation for the fluid problem, which is stiffer in nature. Moreover, with this choice it is ensured that the spectral radius for the second-order system remains well bounded and always stays below 1, hence a stable scheme. For all subsequent examples we consider ρ ∞ = 1/2.
As mentioned earlier, the semi-discrete system of Eqs. (43-45) is inherently nonlinear by nature. The generalized-α method is invoked in the framework of a Newton-Raphson method to linearize and subsequently advance the simulation in time. Using an iterative predictor-multicorrector approach, we expand the residual of Eq. (47) about a previously k th -iterate solution for the residual at (k + 1) th -iterate, which should ideally be zero, i.e., T is a vector composed of the nodal accelerations and fluid pressure dofs. In the predictor-multicorrector framework, the solution procedure is implemented as follows: Predictor: The counter k is set to zero and a reasonable prediction is made for the field variables. We employ a constant velocity predictor such that: Corrector: After the predictor step, the execution of the solver moves into an iterative corrector step for k = 0, 1, 2, . . . , k max or until an acceptable convergence is achieved. The multicorrector step conducts the following sequence recursively: 1: Collocation of field entities at the intermediate time-level: 2: Setup of the linear system of Eq. (60) using the intermediate field entities and the solution for the incremental improvements ∆ā for the k th -iteration: 3: Update of the field variables using the incremental vector ∆ā n+1,k : The convergence of the algorithm for the iteration k is tested by formulating a relevant norm which is a function of the incremental vector ∆ā n+1,k and the residual vector R n+1,k . If this value is less than an acceptable convergence criteria, the time step n+1 is considered converged, otherwise the flow of the program is transfered back to corrector-step 1 for the (k + 1) th corrector iteration.

The linearized monolithic coupled system
The construction and the solution of the linearized system of Eq. (73) constitute the most computational intensive part of the predictor-multicorrector algorithm. Once all the element level contributions have been calculated, the entries are moved to a global system through a finite element assembly operation. For a monolithic solution procedure, the global system prior to the enforcement of the coupling conditions is of the form: where K f and f f are the tangent and the residual force vector associated with the fluid media, while K s and f s are corresponding entities for the solid media. Detail derivations for these term are given in Appendix A. It should be noted that ∆ā f = ∆v T f , ∆p T f T , while ∆ā s = ∆ü s . In order to enforce the interfacial conditions at the fluid-solid interface in a finite element setting, let us isolate the interfacial dofs. Let I be the set of interfacial dofs for the respective media on which it operates. Consequently, we define A as a set that holds all non-interfacial dofs for the fluid, while the set B is comprised of all non-interfacial dofs for the solid media. The coupling conditions of Eq. (16) in a monolithic finite element setting are enforced as Eq. (81) constitute the fully-coupled linearized FSI system. Dirichlet boundary conditions can now be imposed on the coupled system and a solution to Eq. (81) yields the incremental improvements ∆ā for the iterate k.

Grid motion
For FSI simulations, the motion of the fluid-solid interface boundary needs to be accommodated smoothly within the fluid mesh B h f . At the discrete level, the ALE representation of Eq.
(1) needs to be augmented with a mesh update algorithm such that the motion of the interface boundary can be adjusted in B h f without excessive degradation and distortion of the computational mesh. Several mesh update algorithms exist in this regard. These include algebraic schemes (e.g. Yigit et al. [2008], Ahn et al. [2010]), Laplace smoothers (e.g. Löhner and Yang [1996]), spring/elastic medium analogy (e.g. Farhat et al. [1998], Johnson and Tezduyar [1994]) and mesh regeneration schemes (e.g. Johnson and Tezduyar [1999]), among others.  For this study, we employ a simplistic algebraic mesh motion algorithm, similar to the strategy described in Ahn et al. [2010]. With reference to Figure 1, let us consider a rigidly moving solid body B s . We identify three distinct regions in B f . The region Ω 1 contains all those points of B f whose radius from the center of B f doesnot exceed r 1 . Similarly, Ω 2 contains points of B f whose radius is smaller than r 2 . Note that Ω 1 and Ω 2 are overlapping. The remaining outer region of B f can then be identified as B f \Ω 2 . A purely Lagrangian description is imposed on the points included in Ω 1 . This is achieved by moving the points with the same velocity as the rigidly moving solid body B s . At the same time, points in the B f \Ω 2 are kept static. The region between these two domains (i.e., Ω 2 \Ω 1 ) then constitutes the buffer region, where points are moved so as to gradually absorb the largest displacements occuring at the periphery of Ω 1 to zero displacement of points in B f \Ω 2 . Several weighting functions ϑ can be formulated to interpolate the displacement of the points in the buffer region, the simplest of which is of the following linear form: Once the weighting associated with a particular point in the fluid domain has been determined, the new position and the velocity of the point are given as The scheme is easily extended to problems involving solid bodies undergoing deformations. This is achieved by coupling discrete points in B f to the closest point on the fluid-solid interface. The updated configuration of the candidate point can then be obtained as a weighted multiple of the coupled interface point, as done in Eqs. (83)(84). For the numerical examples considered in this study, we employ an explicit representation of the mesh update scheme, i.e., the mesh configuration of the domain B h f for the time step (n + 1) is determined from the interface description of the state at the time step n. Such a treatment allows mesh updates without the need of solving additional equations and is observed to be reliable for the time step sizes considered in this study.

Isogeometric enrichment of the discretized domain
Isogeometric analysis has the potential to increase the accuracy of finite element simulations, but it simultaneously increases the required computational effort (see Section 1 and Section 4). Moreover, pure isogeometric discretization of volumetric domains is a challenging task that is still not solved for arbitrary geometries. Isogeometric enrichment mixes isogeometric elements (spline interpolation) with standard finite elements (Lagrange interpolation) in order to overcome these deficiencies. Within this framework, two general isogeometric enrichment strategies are described in Section 3.1, while the construction of the enriched elements for interfacing standard elements with isogeometric regions is presented in Section 3.2.

Enrichment strategies
The discretization of an entire domain with standard finite elements (see Figure 2a)) leads to an approximative representation of the geometry and to low-order continuity throughout the whole domain. Such drawbacks are addressed by isogeometric enrichment. The key idea is to identify and discretize regions that are most crucial for the accuracy of the physical problem with accurate and smooth isogeometric elements. These crucial regions of interest can be identified from a-priori knowledge of the physical problem being modeled. The remaining domain is discretized with efficient and low-order standard elements. The following two cases can be distinguished: 1 Figure 2: Discretization of an exemplary two-dimensional domain with and without isogeometric enrichment. The enrichment strategies can be applied analogously in three-dimensions.

Surface enrichment
Isogeometric surface enrichment of the standard FEM (see Figure 2b) results in the exact representation of the underlying surface geometry. Moreover, such an enrichment provides a representation, where the field and the solution variables at the surface are continuously differentiable. In order to obtain a surface enriched discretization, a layer of isogemeotrically enriched elements (see Section 3.2) is created at the surface. These special elements interface smooth spline representation at the surface with standard finite elements in the bulk volume. Surface enrichment yields significantly more accurate results and is particularly beneficial for surface-dominated engineering problems such as contact, friction and adhesion (see Sauer [2011Sauer [ , 2013 for two-dimensional and Sauer [2014, 2015] for three-dimensional analysis). At the same time, the number of dofs is only slightly increased compared to pure standard finite element discretizations. For coupled problems, isogeometric surface enrichment additionally offers the possibility of employing formulations that demand continuous surface representations without performing IGA for the entire coupled domain. An example is a wind turbine discretized with rotation-free shell elements (e.g. formulated in Kiendl et al. [2009] and Duong et al. [2017]) in a discretized fluid domain that is composed largely of standard finite elements. However, for problems where solution gradients are dominant within the volume, surface enrichment offers only marginal accuracy gains as demonstrated in Rasool et al. [2016] and Harmel et al. [2017].

Zone enrichment
Zone enrichment (see Figure 2c) denotes isogeometric discretization of certain volumetric regions of the domain, while the remaining domain is discretized with standard elements. A layer of isogeometrically enriched elements (see Section 3.2) interfaces these regions. In physical problems, large gradients of the solution variables typically occur only in certain regions of the domain, such as boundary layers in flow problems or stress-concentrations in structural problems. Accurate representation of these potentially nonlinear gradients with smooth spline basis functions increases the accuracy of the numerical solution significantly compared to standard finite element discretizations, while the computational effort is substantially less than for pure isogeometric discretizations. Such an enrichment promises significant potential for many engineering applications such as fluid flow (see Rasool et al. [2016]) and heat transfer (see Harmel et al. [2017]).

Isogeometrically enriched elements
The isogeometrically enriched finite element is crucial to constructing surface enriched and zone enriched discretizations. It has a spline-based representation at one of its surfaces where it can interface with a pure IGA representation, while the opposite surface has a Lagrange-based representation where it connects with standard finite elements. For simplicity, the construction of two-dimensional elements are presented first, while the extension of the concept to the threedimensional setting is discussed afterwards. An isogeometrically enriched quadrilateral element can be obtained by replacing an edge of a standard quadrilateral element (i.e., Lagrange basis functions of order p) with a spline curve of order q. Without loss of generality, let us locate the spline curve at η = −1 (see Figure 3a and Figure 3b). We denote the respective element type as LpSq, where q denotes the order of the spline curve (B-Spline, NURBS, T-Splines etc.) used to enrich the element at η = −1. An enriched element with linear Lagrange basis functions (i.e., L1Sq) will possess two standard finite element nodes and q + 1 control points of the spline curve. The basis functions of such an element are Here L p i (·) denotes the i th Lagrangian basis function of order p, while S q j (·) represents the j th spline basis function of order q. For the elements in Figure 3a and Figure 3b, i = [1, 2] and j = [1, 2, . . . , q + 1].
A three-dimensional enriched element can be obtained by replacing one face of a standard hexahedral element (e.g. the surface at ζ = −1 surface) with a spline surface (see Figure 3c and Figure 3d). The shape functions corresponding to the control points of the enriched surface will be a tensor product of the spline basis (in ξ-and η-direction) and Lagrangian basis (ζ-direction). Extending the earlier defined nomenclature, three-dimensional isogeometrically enriched elements are denoted as LpSqSr where q and r are the order of the polynomial defining the enriching spline surface. A three-dimensional element with linear Lagrangian basis functions (i.e., L1SqSr) will have (q + 1)(r + 1) control points and four interpolatory nodes.
Isogeometric enrichment of standard elements leads to continuous representation of the field and the solution variables at the enrichment interface. The enriched basis functions retain the partition of unity over the entire element. Moreover, at the enriched surface, the enriched element type fulfills the variation diminishing property and its basis functions are always pointwise non-negative, similar to spline basis functions.

Numerical examples
The monolithic FSI model and the proposed IGA enrichment strategy, discussed in the previous sections, is validated and assessed in this section for different FSI benchmark problems. These include two-dimensional fluid-flow in a box cavity with a flexible base, the propagation of a pressure pulse in a deformable pipe and channel flow past a circular cylinder with a flexible tail. The consistency of the implemented model and the influence of the zone enrichment strategy is assessed through convergence analysis of relevant error norms for successive spatial refinements. Due to the non-availability of an exact solution and dearth of robust mesh convergence reference studies, relative error measures are formulated in relation to the obtained numerical solution of the finest discretization. Only for the third example, the obtained numerical solution is benchmarked with an available reference solution of very fine discretization. The response from the finest discretizations are also provided for future benchmark studies.

Lid driven cavity with a deformable base
Internal fluid flow within a square cavity, due to the motion of the top surface, is a popular benchmark problem for fluid flow solvers. A modification to the pure fluid flow setup was proposed in Wall [1999] to account for a base that could deform under the influence of fluid forces, essentially setting up an FSI example. Since then the problem has been used as a demonstration example for many numerical FSI studies (Mok et al. [2001] Table 1. To avoid numerical singularities in the pressure field associated with typical cavity flow benchmark studies (e.g. see Donea and Huerta [2003]), small segments close to the lid are modeled as traction free boundaries on both lateral faces, consequently allowing fluid inflow and outflow at these segments. Figure 4 illustrates this setup in its entirety.
For the discretization, four enrichment cases are considered. A discretization composed entirely of second-order isogeometric spline elements S2S2 (i.e., with p = q = 2) for both  is continually discretized with S2S2 IGA elements. A further special case of the zone enriched case is the surface enriched discretization, where only the interface is provided with NURBS representation using L1S2 elements, while the remaining bulk is discretized with L1L1 elements. A discretization with no IGA enrichment is termed as the pure L1 case and is composed entirely of L1L1 Lagrange elements. Figure 5 shows the different enriched discretization at the coarsest level, while the respective element and dofs population data for each considered mesh are tabulated in Table 2.
Numerical simulations are performed over the interval t ∈ [0, 50 s] with a time step size of 0.1 s. Figure 6 depicts the fluid pressure field at different stages of the simulation. Additionally, the vertical displacement u 2 of three reference points on the fluid-solid interface, points A, B and C as shown in Figure 4, are recorded over the course of the entire simulation. The obtained response for the evolution of the vertical displacement of these reference points over the entire time interval is shown in Figure 7. This response is similar in character to those reported in Mok et al. [2001] and Gerbeau and Vidrascu [2003]. For subsequent convergence analysis, the numerical solution corresponding to the finest pure IGA discretization (i.e., m = 5) is    considered as the reference solution.
In order to assess the convergence behavior of the obtained dynamic numerical response, we formulate a Fourier approximation of the resulting discrete data and compare the relative differences in the approximated Fourier coefficients. Such a treatment allows us to encapsulate relative differences in the amplitude, the frequency and the phase-shift between the reference and the obtained solution in a single measure. The convergence of this measure is then studied for successive mesh refinements. A Fourier representation of an evolutionary function is given as  where with 2L being the period of the function f (t). For a finite summation of Fourier coefficients a k and b k (i.e., k = 1, . . . , n), only an approximation of f (t) is obtained. The value of n is reflective of the number of distinct frequencies inherent in f (t). Using the discrete data and a predetermined value of n, a Fourier approximation comprising the Fourier coefficients can be obtained using a curve fitting algorithm (e.g., see Bates and Watts [1988]). For the analysis performed in this study, we employ the nonlinear least square method to fit the obtained numerical data to a Fourier approximation. Once the Fourier coefficients have been approximated, the relative difference between the reference and the obtained solution is assessed using the following error measure: where the superscript "ref" denotes Fourier approximation using the reference solution, while the superscript h denotes Fourier approximation using the obtained solution.   Figure 7 were used to construct a Fourier approximation with n = 1. A subsequent increase in the value of n essentially results in plots identical to Figure 8. It is evident from Figure 8 that the relative error decreases consistently for all enrichment cases, reflecting a consistent behavior. The use of S2S2 elements in the enriched zones results in a significantly improved solution for the zone enriched case, as compared to the pure L1 and the surface enriched cases. The region near the fluid-solid interface and the fluid domain near the lid of the cavity experiences rapid movements of the boundary surface and therefore inherently retains strong velocity gradients. An IGA zone enrichment only in these areas improves the analysis considerably as Figure 8 shows.

Pressure pulse in a deformable tube
For the second benchmark example, we consider the case of incompressible fluid flow in a circular tube with a traveling pressure wave generated at the inlet. Although initially proposed as a simplified model for a human cardiovascular system in Formaggia et al. [2001], the problem setup has since been widely employed with different geometrical and boundary condition implementations in many FSI numerical studies (e.g., see Gerbeau and Vidrascu [2003], Küttler and Wall [2008], Gee et al. [2011], Malan and Oxtoby [2013], Mayr et al. [2015], Eken and Mehmet [2016], Eken and Sahin [2017] Table 3. The tube wall at the inlet and the outlet are held clamped over the entire simulation. The inlet and the outlet for the fluid domain are modeled as traction free boundaries, where σ f · n f = 0. For a small time interval, i.e. for t < 0.003 s, a fluid pressure (p pulse ) of 1300 N/m 2 is applied at the inlet. This impulse generates a pressure pulse that continues to travel along the length of the tube, while interacting with the deformable tube. When the pulse reaches the outlet boundary, part of it gets reflected back towards the inlet due to the clamped nature of the outlet boundary.  Due to the geometry and the applied boundary conditions, the problem is expected to yield a symmetric solution within the laminar flow regime. Therefore only a single quadrant of the problem domain is simulated to save computational cost, while retaining a three-dimensional analysis configuration at the same time. Four IGA enrichment cases, as shown in Figure 9, are considered for this benchmark example. The surface enriched case comprises a discretization, where the tube's outer surface and the fluid-solid interface are discretized with L1S2S2 elements to have a NURBS representation of the curved surfaces, while the remaining bulk is composed of linear Lagrange elements (L1L1L1). The zone1 case is an enrichment of the surface enriched case where the solid domain B h s is discretized entirely with second-order IGA hexahedral elements (S2S2S2). For the zone2 enriched case, only a thin layer of IGA hexahedral elements are added in the boundary layer region of the fluid domain, while the major fluid bulk is discretized with linear Lagrange elements. IGA hexahedral elements constitutes the discretization of B h s for zone2 as well. Figure 9 shows the four enrichment cases at the coarsest level, while mesh statistics are tabulated in Table 4. Figure 10: Pressure pulse in a deformable tube: Contours of fluid velocity magnitude ||v f || = v 2 1 + v 2 2 + v 2 3 at different stages for the zone2 enriched discretization (m = 2). The displacement of the solid is increased by a factor of 10 for illustration purposes.
Numerical simulations are performed over the interval t ∈ [0, 0.02 s] with 85 time steps. During Figure 11: Pressure pulse in a deformable tube: Contours of pressure at different stages for the zone2 enriched discretization (m = 2). The displacement of the solid is increased by a factor of 10 for illustration purposes. Figure 12: Pressure pulse in a deformable tube: Evolution of the radial u r and the axial u a displacement of a point on the tube's inner surface at halfway length of the pipe (pure L1 case). this interval, a complete reflection of the pulse from the outlet to the inlet is observed. Figure 10 illustrates the fluid velocity field in the tube's interior together with the deformations of the solid domain, as the pulse travels from the inlet to the outlet. Pressure contours are shown in Figure 11. In order to assess the convergence behavior of the enrichment cases, we measure the radial and the axial displacement of a point (u r and u a , respectively) located on the tube's inner surface at halfway length. The evolution of these entities are shown in Figure 12 for the pure L1 discretizations, in comparison to the findings from Eken and Mehmet [2016]. It is evident from the figure that the numerical solution will converge to a unique finite element solution as the mesh is refined. The difference between the obtained results and the findings from Eken and Mehmet [2016] can be attributed to the different boundary conditions at the outlet. Zero fluid pressure was imposed at the outlet in Eken and Mehmet [2016], whereas we impose a traction free outlet surface. It is also observed that all the considered enrichment cases tend to converge to identical numerical solutions with spatial refinements. This is evident from the response 24 history of u r for different enrichment cases at the finest level, as shown in Figure 13.  To perform a comparative analysis for the different IGA enrichment strategies considered for this example, we use the error measure of Eq. (90) to judge the relative difference between the obtained numerical solutions. For this purpose, we consider the response history of u r and u a obtained from the finest discretization (i.e., m = 3) for the zone enriched case as the reference response. The numerical solutions are fitted to a Fourier series with n = 7, where the first acceptable approximation of the reference solution is obtained. Mesh convergence analysis for the parameters E ur and E ua are shown in Figure 14. It is evident, that the solution obtained from the zone1 enriched discretization, with S2S2S2 IGA elements only in B h s , is much better than the ones obtained from the surface enriched and the pure L1 discretizations. However, the most important observation is regarding the performance of zone2 enriched discretization, where the obtained results are significantly more accurate than all other enrichment cases. This clearly advocates the benefit of incorporating a thin layer of S2S2S2 IGA elements in the boundary layer region of B h f .

Flow past a circular cylinder with flexible tail
Laminar flow past a circular cylinder with a thin flexible tail is a popular two-dimensional numerical benchmark for FSI solvers. The benchmark configuration was first proposed in Turek and Hron [2006], and has since been adopted by many researchers to validate different FSI solvers (e.g., Heil et al. [2008], Kollmannsberger et al. [2009], Bhardwaj and Mittal [2012], Lee and You [2013] among others). Three different test cases -based on the fluid velocity at the inlet and the material parameters of the deformable tail -were proposed in the original paper of Turek and Hron [2006]. Here, we investigate the more challenging unsteady test cases FSI-2 and FSI-3 of the original proposal.
while at the fluid-solid interface The nature of the fluid flow within the domain is essentially controlled through the fluid velocity at the inlet i.e., with v in , while the behavior of the solid is regulated through the material parameters ρ so and µ s . For this example we consider four IGA enrichment levels, which are similar in construction as those of the deformable tube example of Section 4.2. These include the pure L1 case (with only L1L1 elements), the surface enriched case (with L1S2 elements only at B f ∩ B s ), the zone1 enriched case (with S2S2 elements only for B s ) and the zone2 enriched case (with zones of S2S2 elements for B f and B s ). Spatial meshes corresponding to the coarsest refinement level for different enriched discretization are shown in Figure 16, while element and dofs population data is tabulated in Table 6.   Figure 15) is recorded over the entire simulation. Figure 17 shows contours of the velocity magnitude of the fluid at different deflection stages of the tail for the two investigated cases. Contours for fluid pressure at identical time-instances are shown in Figure 18. Although the FSI-3 test case represents a higher Reynolds number flow (i.e., Re = 200), the FSI-2 test case (with Re = 100) is the more challenging one as the solid structure is softer and thus undergoes much larger deformations. The difference in the fluidinduced deformations in the flexible tail for the two test cases is also evident from Figure 19, where the evolution of the displacement of point A is plotted in terms of components u 1 and u 2 .    To gauge the effectiveness of different enrichment strategies used within this benchmark problem, we compare the error in the Fourier norm (see Eq. 90) of the obtained numerical solutions with a reference solution. For this purpose, we use the benchmarking data available at fea [2016] as our reference solution, against which all the obtained numerical solutions are compared with. The response history of the horizontal (u 1 ) and the vertical (u 2 ) displacement of the reference point A corresponding to the last two seconds of the simulation is used for this comparison. Within this interval, the displacement of the point A becomes truly periodic. A Fourier series with n = 2 offers a sound approximation, as any further increase in the number of Fourier coefficients n yields identical approximation of the reference solution, and hence identical error plots. The results for this comparison are shown in Figure 20 for the FSI-2 test case and Figure 21 for the FSI-3 test case.
Figure 20: Flow past a circular cylinder with flexible tail: Mesh convergence behavior of the error norm in the Fourier coefficients, i.e., E u 1 and E u 2 , for the FSI-2 benchmark case. Figure 21: Flow past a circular cylinder with flexible tail: Mesh convergence behavior of the error norm in the Fourier coefficients, i.e., E u 1 and E u 2 , for the FSI-3 benchmark case.
It is evident from the figures that the implemented scheme results in a consistent method. The consistent reduction in the Fourier error norms (E u 1 and E u 2 ) reflects that with each subsequent increment of the refinement level of the spatial mesh, the numerical solutions get closer to the reference solution reported in fea [2016]. This is particularly true for the FSI-2 case, as seen in Figure 20. However, for the FSI-3 test case, the behavior of the error norms for the finer meshes (e.g., at m = 2, 3) reflect a minor difference between the reference solution and the converged solution from the obtained results. For both test cases, the zone2 enriched discretization yields 30 results that are remarkably improved than all other enrichment cases. While all other enriched discretizations results in similar solutions (including the zone1 enriched with S2S2 elements for only B s ), the zone2 enriched results are significantly more accurate and demonstrate higher order of convergence than others, particularly for the FSI-2 case. Hence for some applications, only analyzing the structure domain with IGA elements may not be enough and the inclusion of IGA zones for crucial fluid regions has the potential to offer significant accuracy gains.

Conclusion and outlook
In this paper the efficacy of enriched finite element discretizations -where isogeometric analysis is blended with classical Lagrangian finite elements -is investigated in the context of fluid-structure interaction problems. IGA enriched finite elements offer a mechanism to enrich the classical finite element method with localized IGA. Such a treatment is extremely beneficial for FSI problems, where often large computational domains are employed to mitigate farfield boundary effects for the fluid domain. In this regard, careful IGA enrichment of only the crucial regions will result in improved numerical solutions with minimal increase in computational cost. Moreover, the strategy also provides a mechanism to analyze different interacting physical domains with different discretizations approaches in a conforming mesh setting (e.g., IGA for the solid and classical Lagrangian finite elements for the fluid).
where B v in Eq. (106) is given as, while the finite element divergence operator D is formulated as, with such that The discrete residuals are evaluated through Eqs. (29-30). The divergence of σ h f can be obtained as with C = µ f C 1 , C 2 , . . . , C nn , such that (113) We can now collect the individual terms representing Eq. (97) such that The necessary ingredients for computing the element level representation of Eq. (115) have been defined. The element level contributions are then relocated to a global system using a finite element assembly operation, which ultimately results in where and w contains the values of the test functions w for all the nodes in B h f . Noting that the weak form holds for all admissible choice of test functions, we arrive at Eq. (43), i.e., Adopting a similar procedure for Eq. (27), the discrete finite element equation is of the form where The individual terms can be collected and then subsequently moved to a global system using a finite element assembly operator. Further noting that the global discrete system holds true for all admissible q, we arrive at Eq. (44), i.e., Applying a finite element discretization to Eq. (28), we obtain Here, Ω e o represents a unique finite element obtained from the discretization of the reference configuration. The l th surface element on the boundary ∂ t B h so is represented as Γ l o , while n sl represents the total number of surface elements on ∂ t B h so . The individual terms of Eq. (123) can be further expressed as where B L is the strain-displacement matrix, which in R d is taken as such that 3 F 23 N A,3 F 33 N A,3 F 12 N A,3 + F 13 N A,2 F 22 N A,3 + F 23 N A,2 F 32 N A,3 + F 33 N A,2 F 11 N A,3 + F 13 N A,1 F 21 N A,3 + F 23 N A,1 F 31 N A,3 + F 33 N A,1 F 11 N A,2 + F 12 N A,1 F 21 N A,2 + F 22 N A,1 F 31 N A,2 + F 32 N A,1 while the value of the second Piola-Kirchhoff stress tensor S h s can be obtained using Eq. (8) with E s = (F T F − I)/2. Subsequently performing the finite element assembly procedure, we obtain the global nodal force balance of the solid system, represented by Eq. (45), i.e.,

A.2 The tangent matrices
The nodal force balance system of Eqs. (43-45) consists of a system of nonlinear equations. Using a predictor-multicorrector algorithm, a Newton-Raphson approach is employed to linearize the nonlinear system, and subsequently advanced in time using the generalized-α-method. A major component of this solution approach is the setup and solution of the linearized system of Eq. (73).
In-order to derive the terms associated with this system, let us expand Eq. (73) in terms of its components, i.e., Adopting the above given definition, we express the tangents associated with K f as where C = 2 µ f I, while the matrices For the tangent matrices associated with the solid model, we note that using the chain-rule and the definitions from the generalized α-method, we have Hence the element-level tangent matrix associated with K s is expressed as   (143) is the geometrical stiffness contribution (S IJ is a component of S s ), while is the material stiffness with D being the fourth-order constitutive tensor, such that S s = D : E s . For the considered Saint-Venant material model, the constitutive tensor is D = λ s I ⊗ I + 2 µ s I .