An isogeometric b-rep mortar-based mapping method for non-matching grids in fluid-structure interaction

In this study the isogeometric B-Rep mortar-based mapping method for geometry models stemming directly from Computer-Aided Design (CAD) is systematically augmented and applied to partitioned Fluid-Structure Interaction (FSI) simulations. Thus, the newly proposed methodology is applied to geometries described by their Boundary Representation (B-Rep) in terms of trimmed multipatch Non-Uniform Rational B-Spline (NURBS) discretizations as standard in modern CAD. The proposed isogeometric B-Rep mortar-based mapping method is herein extended for the transformation of fields between a B-Rep model and a low order discrete surface representation of the geometry which typically results when the Finite Volume Method (FVM) or the Finite Element Method (FEM) are employed. This enables the transformation of such fields as tractions and displacements along the FSI interface when Isogeometric B-Rep Analysis (IBRA) is used for the structural discretization and the FVM is used for the fluid discretization. The latter allows for diverse discretization schemes between the structural and the fluid Boundary Value Problem (BVP), taking into consideration the special properties of each BVP separately while the constraints along the FSI interface are satisfied in an iterative manner within partitioned FSI. The proposed methodology can be exploited in FSI problems with an IBRA structural discretization or to FSI problems with a standard FEM structural discretization in the frame of the Exact Coupling Layer (ECL) where the interface fields are smoothed using the underlying B-Rep parametrization, thus taking advantage of the smoothness that the NURBS basis functions offer. All new developments are systematically investigated and demonstrated by FSI problems with lightweight structures whereby the underlying geometric parametrizations are directly taken from real-world CAD models, thus extending IBRA into coupled problems of the FSI type.


Introduction
Computer-based simulations are playing an ever increasing role in the engineering design and production process as they offer reliable predictions based on computational mod-els. These computational models have been traditionally obtained using the standard Finite Element Method (FEM) [1,2] which is typically applied to Computational Structural Dynamics (CSD) or the Finite Volume Method (FVM) [3,4] which is typically used in Computational Fluid Dynamics (CFD). Regardless of whether the aforementioned discretization methods are applied on structured or unstructured meshes, there is a considerable effort involved in the mesh generation which also leads to a discrete representation of the model's geometry. However, the accuracy of the geometric model may play a decisive role in various engineering applications such as, form-finding of membranes [5] and buckling analysis of shells [6].
Isogeometric Analysis (IGA) is a modern numerical method which uses the geometric basis of a parametrized model in Computer-Aided Design (CAD) for the approximation of the unknown fields [7]. Typically, Non-Uniform Rational B-Spline (NURBS)-based IGA is used which means that the underlying CAD model is parametrized using the NURBS basis functions [8]. The use of the NURBS basis functions for analysis typically results in smooth and high order field approximations. Moreover, NURBS-based models are standard in CAD and are therefore favoured for bridging design and analysis [9,10]. It is worth noting that other splines such as the T-Splines [11] can be also used in the context of IGA, see in [12] for more information. NURBS-based IGA was initially applied on single patch untrimmed geometries for the Kirchhoff-Love shell problem in [13] and for the Reissner-Mindlin (RM) shell problem in [14]. The application of IGA on conforming multipatch geometries for the Kirchhoff-Love shell problem is detailed in [15]. It was then extended to non-conforming untrimmed multipatch geometries in [16][17][18][19]. Subsequently, Isogeometric B-Rep Analysis (IBRA) was introduced in [20] for nonlinear shell structures on arbitrary CAD geometries involving non-watertight trimmed multipatches which are standard in real-world CAD models. Subsequently, the application of IBRA in membrane analysis was shown in [21] and generally to lightweight structures in [22]. In this study, the application of IBRA to linear static and modal analysis of the National Renewable Energy Laboratory (NREL) phase VI wind turbine with flexible blades [23,24] is demonstrated to highlight its application range.
Fluid-Structure Interaction (FSI) [25] forms a prominent category of surface coupled problems which govern the mutual interaction between a fluid flow and a flexible structure. A common interface between the fluid and the structural domains has to be identified in order to appropriately define the interface conditions [26,27]. Confining ourselves to partitioned FSI, diverse numerical methods are used for solving each of the CFD and CSD problems whereas the interface constraints are fulfilled in an iterative manner in case of strong coupling which is herein employed. Moreover, the partitioned Gauss-Seidel (GS) approach is herein used, whereby a matching time discretization is assumed and the underlying solvers exchange information within each time step until convergence is achieved, see in [24,28]. This allows for efficient methods targeted for each of the CFD and CSD problems to be employed [27]. As a consequence, the interface discretizations between the fluid and the structural subdomains are typically not matching and therefore methods for transferring fields between these non-matching interfaces have to be developed. Such methods comprise Nearest Neighbor, Nearest Element, Barycentric Interpolation and Mortar-based formulations, see in [29,30]. Amongst these, the mortar-based mapping method is found to be the most accurate and robust method especially for highly diverse interface discretizations which are common in FSI where the fluid interface discretization is typically much finer than the corresponding structural interface discretization due to the resolution of the fluid boundary layer [31]. The mortar method was initially applied to elliptic boundary value problems on multiple subdomains with non-matching discretizations in [32]. It was then employed in the context of FSI problems in [33] with standard finite element discretizations and in [34] with isogeometric structural discretizations.
In this contribution, a mortar-based mapping method on trimmed multipatch NURBS geometries stemming directly from real-world CAD geometries is developed and systematically evaluated which is herein targeted to partitioned FSI, see also in [34]. This method was firstly introduced in the dissertation [35] and this contribution serves as a complement with more detailed elaboration, presentation and evaluation of the method along with additional findings and conclusions. The treatment of the continuity conditions across the NURBS-patch interfaces is done using a Penalty method in a similar fashion as in [20] since multiple trimmed NURBS patches are considered. The proposed B-Rep mortar-based mapping method can be used to perform FSI between a low order discretized CFD problem (FVM) and an IBRA discretized CSD problem. Alternatively, the trimmed multipatch NURBS geometry can be used as an Exact Coupling Layer (ECL) for FSI simulations between a low order discretized CFD problem (FVM) and a low order discretized CSD problem (FEM) taking advantage of the smooth and high order NURBS basis functions for smoothing the fields (displacements and tractions) which are transformed from one low order discretization to the other.
The study is complemented with examples ordered in a sequence of increasing complexity. Firstly, the cavity FSI benchmark case is used since it is standard in literature for the verification of FSI methods. Due to the simplicity of its geometry and given that relatively large deformations are enabled, it allows for a detailed evaluation and verification of the proposed isogeometric B-Rep mortar-based mapping method. Then, the example of an inflatable hangar in numerical wind tunnel is employed in the frame of membrane structural analysis in FSI. Lastly, the NREL phase VI wind turbine with flexible blades in numerical tunnel is demonstrated, thus extending IBRA for Kirchhoff-Love shells in FSI.
The study is organized as follows: "Computer-aided design using non-uniform rational b-splines" section provides a compact but comprehensive introduction to trimmed NURBS multipatch geometries and how these are treated within CAD. "Isogeometric lightweight structural analysis on trimmed multipatches" section serves as a comprehensive review of the latest developments in IBRA where the most important methods concerning the coupling between the multiple trimmed patches are revised. This chapter is complemented with a new large scale numerical example within IBRA to support the theory. "Isogeometric B-Rep mortar-based mapping method on trimmed multipatches" section introduces the proposed isogeometric b-rep mortar-based mapping method which is the core of this study. The method is elaborated in very detail where both its variational statement and its discrete equation system are presented. Moveover, a detailed guide on how the interlying integrals are evaluated is given whereby part of the NREL phase VI wind turbine blade's geometry is used for the demonstration of the algorithm. "Fluid-structure interaction" section serves as a very compact introduction to partitioned FSI with strong coupling regarding the satisfaction of the FSI interface constraints. Lastly, "Numerical results" section presents a series of numerical examples ranging from benchmark to realworld engineering applications demonstrating the applicability of the proposed framework of CAD-integrated analysis for multiphysics problems.
whereq = 0, . . . ,p. Thus, the B-Spline basis functions Np ,i are obtained by a recursive constructionq = 0, . . . ,p where N 0,i = 1 for θ ∈ [θ i ,θ i+1 [ while N 0,i (θ ) = 0 identically elsewhere concerning the constant basis functions. Moreover, the definition 0/0 = 0 is assumed in Eq. (2). The B-Spline and consequently the NURBS basis functions attain C ∞ -continuity within each knot span ]θ i ,θ i+1 [⊂γ and Cp −k i -continuity across knotŝ θ i wherek i stands for the multiplicity ofθ i in . In this way, given a set of pointsX i , i = 1, . . . , n ∈ N in R α Euclidean space, known as the Control Polygon, the corresponding NURBS curve C :γ → R α is parametrically constructed as, at each parametric location θ ∈γ . Within this study, open knot vectors are considered, namely, the first and the last knots,θ 1 ,θ m ∈ respectively, havep + 1-multiplicity so that the curve is interpolated by the control polygon at the beginning and at the end. In the sequel, saying that a NURBS basis Rp ,i , i = 1, . . . , n is of polynomial orderp implies that its underlying B-Spline basis Np ,i attains polynomial orderp.

Non-uniform rational b-spline surfaces
The two-dimensional NURBS basis functions Rp 1 ,p 2 ,i,j , with i = 1, . . . , m 1 ∈ N and j = 1, . . . , m 2 ∈ N, are constructed using the two-dimensional B-Spline basis functions Np 1 ,p 2 ,i,j in a similar fashion as the NURBS curves in Eq. (2), namely, Rp 1 ,p 2 ,i,j (θ 1 , θ 2 ) =ŵ i,j Np 1 ,p 2 ,i,j (θ 1 , θ 2 ) n 1 k=1 n 2 l=1ŵ k,l Np 1 ,p 2 ,k,l (θ 1 , θ 2 ) where m α ,p α and n α stand for the number of knots at each knot vector α , the polynomial orders and the number of one-dimensional B-Spline basis functions in θ α -parametric direction, respectively. Surface's parametric imageˆ is then defined by the square domain spanned in θ α -direction by knot vector α . Additionally,ŵ i,j stands for the weight associated with NURBS basis function Rp 1 ,p 2 ,i,j . The two-dimensional B-Spline basis functions Np 1 ,p 2 ,i,j in turn, are constructed as a tensor product of the underlying one-dimensional B-Spline basis functions, that is, In this way, given a net of pointsX i,j in R 3 known as the Control Point Net, the corresponding NURBS surface S :ˆ → ⊂ R 3 is constructed by, Since open knot vectors α are considered herein, see "Non-uniform rational b-spline curves" section, surface interpolates the four corners of the control point net. The latter property along with the affine covariance of the NURBS basis functions, allows for the application of strong Dirichlet boundary conditions at boundaries of untrimmed patches within NURBS-based IGA. In the sequel, assumed is an one-to-one a map (i, j) → k such that Rp 1 ,p 2 ,i,j → Rp 1 ,p 2 ,k with k = 1, . . . , n 1 n 2 for a sequential ordering of the NURBS (or B-Spline) basis functions which is necessary for the construction of the discrete equation systems.

Trimmed multipatch NURBS surfaces
In CAD, the trimming of a NURBS patch (i) is performed using a set of NURBS trimming curvesγ (i) j defined in its parametric spaceˆ (i) , see in [20,37] for more information. In this way, trimming loops may be defined, each of which consisting of a sequence of unioriented trimming curves. Accordingly, parts lying outside these loops are non-visible and hence a lot of flexibility is added in describing arbitrary shapes in Euclidean space originating from generic models. However, in real world engineering practice, multiple trimmed surfaces are considered to accurately describe large scale models, such as cars, ships, airplanes, etc. Let (i) , with i = 1, . . . , n s ∈ N, be a non-overlapping domain- Fig. 1 Computer-aided description of surfaces: Trimmed NURBS multipatch surfaces [35] decomposition of , meaning that, I being the set of all pairs (i, j) where i, j = 1, . . . , n s with i = j and let n i ∈ N be the number of non-empty sets γ (i,j) i . Each of the surface patch subdomains (i) has a NURBS parametrization S (i) :ˆ (i) → (i) as per Eq. (6) and a set of trimming curvesγ t being the number of curves trimming patch (i) . An example of a trimmed multipatch geometry is depicted in Fig. 1 along with the distinct patch parametric spaces and the corresponding trimming curves. In terms of CAD, each interface γ (i,j) i has a unique representation from each of the neighboring patches, namely,γ (i) j andγ (j) i in the parametric spaces of patches (i) and (j) , respectively. For real world engineering applications, this distinct representation of the patch interfaces is in general not identical in the physical space, that is, γ i . The accuracy of the interface parametrization from each neighboring patch is then controlled by a tolerance which in most CAD software is user-defined. However, when applying IGA on trimmed NURBS multipatches, it is important that the trimming curves γ (i) j and γ (j) i representing the interface γ (i,j) i are identified so that the evaluation of the interface integrals, accounting for the continuity enforcement across the multipatches, can be performed. Most CAD software provide this identification of the patch interfaces through a sew option where the topological information of the geometry is generated.

Isogeometric lightweight structural analysis on trimmed multipatches
Isogeometric B-Rep analysis (IBRA) firstly introduced in [20] allows for performing Isogeometric Analysis (IGA) on real-world CAD models which involve trimmed multipatches. In this section the isogeometric analysis of lightweight structures on trimmed NURBS multipatches is briefly presented. Accordingly, membranes and thin shells of Kirchhoff-Love shell type are used, see in [38] for more information.

Differential geometry of surfaces
Herein, a brief introduction to the differential geometry of surfaces is provided and the underlying notions are used in the sequel, see in [39] for more information. Given is a surface ⊂ R 3 with parametric imageˆ ⊂ R 2 . Given also a parametrization of that surface S :ˆ → which is well-defined almost everywhere (a.e.), that is, every parametric location (θ 1 , θ 2 ) ∈ˆ is mapped onto a unique Cartesian location X = (X 1 , X 2 , X 3 ) through map S a.e. inˆ (see "Non-uniform rational b-spline surfaces" section for the NURBS parametrization of a surface). Accordingly, a covariant basis may be constructed as follows, where (•) ,α = ∂(•)/∂θ α andj = A 1 × A 2 2 . Map S is then well-defined at parametric locations wherej = 0. The components of the metric tensor A = A αβ A α ⊗ A β (also known as the first fundamental form of a surface) are given by, The contravariant components of the metric coefficient tensor, namely, A βγ can be obtained by the relation A αβ A βγ = δ γ α , where δ γ α stands for the Kronecker delta symbol, that is, δ γ α = 1 for α = γ and δ γ α = 0 otherwise. The Einstein's summation convention over repeated indices is assumed in the sequel. In this way, a contravariant basis can be constructed using the contravariant metric coefficients, where surface normal vector A 3 stays the same in both the covariant and the contravariant bases. The components of the curvature tensor B = B αβ A α ⊗A β (also known as the second fundamental form of a surface) are given by, which are linked to the curvature along the parametric directions θ α .

Mechanics of lightweight structures
Lightweight structures are typically represented by their mid-surface which consists of all particles X in the reference configuration, see for example in [40]. Such structures comprise membranes and shells which are considered thin, that is,h/R 20,h andR being the structural thickness and the radius of curvature, respectively. Herein a Lagrangian description of the motion is assumed and the problem is posed on the unknown displacement field d : → t of the mid-surface, where t stands for the current configuration consisting of all particles x = X + d at time t ∈ T where T = [0, T ∞ ], T 0 and T ∞ being the start and the end time of the dynamic process. In this way, assumed is that and t are represented by a parametric domainˆ via the geometric maps ("Differential geometry of surfaces" section) S and S t , respectively, see Fig. 2. Accordingly, the displacement field may be expressed on both the Cartesian basis e i and a curvilinear basis A α , A 3 (see Eq. (8)) as follows, The weak form of dynamic equilibrium for these structures can be written as follows: Find d ∈ H α ( ) for each time instance t ∈ T such that, where H α ( ) stands for the space of all square integrable vector-valued functions with square integrable derivatives up to α-th order in . Moreover, α = 1 and α = 2 for the membrane and the Kirchhoff-Love shell problem, respectively. This is because the curvature tensor involves second derivatives on the displacement field in Kirchhoff-Love shell analysis, see also in [13]. The first and the second terms in Eq. (13) stand for the inertia and damping of the structure, where ρ and c stand for the structural density and the damping coefficient, respectively. The form a is specialized for the membrane and the Kirchhoff-Love shell structural analysis in the following sections. The linear functional l : H α ( ) → R is defined as, where b stands for the body forces acting in . Especially for the membrane and the Kirchhoff-Love shell problems more types of external loads can be considered, see in [17,41] for more information. The inner product •, • 0, in the L 2 ( ) space (space of square integrable vector-valued functions in ) in Eq. (14) is defined as, Note that the weak form of dynamic equilibrium in Eq. (13) is formulated at each patch (i) independently, not accounting for the Dirichlet boundary conditions along a portion of the domain's boundary d ⊂ ∂ . The continuity across the multipatches and the weak application of the Dirichlet boundary conditions are specialized for the membrane and the Kirchhoff-Love shell in the following sections. The weak enforcement of the these constraints is essential as the multiple patches are not conforming along their common interfaces and the Dirichlet boundary conditions are typically enforced along trimming curves where the basis functions are not interpolatory. Thus, the strong enforcement of the interface and boundary constraints is in general inapplicable within IBRA. In the sequel, the dynamic form of weak equilibrium in Eq. (13) is posed on the decomposed open domain d defined in Eq. (7b) and the interface continuity conditions are discussed in the sequel.

Membrane structural analysis on multipatches
The isogeometric membrane structural analysis on multipatches employed in this work is based on the Penalty and Nitsche-type formulations presented in [41], where also weak application of the Dirichlet boundary conditions is considered. The Nitsche-type formulation is considered as a consistent extension of the Penalty method. This is because the Nitsche-type formulation in its original forms lacks coercivity and Penalty-like stabilization terms are added to restore coercivity. The corresponding stabilization parameters can be estimated by the solution of interface and boundary eigenvalue problems at each time step [16,41]. On the other hand, one obtains a pure Penalty formulation when leaving only the Penalty-like stabilization terms by excluding the additional Nitsche terms. Therefore the statement of the Nitsche-type formulation includes that of the Penalty formulation and thus both are herein presented in a unified manner. The presented numerical examples of multipatch isogeometric membrane structural analysis using the Penalty method are computed using the IBRA implementation in Carat++ in-house software [42] whereas the ones using the Nitsche-type method are computed using a MATLAB ® based framework freely available in [43]. Three-dimensional membranes can not in principle withstand compression without any form of stabilization due to wrinkling which is a type of zero energy mode. Wrinkling enhanced models have been extensively studied in the literature, see also in [44]. Additionally, membranes typically need to be under prestress in order to avoid wrinkling and be rendered stable. The latter results in a non-trivial design in that not every free-form shape may render a shape of static equilibrium. For this purpose form-finding methods have been developed [45] and in particular the Updated Reference Strategy, (URS) see in [46,47]. In this study, membranes in their original design are considered and moreover no form-finding is used for the herein presented numerical examples as the chosen geometries are by construction compatible with the applied prestress while no cables are embedded, see also in [21,48] for more information. The Green-Lagrange (GL) strain tensor of the mid-surface ε = ε αβ A α ⊗ A β is employed and its components are given by, a αβ = a α · a β being the covariant metric coefficients and a α the base vectors of the current configuration. Moreover, A α stand for the contravariant base vectors of the reference configuration. The components of the energetically conjugate 2nd Piola-Kirchhoff (PK2) stress-resultant force tensor n = n αβ A α ⊗ A β of the mid-surface are defined by means of the linear Hooke's law (Saint-Venant material), namely, where the components of the material tensor C = C αβγ δ A α ⊗ A β ⊗ A γ ⊗ A δ are given by, E and ν being the Young's (elastic) modulus and the Poisson's ratio, respectively. The traction along any curve γ on surface is defined by [38], where u α stand for the covariant components of the curve's γ normal vector u on surface and where n αβ 0 stand for the contravariant coefficients of the prestress tensor n 0 . Concerning the multipatch formulation, the solution space for the Nitsche-type and the Penalty methods is respectively. Fields restricted in a patch and along an interface are represented in the sequel by a superscript that is, • | (i) = • (i) and • | γ (i,j) = • (i,j) , respectively. Letχ stand for the interface displacement jump, that is, (j) . The mean interface traction field is given by, Accordingly, form a : V × V → R in Eq. (13) is defined as follows for the multipatch isogeometric membrane analysis using the Nitsche-type method, whereα : γ i → R andᾱ : d → R stand for the stabilization parameters (Nitschetype method) or the Penalty parameters (Penalty method) when the additional terms stemming from the Nitsche-type method are omitted. In case the Nitsche-type method is employed, the corresponding stabilization parameters are estimated automatically by solving a sequence of interface and boundary eigenvalue problems, see in [41] for more information. There are defined as piecewise constant along each interface γ (i,j) i and each Dirichlet boundary (i) d . On the other hand, in case the Penalty method is employed the Penalty parameters are discretization-dependent and are computed similar to the rule proposed in [17], namely, where h (i,j) −1 and h (i) −1 stand for the inverse of the smallest knot span length in the physical space along the interface trimming curve γ (i,j) i and along the trimming curve (i) d defining the portion of the Dirichlet boundary having an intersection with ∂ (i) within an isogeometric discretization, respectively. The norm of the material tensor in Eqs. (22) is , that is, the square root of the sum of its squared components.

Kirchhoff-Love structural analysis on multipatches
Similar to "Membrane structural analysis on multipatches" section, the herein employed isogeometric Kirchhoff-Love shell structural analysis on multipatches accounting for weak Dirichlet boundary conditions is based on a Penalty formulation as presented in [17,20]. The employed numerical example of multipatch isogeometric Kirchhoff-Love shell structural analysis using the Penalty method, that of the NREL phase VI wind turbine [23], is computed using the IBRA implementation within the in-house software Carat++. Moreover, small strains are herein assumed and thus the corresponding linearised theory is briefly presented.
The linearised GL strain strain tensors ε = ε αβ A α ⊗ A β and κ = κ αβ A α ⊗ A β for the membrane and the bending strain are defined as [49], where (•) ,αβ = ∂(•) ,α /∂θ β . The PK2 stress-resultant force tensor for the in-plane stiffness of the Kirchhoff-Love shell is defined as in Eq. (17). Similarly, the PK2 stress-resultant tensor for bending stiffness of the Kirchhoff-Love shell m = m αβ A α ⊗ A β is defined using also the linear Hooke's law, that is, The rotation field ω = ω ζ A ζ needs to be in this case defined, namely, where αζ is the Levi-Civita symbol. For the multipatch formulation using the Penalty method, the solution space is in this case Letχ stand for the jump on the rotation field across the multipatches, that is, (j) . In this way, the form a : V × V → R in Eq. (13) for the multipatch isogeometric Kirchhoff-Love shell analysis using the Penalty method is defined as, Additionally,α : γ i → R stands for the Penalty parameter associated with the imposition of the rotation continuity across the interfaces. It is chosen also piecewise constant and is defined similar to Eqs. 22, that is, along interface boundary γ (i,j) i .

Isogeometric spatial discretization on trimmed multipatches
Concerning the discretization of the aforementioned weak forms, the Isogeometric B-Rep Analysis (IBRA) is employed, see also in [20]. In this way, the finite dimensional subspace V h ⊂ V is constructed using the parametric description of each patch (i) as R( (i) ) being the space of all vector-valued piecewise rational polynomials for which the NURBS basis functions of the geometric parametrization constitute a basis in each patch (i) . Letφ Herein, the vector-valued NURBS basis functions are constructed as, where k = r/3 and l = r−3 r/3 +3 for all r = 1, . . . dim V (i) h stand for the indices of the control points and the Cartesian directions, respectively. Additionally, R In this way, projection of variational problem in Eq. (13) onto V h results into the following discretized in space equation system, whered,ḋ andd stand for the vectors of acceleration, velocity and displacement DOFs, respectively. In addition, M and D stand for the mass and damping matrices resulting from the spatial discretization of the first and second terms of variational problem in Eq. (13), respectively. Moreover, R(d) stands for the steady-state residual vector whose linearization results in the steady-state tangent stiffness matrix K(d) and whose entries are given by K ij (d) = ∂R i (d)/∂d j , R i (d) andd j being the i-th component of the residual vector and the j-th DOF, respectively. The definition of the tangent stiffness matrices for the membrane BVP can be found in [21,41] and for the Kirchhoff-Love shell BVP in [17,20]. In this study, the damping matrix is approximated using the Rayleigh damping method, that is, where α r and β r stand for the so-called Rayleigh damping parameters andd 0 stands for the initial condition on the displacement field, see in [50] for more information.

Time discretization and modal analysis
The Newmark method [51] is used in this study for the time discretization of linear equation system in Eq. (31). Accordingly, the continuous time domain T is discretized into a set of time steps tn. The system is linearised using the Newton-Raphson iterative method, that is, where îdn =dn ,î+1 −dn ,î ,dn ,î being the vector of DOFs at then-th time step and at i-th Newton-Raphson iteration. The dynamic stiffness matrixKn ,î and residual vectorRn ,î at then-th time step and atî-th Newton-Raphson iteration are defined by means of the corresponding steady state tangent stiffness matrix Kn ,î and residual vector Rn ,î , namely [35,41,52], where β n and γ n stand for the Newmark parameters and wheredn −1,î ,ḋn −1,î anddn −1,î stand for the displacement, velocity and acceleration DOFs at time step tn −1 , see in [41] for more information on the discrete equation systems. Concerning modal analysis, this is performed on the linearised system using the linear stiffness matrix K(d 0 ) and by solving the following eigenvalue problem, where ω i = 2πf i are the circular eigenfrequencies and where f i stand for the natural eigenfrequencies of the system.

Isogeometric B-Rep analysis of the NREL phase VI wind turbine
In this section, the NREL phase VI wind turbine with flexible blades [23] is employed as demonstration of isogeometric analysis on multipatch surfaces in industrial scale applications, see Fig. 3c. This numerical example is herein employed for the demonstration of IBRA on a real-world engineering structure in multiphysics environment and for validating the underlying computational models which are later on used in the context of FSI with the proposed isogeometric B-Rep mortar-based mapping method. A picture of the actual turbine can be seen in Fig. 3a. The corresponding CAD model consisting of rigid parts and the two flexible blades whose stiffness is enhanced using two longitudinal spars along the longitudinal trimming curves on the blades' surfaces, is shown in Fig. 3b. The problem is solved using the linearised Kirchhoff-Love shell theory within IBRA presented in "Isogeometric lightweight structural analysis on trimmed multipatches" section. The results of this simulation were firstly presented in the dissertation [35] and are repeated herein for the sake of completeness of this study. The original computational model in [53] involves a composite material model with varying thickness by analysing the data provided in [23]. Herein a simplified model is used with a Saint-Venant Kirchhoff material. The homogenized Young's modulus, density and thickness of the flexible blades are obtained by a calibration using a geometrically linear static and a modal analysis against the maximum displacement and the first eigenfrequency, respectively, computed in [53]. In this way, the Young's modulus, the density and the thickness of the flexible blades assumed to be E = 6 × 10 10 Pa, ρ = 1.515 × 10 3 Kg/m 3 andh = 7 mm, respectively. The Poisson ratio is then chosen as ν = 0.2. Regarding the static analysis, the flexible blades are subject to self weight, namely, b = −ρh e 3 . The results of IBRA for this example are compared with the results obtained using a standard finite element discretization of the flexible blades, see Fig. 3. Accordingy, the FEM model consists of 48630 triangular elements (Fig. 3(c)) based on a shell model with Reissner-Mindlin (RM) kinematics within Carat++ software ( [42]). Then, the corresponding h-refined multipatch NURBS computational model of the flexible blades is shown in Fig. 3(d). Subsequently, the NURBS computational model of the right blade is shown both intact and decomposed into its underlying trimmed patches in Fig. 4 where the geometric complexity and the large number of the underlying trimmed NURBS multipatches comprising the geometry is highlighted. It is worth mentioning that the spars and the tip of the NURBS computational model are connected to the rest of the blades' skin with a C 0 -parametric continuity forming geometric kinks, thus adding another complexity to the NURBS multipatch model. Each blade consists of 37 trimmed patches with 170 interface boundaries of highly diverse sizes and parametrizations. The scaling associated to the Penalty parameters is then chosen as the inverse of the minimum element edge size along each interface and Dirichlet boundary (see "Kirchhoff-Love structural analysis on multipatches" section) forα,α andᾱ, respectively.
The contour of the 2-norm of the displacement field d 2 across the blades in the current configuration due to self-weight for both the standard finite element analysis and IBRA is shown in Fig. 5 demonstrating excellent accordance of the results. Moreover, an eigenfrequency analysis for both models is performed, see Eq. (35), and the first three eigenfrequencies of both the standard FEM and IGA models are shown in Fig. 6 demonstrating once more an excellent accordance of the results also in this context.

Isogeometric B-Rep mortar-based mapping method on trimmed multipatches
In this section the isogeometric B-Rep mortar-based mapping method is described in which fields are transformed between a low order faceted discretization and a NURBS multipatch description of a surface. Additionally, the NURBS multipatch description of the surface can be used as a mediator Exact Coupling Layer (ECL). The ECL is used to smooth the transformed fields between two low order representations of the FSI interface and it is represented using the exact CAD model of the common interface. In the sequel of this chapter it is assumed that and h are the exact surface representations stemming from CAD as described in "Computer-aided design using non-uniform rational b-splines" section and a low order faceted representation of the surface with a finite number of polygonal elements, respectively.

Theory
This section presents the problem placement along with the weak formulation and the discrete equation system governing the isogeometric B-Rep mortar-based mapping method which is formulated as extension of the isogeometric mortar-based mapping method proposed in [52]. This isogeometric B-Rep mortar-based mapping method was firstly presented in [35] and it is herein repeated in more detail by including additional implementation and methodological aspects. All formulas are provided for the special case where fields are transformed between a low order discretized and a multipatch NURBS surface, however the following principles might well apply for any mortar-based mapping method. Accordingly, let T i , i = 1, . . . , n e ∈ N stand for the set of standard low order finite elements in h . Let q h ∈ V h q be a field defined isoparametrically on the low order discretized surface h where, and where P α (T i ) stands for the linear (α = 1) or bilinear (α = 2) basis functions in each finite element T i . Let also (i) , i = 1, . . . , n s be a non-overlapping decomposition of as defined in Eqs. (7). The goal is to find that field q ∈ L 2 ( d ) which is the closest to q h in the L 2 ( d )-space, namely, where field q is discontinuous along the interface γ i and where d is defined in Eq. (7b). The problem in Eq. (37) is herein also subject to the following interface and boundary Problem placement: Low order FEM discretization and CAD representation of a surface [35] conditions, where q (i) = q | (i) as described in 'Membrane structural analysis on multipatches" section and ω (i) t = ω t (q (i) ) using only the rotation around the tangent to each boundary γ (i,j) i vector following the definition introduced in Eq. (25), that is, q γ and q 3 being the curvilinear components of q defined similar to Eq. (12) and u α the contravariant components of the normal vector u along γ i , see in [17] for more information. The aforementioned problem is depicted in Fig. 7. Similar to "Kirchhoff-Love structural analysis on multipatches" section, letχ andχ represent the interface jump on q and its rotation around the tangent to the interface vector ω t , respectively. With the aforementioned condition one can restrict the transformed field along d while simultaneously maintaining a solution in H 1 (γ i ). The solution of problem (37) subject to the interface and boundary conditions (38) can be obtained by the minimization of the following, augmented with Penalty terms, variational formulation: Given a q h ∈ V h q , find a q ∈ V q , such that, for all δq ∈ V q . Space V q is defined in a similar manner as V h in "Isogeometric spatial discretization on trimmed multipatches" section, that is, The weak form in Eq. (40) has a unique solution according to Lax-Milgram theorem [54] since the bilinear form defining the left-hand side of Eq. (40) is coercive and continuous in V q × V q whereas the functional defining the right-hand side of Eq. (40) is linear in V q . However, the quality of the solution depends on the choice of the Penalty parameters as typical for the Penalty methods, see also in [55]. Since spaces V h q and V q are by construction finite dimensional, given the vector-valued standard finite element and NURBS basis functions at each patch (i) (see also "Isogeometric spatial discretization on trimmed multipatches" section), j stand for the DOFs of the finite element discretization and the DOFs of the isogeometric discretization within each patch (i) , respectively. These can be grouped into vectors, Accordingly, the discrete equation system corresponding to the weak form in Eq. (40) reads, where, Additionally, the entries of matrices C α are given by, and where the ± sign in Eq. (46f) depends on the ordering of the neighboring patches (i) and (j) in I. The matrix containing the Penalty terms, namely, Cα ,α,ᾱ in Eq. (45c), is optional and can be selectively used. In case some of the corresponding Penalty contributions are not considered, the corresponding Penalty parameter at the subscript of Cα ,α,ᾱ is replaced by zero. Its application depends on the numerical example, see "Numerical results" section. As aforementioned, problem in Eq. (44) is well defined provided that the corresponding Penalty parameters are carefully chosen. Herein, the choice of the Penalty parameters is made similar to Eqs. (22), that is, where the definitions of h −1 can be found in "Membrane structural analysis on multipatches" section. Moreover, B t stands for the B-operator vector B t resulting from the discretization of the rotation field ω t = ω ·ê n (see also in Eq. (25)),ê n being the outward normal vector to the boundary (yet tangent on the surface) along which the bending component ω t of the rotation ω is defined. In other words, the i-th component h ⊂ h which has a projection on (i) [35] of B t is given by, e α n being the contravariant components of unit vectorê n , for more information see in [35]. For the mortar-based transformation of a field q defined isogeometrically over a multipatch NURBS surface to a field q h defined on a standard finite element discretized surface, problem in Eq. (44) simply reverses, that is, where C nn is defined similar to C rr in Eq. (45a) with entries, and where C nr = (C rn ) T , see in Eq. (45b). The fact that the variational problem of the isogeometric B-Rep mortar-based mapping method is well-posed is reflected onto the fact that matrices C rr and C nn , are square, symmetric and positive definite, see Eqs. (44) and (49), respectively.

Realization
In this section the implementation and methodological aspects of the isogeometric B-Rep mortar-based mapping method are discussed in detail. As already mentioned in "Theory" section, fields are to be transformed between a trimmed multipatch NURBS and a low order discretized surface. In the following, the corresponding algorithms and their properties are discussed in detail.

Projection of the finite element mesh on multipatch surface
Firstly the numerical evaluation of the integrals on (i) is discussed, see Eqs. (45a), (45b) and (50), respectively. Within this study, the exact geometry (i) for each patch is chosen as the integration surface. Accordingly, the finite element mesh is projected onto the NURBS surface, by projecting each node X i , i = 1, . . . , n n ∈ N onto (i) through the nonlinear

Fig. 9
Realization: Projection of the elements of (i) h , e.g. X j -X k -X l (see Fig. 8), onto (i) [35] map θ (i) j = (S (i) ) −1 X j for all nodes in the finite element mesh using a Newton-Raphson scheme. In order to accelerate the projection process, the Newton-Raphson algorithm for the projection of a node onto each patch (i) is only performed for these nodes X i which are contained within the patch bounding box scaled by a small tolerance. The latter scaling of the patch bounding box is necessary in order to include as candidates these nodes which have projection very close to the patch boundary. Consider the finite element X j -X k -X l in right part of Fig. 8. Each of its nodes are projected onto (i) to obtain the corresponding parametric coordinates in the parametric spaceˆ (i) of the patch. Subsequently, a linear connection in the parametric spaceˆ (i) is made in order to obtain the image of the element in the parametric space of the patch, namely, θ Fig. 9. The latter linear approximation of the element edges in the parametric space of the patch is a consistent approximation since the projection error tends to zero as the element size gets smaller. Then, the image of the element back onto the geometric space can be obtained using the geometric transformation X (i) = S (i) (θ (i) ), defined in Eq. (6), see right part of Fig. 9.

Projection on patch boundary
Consider one NURBS surface patch and the corresponding part (i) h ⊂ h of the finite element mesh h which has a projection on (i) , see Fig. 8. Occasionally, there exist finite elements which are partially projected inside and partially outside the boundaries of the patch parametric space (edge X n -X m in Fig. 10a). To obtain the parts of these elements which lie within the computational domain of the patch, the corresponding finite element edges are clipped with the patch boundary using two methods: A Newtontype and a bisection method. Firstly the Newton-type method is employed due to its quadratic convergence behavior and in case convergence to the solution is not achieved, then the bisection method is used which is slower in terms of convergence but robust, providing always a solution.
For the Newton-type method, as closest point of the patch boundary ∂ (i) to edge X n -X m is considered the pointX nm ) which lies on the plane spanned by the surface unit normal of (i) on the patch boundary and edge X n -X m , which in addition passes from vertex X n . The surface normal vector defining the aforementioned plane is given by, 2 ) along boundary ∂ (i) since patch boundaries are always aligned with either parametric direction. In this way, the residual equation writes, For the case depicted in Fig. 10a it clearly holds ξ (i) = θ (i) 1 (see Fig. 9). The Newton-Raphson linearisation of residual in Eq. (52) results in, where ξ nm,j and where Jacobian ∂r (i) /∂ξ (i) is given by, where A . Solving the iterative system in Eq. (53) yields the solutionX . The corresponding convergence criterion is based on the residual magnitude, namely |r (i) (ξ (i) )|. To obtain the part of the edge X n -X m which is projected on patch (i) , namely X n -X (i) nm , one has, thus obtaining the closest point of edge X n -X m to patch boundary ∂ (i) namely X nm − X n 2 / X m − X n 2 . As initial guess for the nonlinear problem in Eq. (53) the middle parametric location in either parametric direction θ (i) 1 or θ (i) 2 is in this study chosen depending on which parametric direction the patch boundary is aligned. The Newton-type algorithm is depicted in Fig. 10b.
Concerning the bisection method, two sequences of points X i,j and X o,j on edge X n -X m are generated, where X i,j have a projection and where X o,j do not have a projection on patch (i) . We seek the pointX (i) nm which the closest to edge X n -X m on patch boundary ∂ (i) and its image X (i) nm on edge X n -X m . The initial condition of the bisection algorithm are the vertices of the edge itself which is X i,0 = X n and X o,0 = X m for the case depicted in Fig. 10a. Then, the midpoint of segment X i,j -X o,j is iteratively computed as X m,j = (X i,j + X i,j )/2 and assigned to either X i,j+1 or X o,j+1 depending on whether or not it has a projection on patch (i) while its projection X nm,j 2 . In case the obtained projection X nm is not sufficiently close to the boundary then an additional projection step is performed whereθ nm is projected on the a Finite element edge crossing patch's boundary.
Newton-type algorithm for projecting an edge onto a patch boundary.
c Bisection algorithm for projecting an edge onto a patch boundary. Fig. 10 Realization: Projection of a finite element edge onto patch's boundary using the bisection and the Newton-type method closest patch boundary via a two-dimensional point-to-curve projection in the parametric space of patch (i) , namelyˆ (i) . Accordingly, the ratio of the projected part of edge X n -X m on patch boundary ∂ (i) can be obtained as λ nm − X n 2 / X m − X n 2 similar to the Newton-Raphson method described before. The bisection algorithm is depicted in Fig. 10c.
The Newton-type and bisection methods developed for the projection of a finite element edge onto a patch boundary may not yield the exact same results because of the inherent non-convex nature of the closest projection. This is because the projection problem might be locally non-convex and thus a unique solution may not be available. However, this algorithmic step is only necessary for obtaining a partial image of the element in the parametric spaces of the different patches it has a projection for evaluating the integrals presented in "Theory" section. The combination of both aforementioned methods for the projection of an element edge onto a patch boundary provides a very reliable and robust solution for the projection of an edge on a patch boundary.

Reconstruction of the finite element in the parametric space of the patch
Having the projection of the element edges on the patch boundary ("Projection on patch boundary" section) the next step is to reconstruct the partially projected finite element in the parametric space of the patch. This is necessary as the basis functions of the element in the finite element mesh are computed based on the shape of the element in the parametric space of the patch for the evaluation of the integrals in Eqs. (46b) and (50) as the integration in this study is always performed over the NURBS parametric description of the interface. Since herein finite element meshes with including elements of triangular and quadrilateral types are employed, this study is confined into the three cases depicted in Fig. 11, without loss of generality, while all other cases may be treated by a combination of the latter three cases.
For case 1 (Fig. 11a) consider a triangular element X m -X n -X l for which nodes X n and X l have a projection on patch (i) and where node X m does not have a projection on patch (i) . Given the parametric coordinates θ (i) k and θ (i) r of the projections of edges X n -X m and X l -X m with the patch boundary ∂ (i) , respectively, these are extended in the fictitious part of parametric spaceˆ (i) and their intersectionθ (i) m is considered as the reconstructed fictitious image of node X m in the parametric space of patch (i) .
For case 2 (Fig. 11b) consider a triangular element X m -X n -X l for which node X m has a projection on patch (i) and where nodes X n and X l do not have a projection on patch (i) . The fictitious parametric imagesθ (i) n andθ (i) l of nodes X n and X l , respectively, are reconstructed using the projections θ (i) mn and θ (i) ml of edges X m -X n and X m -X l with patch boundary ∂ (i) , respectively, and by also using the corresponding projection ratios λ mn = 1 − λ nm and λ ml = 1 − λ lm . The computation of projection ratios λ nm and λ lm is detailed in "Projection of the finite element mesh on multipatch surface" section.
For case 3 (Fig. 11c) consider a quadrilateral element X m -X n -X p -X l for which node X m has a projection on patch (i) and where nodes X n , X p and X l do not have a projection on patch (i) . The fictitious parametric coordinatesθ (i) n andθ (i) l of nodes X n and X m are computed similar to case 2. The fictitious parametric coordinateθ (i) p is obtained by means of the projection θ (i) mp of edge X m -X p with the patch boundary ∂ (i) and by also using projection ratio λ pm = 1 − λ mp .

Numerical integration
The set of trimming curves γ (i) k , k = 1, . . . , n (i) t ∈ N which trim patch (i) is subsequently linearised, see also in [35]. For the linearisation, the union ofp k , see in "Computer-aided design using non-uniform rational b-splines" section. Next, the projected elements onto the NURBS parametric space are clipped with the a Case 1: A node without projection on patch Ω (i) reconstructed using two neighboring nodes which have a projection on Ω (i) (triangular element).
Case 2: Two nodes without projection on patch Ω (i) reconstructed using the same neighboring node which has a projection on patch Ω (i) (triangular element).
c Case 3: Three nodes without projection on patch Ω (i) reconstructed using the same neighboring node which has a projection on patch Ω (i) (quadrilateral element) Fig. 11 Realization: Reconstruction of nodes without projection on patch (i) with fictitious parametric coordinatesθ using neighboring nodes from the same finite element with projection on patch (i) aforementioned linearised trimming curves, in order to exclude parts of the elements which are projected outside the computational domain of the patch. In this way, the computational domain for the evaluation of the integrals in Eq. (40) is obtained, see the shaded area in Fig. 9. Then, the projected finite elements in the parametric space of patch (i) are clipped with the knot lines of the parametric space of patchˆ (i) , thus obtaining subregions where the integrands in Eqs. (46) are C ∞ -continuous and where the Gauss integration can be performed. The aforementioned resulting regions may attain an arbitrary polygonal shape and thus they are subsequently triangulated using a very simple triangulation rule since this is merely needed for the distribution of the Gauss points. Then, the Gauss points are generated on each resulting sub-triangle where the integrals in Eqs. (46) are numerically evaluated.

Fluid-structure interaction
In this section the employed approach to FSI is demonstrated. Accordingly the equations governing the incompressible Navier-Stokes equations for the Computational Fluid Dynamics (CFD) problem on a moving frame are presented, as a body fitted approach is herein utilized and subsequently the staggered (partitioned) Gauss-Seidel (GS) FSI approach is provided.

Computational fluid dynamics
In this study, the incompressible Navier-Stokes equations for a Newtonian fluid are numerically solved by means of the Finite Volume Method (FVM), see in [3], within the opensource software OpenFOAM ® , see also in [56]. The underlying equations are common in the literature, see for instance in [57][58][59][60][61] among others and therefore are not herein repeated. In the sequel, a set of variables which are necessary for the comprehensive presentation of the study are introduced. Accordingly, letṼ be the computational fluid domain. The primary unknown fields of the incompressible Navier-Stokes BVP are the velocity and the pressure fields, u and p, respectively. Prescribed fluid inlet velocity u =υ is assumed along a portion of its boundary˜ d ⊂ ∂Ṽ whereas applied tractionst =τ are assumed along another portion of the domain's boundary˜ n ⊂ ∂Ṽ . The incompressible Navier-Stokes BVP is posed in a moving frame in terms of a Arbitrary Lagrangian-Eulerian (ALE) description of the fluid motion to accommodate the mesh displacements along the FSI interface. Among the various solution procedures and adaptations offered by OpenFOAM ® , herein a laminar solver for the cavity FSI benchmark, a Large Eddy Simulation (LES) for the hangar FSI simulation and an Unsteady Reynolds Averaged Navier-Stokes (uRANS) for the NREL phase VI wind turbine FSI simulation are employed, see in [3,27,62] for more information. These diverse solution approaches are chosen in order to show the applicability of the proposed FSI methodology for the various fidelities of the CFD problem.

Partitioned fluid-structure interaction
In this section the partitioned FSI approach is briefly introduced, see also in [26]. Herein assumed is that the structural and the fluid domains share a common interface S = ∂Ṽ ∩ , where is the domain of definition for the structural IBVP, see "Kirchhoff-Love structural analysis on multipatches" and "Membrane structural analysis on multipatches" sections for the Kirchhoff-Love and the membrane IBVP, respectively. In this way, the structural and the fluid IBVPs are subject to the following Dirichlet and Neumann interface conditions across S , respectively, to account for a continuous solution across . Accordingly, the Computational Structural Dynamics (CSD) problem governed by either the Kirchhoff-Love or the membrane theory, see "Kirchhoff-Love structural analysis on multipatches" and "Membrane structural analysis on multipatches" sections , respectively, is discretized using IBRA or standard FEM. Additionally,ḋ and u stand for the structural and fluid velocity fields, respectively, which have to be equal across the common interface as per Eq. (56a). Moreover, Eq. (56a) enforces also the continuity of the interface displacements across the common interface, namely, when integrating Eq. (56a) on time, U being the fluid displacement field across the FSI interface. Concerning the traction equilibrium across S , see Eq. (56b), the surface traction vectorb of the structural FSI interface contributes to the body force vector b on the right-hand side of structural weak forms in Eqs. (13) and the fluid traction vectort (see "Computational fluid dynamics" section).
Since the FVM scheme within OpenFOAM ® is chosen for the solution of the CFD problem, the fluid FSI interface is represented by a low order faceted surface, see "Computational fluid dynamics" section. Additionally, the herein presented FSI simulations involving IGA discretizations for the structure or the ECL, employ the isogeometric B-Rep mortar-based mapping method introduced in "Isogeometric B-Rep mortar-based mapping method on trimmed multipatches" section and are compared against FSI simulations of the same problems involving standard FEM discretizations of the structure using the standard mortar-based mapping method, see in [28]. Let S and S h denote in the sequel the exact CAD representation and a low order discretization of the FSI interface. The restriction of the FSI interface at each patch is denoted by S (i) = S ∩ (i) . Moreover, S h may represent the FSI interface of the fluid domain since the FVM is employed as discretization method for the CFD problem or the FSI interface of the structural domain whenever FEM is used for the CSD problem. A distinction should be then made clear from the context.
The solution of the CSD and CFD problems subject to interface conditions in Eqs. (56) and (57) is achieved using a fixed-point iteration approach known as Gauss-Seidel (GS) iterative method [26]. In this study, the underlying framework for the utilization of the aforementioned partitioned GS approach is hosted in EMPIRE software [63]. Since this coupling scheme has been extensively explained in the literature [30,35,52,53,64] it is herein presented only an outline of the general algorithm along with the variables necessary for the sequel of this study. Accordingly, the CFD problem is solved as a Dirichlet problem by complying with the resulting interface displacements from the solution of the CSD problem, whereas the CSD problem is in turn solved as a Neumann problem subject to the interface traction field emanating from the solution of the CFD problem. In this way, the displacement field of the structural FSI interface d | S is transformed onto the fluid FSI interface displacement field U and accordingly the fluid traction field at the FSI interfacet | S h is transformed onto the traction field on the structural FSI interface. This interaction takes place at each time step tn, assuming a matching time discretization for both the CSD and the CFD problems, until a specified termination criterion based on the relative change of the structural displacement across the FSI interface in the 2-norm at each FSI iterationk is met, given a user-defined tolerance˜ , namely, dn ,k | S being the vector of structural displacement DOFs on the FSI interface at time step tn and at k-th FSI iteration. Moreover,dn −1 | S stands for the vector of structural displacement DOFs on the FSI interface at the converged coupled FSI state and at time step tn −1 . In Eq. (58) the indexî introduced in "Time discretization and modal analysis" section ond representing the Newton-Raphson iteration is omitted sincedn ,k is assumed to be the set of displacement DOFs at the converged state of the CSD problem for the geometrically nonlinear analysis at time step tn and atk-th FSI iteration. The aforementioned interface fixed point iterations are stabilized and accelerated using the Aitken relaxation method, see in [65]. Moreover, the displacement and the traction fields can be transformed using Eq. (44) or Eq. (49), depending on the transformation direction, which is known as consistent transformation, see also in [65]. In case Eq. (44) is used for the transformation of displacements, all additional Penalty terms defined in Eq. (45c) are considered, since the displacement and the rotation continuity across the multipatches along with the Dirichlet boundary conditions need in this case to be weakly enforced. On the other hand, only the additional Penalty matrix Cα ,0,0 is employed in case Eq. (44) is used for the transformation of the tractions, excluding continuity of the rotation field and weak application of the Dirichlet boundary conditions. LetF t andF b be the forces acting on the structural and fluid FSI interfaces, respectively. In this case, it is also possible to transform the force vector in a conservative manner, meaning that the discrete interface work is exactly satisfied. This can be achieved by [35], The relation in Eq. (59) above is also known as conservative mapping [28,65]. Similar to the consistent mapping, in case that the conservative transformation of forces takes place from the CAD surface to the low order discretized surface in the frame of the ECL approach, relation in Eq. (59) simply inverses, namely, In this study, the additional Penalty terms are excluded for the conservative transformation of the force vectors. Thus, a direct transformation of the consistent force vectors is offered by the conservative mapping as per Eqs. (59) and (60) bypassing the computation of the traction fields. A thorough comparison of both the consistent and the conservative transformation of tractions and forces, respectively, is provided in [28] and in practice either can be used.

Numerical results
In this section, one benchmark and two real-world application examples are presented. As benchmark the lid-driven cavity FSI benchmark is herein chosen as standard in literature, see in [66]. The real-world applications comprise the FSI simulations of the inflatable hangar structure presented in [41] and the NREL phase VI wind turbine with flexible blades (see "NREL phase VI wind turbine in numerical wind tunnel" section) in corresponding numerical wind tunnels. The results are validated using standard structural finite element discretizations and the standard mortar-based mapping method as presented in [28]. Moreover, for each numerical example, results of given time steps are selected and transformed on either the standard finite element or on the NURBS surface using the proposed B-Rep mortar-based mapping method.
The consistent mapping approach is chosen for the numerical examples of the liddriven cavity FSI benchmark and the NREL phase VI wind turbine with flexible blades, whereas the conservative mapping approach is used for the inflatable hangar in numerical wind tunnel numerical example. In order to quantify the residual energy that occurs using the consistent mapping approach (since the interface energy conservation is not by construction satisfied in this case), the interface work from both the fluid and the structural sides as well as the residual interface energies are used. The important role of the energy transfer within partitioned FSI simulations via the mapping employed method has been demonstrated in [67] and [68] among others. Therefore, the energy transfer using the isogeometric B-Rep mortar-based mapping method is also herein considered as an extension of study in [52] on real-world trimmed multipatch geometries. For the sake of completeness, the underlying formulas which are used for the evaluation of the energy transfer in [52] using the proposed isogeometric B-Rep mapping method are repeated herein. Let W S be the work done by the fluid forces on the moving FSI interface S at a given time. This is defined as follows, where t and d stand for the traction and displacement fields along the FSI interface S at a given time defined on either the fluid or the structural subdomain. The superscript • (f) or • (s) is accordingly used. Since in this study the traction and the displacement fields are transformed using the mortar-based mapping method from one interface to the other (see "Isogeometric B-Rep mortar-based mapping method on trimmed multipatches" section), their discrete representations (see Eq. (42)) can be substituted in Eq. (61). In the case when FVM and IGA are used for the discretization of the fluid and structural fields, the following expressions for the interface work are obtained: If the ECL is applied, the structure is discretized by FEM and thus W (s) S is computed using Eq. (62a) using the corresponding structural C nn matrix. The residual interface energy E S is defined by the difference between the structural and the fluid interface work:

Lid-driven cavity
In this section the lid-driven cavity FSI benchmark is employed (Fig. 12), see also in [66], for the demonstration and evaluation of the isogeometric B-Rep mortar-based mapping method and the application of isogeometric membrane analysis on multipatches in multiphysics problems, due to the simplicity of its geometry while relatively large deformations are allowed. This example was firstly presented in [35] and it is herein complemented with an additional study on the evolution of the interface work from each subdomain and their difference against the time. Accordingly, the CFD domain is a unit square with thickness of 10 cm. The CFD problem is solved as two-dimensional in X 1 -X 2 plane, meaning that the velocity and pressure are constant in the width direction. Accordingly, the CFD domain is discretized using a 30×30 grid for all employed simulations, see Fig. 12b. On the other hand, the membrane structure is discretized using a reference finite element mesh with 100 bilinear elements (FEM100), a coarse finite element mesh with six bilinear elements (FEM6), a single patch geometry with twenty quadratic in X 1 -direction and linear in the X 3 -direction elements (IGA1) and a trimmed two-patch geometry where the interface is an arc of a circle (IGA2) see Fig. 13. The CSD problem is also two-dimensional and accordingly the displacement field d 0 3 in the X 3 -direction is set to zero.
A set of FSI simulations is accordingly performed involving FEM100 mesh for the structure as reference (FSI-FEM100), IGA1 structural discretization (FSI-IGA1), IGA2 structural discretization (FSI-IGA2), FEM6 mesh for the structure (FSI-FEM6), FEM6 mesh for the structure with IGA1 parametrization for the ECL (FSI-FEM6-IGA1) and FEM6 mesh for the structure with IGA2 parametrization for the ECL (FSI-FEM6-IGA2). The streamlines of the FSI simulations at time t = 19 s using the IGA structural discretizations  against the reference solution using FEM100 is shown in the set of Fig. 14 Fig. 16 Lid-driven cavity: Time-displacement curves for the structural displacement at X 1 = 0.5 m [35] see Fig. 16. Accordingly, the time-displacement curves for the FSI simulations and the FSI simulations with the ECL are depicted in Figs. 16a and 16b, respectively. As it can be observed, the FSI simulations using single and multipatch IGA for the structural discretization deliver highly accurate results when compared to the reference FSI simulation involving a highly refined finite element structural discretization. Please note that the cavity FSI benchmark is a two-dimensional benchmark. To be able to demonstrate the applicability of the trimmed multipatches for this benchmark, the underlying patches are trimmed in the X 3 -direction, see Fig. 13. The latter leads in a non-uniform placement of the integration points at each patch about the trimming, triggering some additional dynamics of the membrane X 3 -direction that are not seen in the finite element and single patch discretizations, where the integration points are uniformly placed with respect to the X 3 -direction. The latter explains the slight deviation of the results obtained by the two-patch discretization in Fig. 16, which is nevertheless minimal.
On the other hand, the application of the ECL improves the quality of the interface displacement field, see Figs. 15b and 15c for the single and multipatch representation of the ECL, respectively, whereas it produces highly accurate results given that the CSD problem is only discretized using six elements, see Fig. 16b. It can be seen in Fig. 16b   that the response of the membrane in this case is slightly stiffer when using the ECL. The underlying reason is that the smoothing induced by the ECL adds a constraint to the coupling interface. This constraint renders the membrane in this case slightly stiffer, especially because the membrane is purposely discretized using only six elements to produce a rough displacement field and to highlight the effect of the smoothing by means of the ECL, see Fig. 15. The isogeometric mortar-based mapping method is subsequently evaluated based on the work transferred through the FSI interface, see Fig. 17. An important finding herein is that the evolution of the interface work from both the structural and the fluid sides is smoother when using the isogeometric mortar-based mapping method as opposed to the standard mortar-based mapping method, see the magnifying window in Fig. 17a. The rest of the patterns in the latter figure correspond to the interface displacement field that was observed previously in Fig. 16. In Fig. 17b it can be seen that the residual interface energy as defined in Eq. (63) remains in insignificant levels for all the underlying discretizations and for both mapping technologies (standard and isogeometric).
Next, the isogeometric B-Rep mortar-based mapping method is evaluated with regard to its convergence behaviour. In the following convergence graphs the observed convergence rates are mentioned, for which mathematical proofs are however pending. Accordingly, the quantities of interest are the displacement and the traction fields along with their corresponding transformations. The first set of graphs in Fig. 18 shows the convergence of the consistent mapping in the displacement and traction fields onto S h , when these are originally defined on the CAD surface.
More specifically, the convergence graph in Fig. 18a shows the convergence based on the relative error in the L 2 (S )-norm for the displacement field defined on the CAD surface for the FSI-IGA1 and the FSI-IGA2 simulations at time t = 19 s against its transformed field on h i , for various mesh densities with 5, 10, 20, 40, 80 and 160 elements, respectively, where h stands for the minimum finite element edge in h i . Similarly, Fig. 18b shows the convergence based on the relative error in the L 2 (S )-norm of the traction field defined on the CAD surface for the FSI-FEM6-IGA1 and FSI-FEM6-IGA2 simulations at time t = 19 Relative error in the L 2 (S )-norm for the traction transformation.

Fig. 18
Lid-driven cavity: Relative error in the mapping from IGA to FEM using a reference field defined on IGA1 and IGA2 surface representations at t = 19 s of the corresponding FSI simulation [35] Fig. 19 Lid-driven cavity: Mapping of the displacement field defined on IGA1 and IGA2 surface representations at t = 19 s from FSI-IGA1 and FSI-IGA2 simulations, respectively, to FEM for a given discretization level [35] s against its transformed field on h i . Both graphs show excellent convergence rates. The displacement and traction fields defined on the IGA1 and IGA2 surface representations from FSI-IGA1, FSI-IGA2 and FSI-IGA1-ECL, FSI-IGA2-ECL simulations at time t = 19 s along with their corresponding transformations on the low order discretized surface are then shown in Figs. 19 and 20 , respectively.
In this case, the square root of the smallest area among the isogeometric elements in the multipatch model is chosen as characteristic measure for the discretization density h Fig. 20 Lid-driven cavity: Mapping of the traction field defined on IGA1 and IGA2 surface representations at t = 19 s from FSI-IGA1-ECL and FSI-IGA2-ECL simulations, respectively, to FEM for a given discretization level [35]

Fig. 21
Lid-driven cavity: Error in the transformation of tractions from FEM to IGA using a reference traction field defined on the FE surface at t = 10 s of the FSI simulation [35] for the relative errors in the domain S , see Fig. 21a. Regarding the interface error, the smallest element edge length along the interface γ i is chosen as characteristic measure h of the mesh density, see Fig. 21b. The relative error in the L 2 (S )-norm for the traction field taken from the FSI-FEM100 simulation at time t = 19 s and transformed into the different CAD surface representations with the aforementioned refinement is shown in Fig. 21a, for both the low-and the high order bases. Accordingly, the jump in the traction field along the interface γ i for the multipatch CAD representations of the surface in the L 2 (γ i )norm is shown in Fig. 21b. It can be observed that the single patch solution demonstrates quadratic convergence order for both polynomial order settings, whereas the convergence order drops to linear concerning the multipatch model with Penalty. The latter is expected as the optimal convergence rates for the Penalty methods are typically bounded by both the mesh size and the underlying Penalty parameters themselves, in contrast to the single patch discretizations. For more information on the convergence rates please refer to [69]. The traction field defined on the fluid FSI interface from FSI-FEM100 simulations at time t = 19 s and its transformation into the various CAD representations of the surface are then shown in Fig. 22 for a qualitative assessment of the isogeometric B-Rep mortar-based mapping from a low order surface discretization to CAD surface representations.

Fig. 24
Lid-driven cavity: Transformation of the displacement field defined on the low order surface discretization at t = 19 s from FSI-FEM100 simulation to the various CAD surface representations of for given refinement levels [35] Lastly, convergence graphs for the isogeometric B-Rep mortar-based mapping of the displacement field defined on the fluid FSI interface mesh from simulation FSI-FEM100 at time t = 19 s onto IGA1 and IGA2 surface representations are drawn, corresponding to the ECL concept (Fig. 23). The refinement studies of IGA1 and IGA2 surface models are the same as previously, where herein the complete Penalty matrix Cα ,α,ᾱ is taken into account. The relative error on the displacement field in the L 2 (S )-norm is shown in Fig. 23a where as before quadratic rates of convergence are observed for the single patch model and linear convergence rates are observed for the multipatch model. The solution accuracy however is not significantly improved for the high order bases when compared to the low order bases. The interface error on the jump of the displacement field and its rotation around the tangent to the patch boundary vector for the different refinement levels in the L 2 (γ i )norm is shown in Fig. 23a regarding the trimmed two-patch model IGA2. One can observe here an improvement of the fulfillment of the interface conditions for the high order bases and an almost linear convergence rate. Then, Fig. 23b shows the error in the fulfillment of the Dirichlet condition in the L 2 ( d )-norm for both CAD models where quadratic convergence rates are observed. Additionally, an improvement of the solution for the high order bases is also observed in this case. In the latter case, the minimum element edge size along d is chosen as a measure of the discretization density h. Lastly, the 2norm of the displacement field defined on the fluid FSI interface and the corresponding 2-norms of its transformation on the IGA1 and IGA2 surface representations is shown in Fig. 24. Once more, an excellent transformation can be seen which in addition respects the interface and boundary conditions when using the isogeometric B-Rep mortar-based mapping method in combination with Penalty.

Inflatable hangar in numerical wind tunnel
In this section the FSI simulation of an inflatable hangar [41] in numerical wind tunnel is investigated. A similar FSI simulation is also presented in [35]. However, the stiffness of the hangar in terms of the applied prestress and inner pressure in that numerical investigation was chosen such that significant deformation occurs in order to highlight the advantages of the ECL approach. That led in many cases to wrinkling of the membrane, which is a type of rigid body mode, as discussed in "Membrane structural analysis on multipatches" section, rendering the results inaccurate. In the present study, the stiffness of the hangar in terms of the applied prestress and inner pressure is chosen higher such that  [35] no significant wrinkling occurs allowing for an appropriate comparison of the standard mortar-based mapping method and the newly proposed IBRA mortar-based mapping method. Accordingly, the material properties of the hangar are described in [41], with the only difference that herein three different magnitudes for the inner pressure are chosen, namely, b 2 = 10 3 Pa (p1000), b 2 = 2×10 3 Pa (p2000) and b 2 = 4 ×10 3 Pa (p4000) and the prestress is adapted so that the hangar remains in equilibrium with respect to its shape, see in [41] for more details. The wind is modelled through the incompressible Navier-Stokes ("Computational fluid dynamics" section), with kinematic viscosity of the airν = 1.5451 × 10 −5 m 2 /s. Regarding the CFD domain a mesh with approximately 100,000 polyhedral elements is employed, which is successively refined towards the hangar region.
The locally refined region around the hangar is made using the snappyHexMesh mesh generator of OpenFOAM ® , see Fig. 26. Accordingly, an LES solution approach is employed using a one equation eddy-viscosity Subgrid-Scale Model for the turbulence modelling, see in [70,71] for more details. The inlet velocityυ is chosen using a 1/7-power law from the bottom up to the height of the hangar and then is kept constant with an amplitude of 13 m/s. The behavior of the structural deformation due to wind and for the various internal pressure magnitudes with accordingly adjusted prestress is investigated. No-slip conditions are assumed at the two side and the bottom walls whereas slip conditions 3 are assumed at the top wall of the wind tunnel. Lastly, the pressure is set to zero at the outlet  [35] and the problem setup is depicted in Fig. 25. The time domain is chosen as T = [0, 5] s with time step size t = 2 × 10 −3 s. Before the beginning of the FSI simulation, the CFD problem is solved irrespective of the structure for 20 s with time step size equal to 5 × 10 −3 s in order to start the FSI simulation with a divergence free velocity field.
Concerning the partitioned FSI approach, the consistent mapping method is chosen for the transformation of the displacement fields whereas the conservative mapping method is chosen for the transformation of the consistent force vectors as described in "Fluid-structure interaction" section. Accordingly, nine simulations are performed: FSI simulations using a standard FEM structural discretization for all three inner pressure magnitudes (FSI-FEM-p1000, FSI-FEM-p2000 and FSI-FEM-p4000), FSI simulations for the standard FEM structural discretization with an ECL for all three inner pressure magnitudes (FSI-FEM-ECL-p1000, FSI-FEM-ECL-p2000 and FSI-FEM-ECL-p4000) and FSI simulations with the multipatch IGA structural discretization (FSI-IGA-p1000, FSI-IGA-p2000 and FSI-FEM-p4000) once more for all three inner pressure magnitudes. The geometric parametrization of the computational model corresponds to the fine setting investigated in [41] and it is used for both the ECL and the IGA discretization.
According to the methodological procedure concerning the isogeometric B-Rep mortarbased mapping method described in "Isogeometric B-Rep mortar-based mapping method on trimmed multipatches" section, the low order discretized surface is projected on the NURBS multipatch geometry. The projected structural finite element mesh onto the multipatch NURBS geometry in the frame of simulations FSI-IGA-p1000, FSI-IGA-p2000 and FSI-IGA-p4000 and the projected fluid FSI interface mesh onto the multipatch NURBS geometry concerning the simulations FSI-FEM-ECL-p1000, FSI-FEM-ECL-p2000 and FSI-EFM-ECL-p4000 are depicted in Figs. 27 and 28 , respectively, highlighting the boundary projection algorithm for elements with projection on more than one patch (see Fig. 10). Moreover, the fluid FSI interface mesh does not follow the exact torus shape and the tori comprising the hangar geometry from the fluid FSI interface mesh side are not connected with each other through a shared curved interface as for the finite element and the mul- tipatch IGA structural models but with straight planes allowing for a better fluid mesh between the tori. The latter choice is critical for the fluid mesh at these locations which otherwise would be highly distorted. Therefore, a gap between the tori can be observed in Fig. 28 as the elements comprising the planes between the tori have no unique projection on the multipatch NURBS surface. This however causes no problem to the displacement transformation as it can be seen in the sequel, but it does not allow for the consistent transformation of tractions. Therefore, the conservative mapping approach is herein chosen for the transformation of the consistent force vectors, see also in "Fluid-structure interaction" section. For the forthcoming investigations based on the time-displacement curves, the material point X m = 2.5 e 2 + 12.5 e 3 in the middle of the hangar is used.
The streamlines on and around the hangar are shown in the set of Fig. 29 for all employed simulations at t = 2 s demonstrating good qualitative accordance of the results when using standard finite element structural discretization, standard finite element structural discretization with the ECL and isogeometric structural discretization. The time displacement curves for all employed simulations at point X m along with the corresponding relative displacement error due to mapping are shown in Fig. 31. It can be seen that increasing the internal pressure of the tori (with according increase in the prestress) leads to structural displacements with lower amplitude as expected, see Fig. 31a. This shows that the structural response can be controlled through the internal pressure-prestress relationship without the necessity of adjusting the material properties for this kind of structures. In what concerns the displacement mapping error at X m , this stays very low throughout the The partitioned FSI scheme for this case is detailed in set of Fig. 30 whereas the depicted results are taken from the simulations with internal pressure of the tori equal to 1000 Pa. Accordingly, Figs. 30a and 30b demonstrate the concept of the partitioned FSI with consistent displacement mapping and conservative force mapping for the FSI simulations using standard finite element structural and isogeometric discretizations, respectively. The only difference between these two types of FSI simulations is the computation of the mortar based matrices, which for the former case are computed as in study [30] and for latter    case their computation is detailed in "Isogeometric B-Rep mortar-based mapping method on trimmed multipatches" section. Fig. 30c demonstrates the ECL concept where both the consistent displacement and the conservative force mapping are taking place through the ECL parametrized with the same CAD model used for the isogeometric structural discretization. Lastly, the results from the isogeometric B-Rep mortar-based mapping method for the transformation of the displacement field defined on the multipatch NURBS surface onto the fluid FSI interface mesh and from the structural finite element mesh to the multipatch NURBS surface are demonstrated. Accordingly, the FSI solution in terms of the displacement field from FSI-IGA-p1000 simulation at time t = 0.4 s defined on the multipatch NURBS geometry is transformed on the fluid FSI interface mesh, see Fig. 32, demonstrating an excellent accuracy even for highly non-matching surface representations. The corresponding relative error of the displacement in the L 2 (S )-norm is found in this case 8.462112E-03%. In the same way, the displacement field defined on the finite element discretized structural domain from the FSI-FEM-ECL-p1000 simulation also at time t = 0.4 s is transformed onto the multipatch NURBS surface, see Fig. 33. The relative displacement error in the L 2 (S )-norm is found in this case to be 3.277540E-03. Moreover, concerning the displacement and rotation interface jump in the L 2 (γ i )-norm, these are χ 0,γ i = 3.438509E − 04 m and χ 0,γ i = 1.811682E − 04 rad, respectively, whereas the L 2 ( d )-norm of the displacement field along the Dirichlet boundary d is computed

NREL phase VI wind turbine in numerical wind tunnel
In this section the FSI simulation of the NREL phase VI wind turbine with flexible blades introduced in "Isogeometric B-Rep analysis of the NREL phase VI wind turbine" section in numerical wind tunnel is investigated, see also in [23,53]. This example was also firstly presented in [35] and it is herein complemented with an additional study on the evolution of the interface work from each subdomain and their difference against the time only for the simulation involving the IBRA discretization of the blades. This is so because the coupling matrices are used in their original form for the computation of the interface work, see the preamble of "Numerical results" section, and the coupling matrices resulting from the standard mortar mapping method employed herein are modified for stability purposes as presented in the original work [63]. Therefore these matrices can not be used for this numerical example in order to evaluated the evolution of the interface work and therefore the corresponding results are not shown here. The selected material parameters for the flexible blades are presented in "Isogeometric B-Rep analysis of the NREL phase VI wind turbine" section and the kinematic viscosity for the air is given byν = 1.5451×10 −5 m 2 /s. The inlet velocity is chosen constant asυ = −7 e 2 m/s and the pressure is set equal to zero at the outlet, see Fig. 34. Accordingly, the side walls, the top and the bottom walls of the wind tunnel are set to slip boundary conditions. Moreover, the domain consists of two parts 4 : One non-rotating outer part and a rotating inner cylindrical part containing the rotor hub and the wind turbine blades which is rotating around X 2 -axis with constant angular velocity ω = 7.5398 rad/s. The fluid FSI interface is assigned to be the part containing the flexible blades. Note that S contains only the flexible blades across which the aerodynamic forces are acting, that is, the inner spars of the flexible blades are not part of S , see Fig. 4. The time interval for the coupled problem is chosen as T = [T 0 , T ∞ ] = [0, 5] s with a time step size of t = 10 −3 s. As in the previous numerical example, herein also the CFD problems is solved independently of the structure for 5 s with the same time step size as for the FSI simulation in order to have a divergence free fluid velocity field at the start of the FSI simulation. Concerning the CSD problem, the multipatch NURBS Kirchhoff-Love shell model with Penalty is employed, see "Isogeometric B-Rep analysis of the NREL phase VI wind turbine" section. Since the flexible blades are subject to constant angular velocity, the corresponding CSD problem is solved with time varying gravitational b g and constant centrifugal b c body forces given by, respectively, where X 1 is the distance from the center of rotation. The rotation tensor where ij 2 = αβ and δ ij 2 = δ αβ for α, β = 1, 3 stand for the components of the permutation and delta Kronecker tensor on the X 1 -X 3 plane, respectively, meaning that i2 2 = 2i 2 = δ i2 2 = δ 2i 2 = 0 for all i = 1, . . . , 3. Concerning the aerodynamic body forcesb acting along the flexible blades on S , these are computed similar to the gravitational body forces b g in Eq. (64a), namely, given that the fluid tractionst are referred to the current rotated configuration of the flexible blades at each time instance t. The latter approach allows for solving the CSD problem without considering inhomogeneous Dirichlet boundary conditions which accelerates the solution process. A limitation is however that only flexible blades rotating with constant angular velocity can be confronted with this approach, where in another case additional rotational inertial effects need to be addressed. Concerning the CFD setup, this is taken from the study in [53] and corresponding views of the CFD mesh with a close-up on the Fig. 35 NREL phase VI wind turbine in numerical wind tunnel: CFD computational domain [35] right blade are depicted in Fig. 35 where a mesh refinement in the neighborhood of the wind turbine can be observed and the CFD mesh consists of approximately ten million cells. Moreover, the sliding mesh interface using the Arbitrary Moving Interface (AMI) method provided in OpenFOAM ® is employed in order to couple the solution between the steady and rotating parts of the fluid domainṼ , see also in Fig. 34. For the CSD the standard finite element mesh of a shell with Reissner-Mindlin kinematics and the multipatch IGA model with Penalty and Kirchhoff-Love kinematics introduced in "Isogeometric B-Rep analysis of the NREL phase VI wind turbine" section are herein employed and evaluated. For the FSI simulation using the standard finite element mesh of the flexible blades, the standard mortar-based mapping method elaborated in [28] is used whereas for the FSI simulation using the multipatch IGA discretization of the flexible blades, the isogeometric B-Rep mortar-based mapping method introduced in "Isogeometric B-Rep mortar-based mapping method on trimmed multipatches" section is used. Accordingly, the mapped elements from the fluid FSI interface onto S of the multipatch NURBS surface are shown in Fig. 36 highlighting once more the excellent performance of the proposed methodology especially across the patch boundaries (see Fig. 10).
The Q-criterion 5 colored with the corresponding fluid velocity magnitude at exemplary time instances for both FSI simulations with the standard finite element and multipatch isogeometric discretizations is shown in the set of Figs. 37 where the deformation of the flexible blades is herein scaled by 170. The results demonstrate excellent qualitative accordance regardless of the highly diverse structural discretizations and mapping techniques, thus extending the isogeometric B-Rep mortar-based mapping method to real-world engineering applications.
Next, a quantitative comparison of the results in Fig. 38 is provided. Accordingly, the time displacement curves at the tip of the right blade X t = 5.029 e 1 − 0.013007 e 2 + 0.24821 e 3 [m] and the rotor shaft torque are depicted in Figs. 38a and 38b, respectively. The magnitude of the displacement field at the tip of the right blade from the CSD solution shows excellent accordance between the standard finite element mesh and the multipatch isogeometric discretization of the flexible wind turbine blades in terms of the pattern and frequency of the oscillations. However, the FEM solution exhibits slightly larger displacements which can be attributed to the underlying Reissner-Mindlin kinematics of the employed model for the standard FEM discretization in contrast to the Kirchhoff-Love   [35] shell kinematics associated with the multipatch isogeometric discretization of the flexible wind turbine blades. The error of the transformed displacement fields onto the fluid FSI interface at the tip is found negligible for this case. Concerning the rotor shaft torque, it can be observed that the pure CFD simulation produced the largest values and the FSI simulation with the FEM discretization of the flexible blades the lowest ones.
Subsequently     NREL phase VI wind turbine in numerical wind tunnel: Traction magnitude and the corresponding relative error of the transformation at X h on the upstream side versus time [35] respectively. It can be observed that taking into consideration the FSI coupling has an effect on the fluid traction field, see Fig. 40a, even for relatively small displacement fields as in this numerical example. Additionally, the relative error of the traction transformation at X h and X l versus time is shown in Figs. 39b and 40b , respectively, where it can be observed that both the standard and the isogeometric B-Rep mortar-based mapping methods produce excellent transformations for the traction fields on the flexible wind turbine blades.
Next, the transformation of fields using the isogeometric B-Rep mortar-based mapping method developed in "Isogeometric B-Rep mortar-based mapping method on trimmed multipatches" section is quantified. Firstly, the displacement field taken from the FSI simulation with the multipatch isogeometric discretization for the flexible blades at time t = 3 s are transformed onto the fluid FSI interface using Eq. (44), see Fig. 41, demonstrating excellent performance of the proposed method. For this case, the relative transformation error of the displacement field in the L 2 (S )-norm is found 0.038% which is highly satisfactory given the complexity and the size of the geometry. Moreover, the traction field    [35] from the same FSI simulation defined on the fluid FSI interface is taken at time t = 3 s and transformed onto the NURBS multipatch geometry S using Eq. (49), see Fig. 42. For the sake of clarity, both the high and the low pressure sides are herein depicted, see Figs. 42a and 42b, respectively. The results show once more excellent accordance even for a highly oscillatory field, such as the traction field in this case, and the corresponding transformation error in the traction field is in this case found 4.31% based on the L 2 (S )-norm whereas the interface jump of the traction field between the multipatches is equal to χ 0,γ i = 0.215 m. Lastly, the evolution of the interface work from both the fluid and the structural sides versus the simulation time and the corresponding residual energy are shown in set of Fig. 43 for the simulation with the IBRA structural discretization. The results show a satisfactory energy transfer across the interface and the corresponding residual energy at the interface stays at low levels for this large scale application. The corresponding results for the standard structure FEM discretization where the standard mortar-based mapping method is used are not shown herein. The reason is that the corresponding implementation of the standard mortar-based mapping method used herein is the one developed in [28] where many robustness enhancements are taken into consideration in order to render the methodology applicable for real-world engi-  neering problems. Accordingly, whenever a projection of a finite element fails the mortar transformation matrices for the standard mortar-based mapping method are enhanced with components stemming from a nearest neighbor method and additionally consistency enforcement is applied when the two interface meshes do not exactly overlap which is standard for real-world applications. Therefore, the matrices resulting from this standard mortar-based mapping method can not be used for the evaluation of the interface work and the corresponding results are omitted. However, the purpose of this study is to highlight the advantages of the newly proposed IBRA mortar-based mapping method and its consistency in terms of the satisfaction of the interface work.
It is important to note that the computational footprint of the mortar-based mapping method is negligible in the context of partitioned FSI, regardless of whether the standard finite element or the isogeometric B-Rep mortar-based mapping method is used. That is so, because the coupling matrices in terms of the mortar-based mapping method are computed only once at the beginning of the coupled simulation and then used throughout the simulation by means of simple matrix-vector multiplications, see for instance Eqs. (44), (49), (59) and (60). The overall additional computational overhead of the mortarbased mapping method is then negligible in this context, as the coupled simulation may take days or even months to complete, given the computational expense that 3D CFD simulations involve. It is demonstrated in study [52] that it takes about seven times longer to generate the coupling matrices for the isogeometric as compared to the standard finite element mortar-based mapping method based on the FSI simulation of a flexible hemisphere discretized with a single untrimmed patch. The computational overhead when computing the coupling matrices in the frame of the ECL is simply the sum of the computational overheads involved by the standard finite element and the isogeometric mortar-based mapping, see Fig. 30c for an illustration or Eq. (31) in study [52]. In what concerns the mapping over the NREL Wind Turbine blades presented in this section, it was found that it takes about 3 and 20 minutes for the generation of the coupling matrices regarding the standard finite element and the isogeometric mortar-based mapping methods, respectively, which is in accordance to the findings in study [52]. This is to be expected, since the isogeometric mortar-based mapping method involves additional algorithmic steps for the generation of the coupling matrices, see "Realization" section. However, these findings cannot be used conclusively concerning the comparison of the computational efficiency between the two mortar-based mapping methods. It is therefore encouraged to develop an efficient implementation of the isogeometric mortar-based mapping method in a future study, in order to be able to provide conclusive results regarding the comparison of the computational efficiency of both methods, given that such a comparison is out of the scope of the present study. Overall it can be said that FSI with IGA is highly efficient by means of the proposed isogeometric B-rep mortar-based mapping method for real-world engineering problems.

Conclusions
In this contribution, a mortar-based mapping method is elaborated and assessed for its application to field transformations between trimmed NURBS-based CAD models and standard low order discretizations (FEM, FVM etc.) of a surface. The application of the aforementioned method considered herein is that of the partitioned FSI simulations, either directly involving isogeometric structural discretizations or using the geometric parametrization of an Exact Coupling Layer (ECL) for smoothing the description of the interface fields. However, the herein proposed isogeometric B-Rep mortar-based mapping method can be also applied to the regeneration of CAD B-Rep models, see the work done in [72]. Three numerical examples are used in order to validate and assess the proposed methodology. The lid-driven cavity is employed as benchmark example, whereas the FSI simulations of an inflatable hangar and the NREL phase VI wind turbine with flexible blades are used in the context of real-word applications. The results clearly show, that real-world CAD models involving trimmed multipatches can be efficiently used in the context of partitioned FSI by means of the proposed isogeometric B-Rep mortar-based mapping method. Additional focus is given to the use of the isogeometric B-Rep mortarbased mapping method in conjunction with the CAD description of the interface as an ECL in order to obtain smooth results that are highly desirable in the context of surface coupled problems and in particular in the context of FSI. The results section is complemented with convergence studies and energy transfer evaluations of the proposed method. The convergence studies are concerned with the assessment of the errors in the L 2 -norm regarding the mapping of given fields from the trimmed isogeometric multipatch surface to the finite element mesh and the other way around, demonstrating consistence of the proposed method. The energy transfer evaluations are concerned with the transferred work through the interface for the corresponding FSI simulations, demonstrating inherent physical relevance of the method. It can be seen that although the proposed method is not by construction energy-conservative, the energy gain or loss when transferring fields by means of the proposed isogeometric B-Rep mortar-based mapping method is minimal. Therefore, it can be concluded that the herein proposed isogeometric B-Rep mortarbased mapping method extends IBRA to FSI in a consistent and efficient computational framework.