# Systematic evaluation of the interface description for fluid-structure interaction simulations using the isogeometric mortar-based mapping

Andreas Apostolatos, Guillaume de Nayer, Kai-Uwe Bletzinger, Michael

Breuer, Roland Wüchner

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Andreas Apostolatos, Guillaume de Nayer, Kai-Uwe Bletzinger, Michael Breuer, Roland Wüchner. Systematic evaluation of the interface description for fluid-structure interaction simulations using the isogeometric mortar-based mapping. International Journal of Heat and Fluid Flow, 2019, 86, pp.368-399. 10.1016/j.jfluidstructs.2019.02.012. hal-02132323

HAL Id: hal-02132323

https://hal.science/hal-02132323

Submitted on 17 May 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Systematic evaluation of the interface description for fluid-structure interaction simulations using the isogeometric mortar-based mapping

<figure><img src="image_3.png" /></figure>

Check for
updates

Andreas Apostolatos<sup>a,*</sup>, Guillaume De Nayer<sup>b</sup>, Kai-Uwe Bletzinger<sup>a</sup>, Michael Breuer<sup>b</sup>, Roland Wüchner<sup>a</sup>

## ARTICLE INFO

**Article history:**
Received 11 October 2018
Received in revised form 5 February 2019
Accepted 12 February 2019
Available online 13 March 2019

**Keywords:**
Fluid-structure interaction (FSI)
Isogeometric analysis (IGA)
Isogeometric mortar-based mapping
method
Membrane structural analysis
Large-eddy simulation (LES)

## ABSTRACT

Within this study the influence of the interface description for partitioned *Fluid-Structure Interaction* (FSI) simulations is systematically evaluated. In particular, a *Non-Uniform Rational B-Spline* (NURBS)-based isogeometric mortar method is elaborated which enables the transfer of fields defined on low-order and isogeometric representations of the interface along which the FSI constraints are defined. Moreover, the concept of the *Exact Coupling Layer* (ECL) using the proposed isogeometric mortar-based mapping method is presented. It allows for smoothing fields that are transferred between two standard low-order surface discretizations applying the exact interface description in terms of NURBS. This is especially important for highly turbulent flows, where the artificial roughness of the low-order faceted FSI interfaces results in spurious flow fields leading to inaccurate FSI solutions. The approach proposed is subsequently compared to the standard mortar-based mapping method for transferring fields between two low-order surface representations (finite volume method for the fluid and finite element method for the structure) and validated on a simple lid-driven cavity FSI benchmark. Then, the physically motivated 3D example of the turbulent flow around a membranous hemisphere (Wood et al., 2016) is considered. Its behavior is predicted by combining the large-eddy simulation technique with the isogeometric analysis to demonstrate the usefulness of the isogeometric mortar-based mapping method for real-world FSI applications. Additionally, the test case of a bluff body significantly deformed in an eigenmode shape of the aforementioned hemisphere is used. For this purpose, both “standard” low-order finite element discretizations and a smooth IGA-based description of the structural surface are considered. This deformation is transferred to the fluid FSI interface and the influence of the interface description on the fluid flow is analyzed. Finally, the computational costs related to the presented methodology are evaluated. The results suggest that the proposed methodology can effectively improve the overall FSI behavior with minimal effort by considering the exact geometry description based on the *Computer-Aided Design* (CAD) model of the FSI interface.

© 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

<footnote><sup>a</sup> Chair of Structural Analysis, Technical University of Munich, D-80333 Munich, Germany</footnote>

<footnote><sup>b</sup> Department of Fluid Mechanics, Helmut-Schmidt-University, D-22043 Hamburg, Germany</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: andreas.apostolatos@tum.de (A. Apostolatos).</footnote>
Nomenclature

<table>
  <tr>
    <td>Indices α, β, ...</td>
    <td>Range from 1 to 2</td>
  </tr>
  <tr>
    <td>Indices i, j, ...</td>
    <td>Range from 1 to 3</td>
  </tr>
  <tr>
    <td>X<sub>i</sub>, x<sub>i</sub></td>
    <td>Cartesian coordinates in the reference and current configuration</td>
  </tr>
  <tr>
    <td>Ω, Ω<sub>t</sub></td>
    <td>Surface domain in the reference and current configuration</td>
  </tr>
  <tr>
    <td>Ω̂</td>
    <td>Parametric space of the surface</td>
  </tr>
  <tr>
    <td>Ω<sub>h</sub>, Ω<sub>FSI</sub></td>
    <td>Discrete representation of the surface and FSI interface</td>
  </tr>
  <tr>
    <td>Γ<sub>d</sub></td>
    <td>Dirichlet boundary of the structure</td>
  </tr>
  <tr>
    <td>S, s</td>
    <td>Parametrization of the surface in the reference and current configuration</td>
  </tr>
  <tr>
    <td>θ<sub>i</sub></td>
    <td>Parametric coordinates</td>
  </tr>
  <tr>
    <td>A<sub>α</sub>, a<sub>α</sub></td>
    <td>Covariant base vectors in the reference and current configuration</td>
  </tr>
  <tr>
    <td>A<sup>α</sup></td>
    <td>Contravariant base vectors in the reference configuration</td>
  </tr>
  <tr>
    <td>ε</td>
    <td>Green–Lagrange strain tensor of the structure</td>
  </tr>
  <tr>
    <td>ε<sub>αβ</sub></td>
    <td>Covariant components of the Green–Lagrange strain tensor of the structure</td>
  </tr>
  <tr>
    <td>S<sup>αβ</sup></td>
    <td>Contravariant components of the 2nd Piola–Kirchhoff stress tensor of the structure</td>
  </tr>
  <tr>
    <td>n</td>
    <td>2nd Piola–Kirchhoff stress resultant force tensor of the structure</td>
  </tr>
  <tr>
    <td>n<sub>0</sub></td>
    <td>Pre-stress force tensor of the structure</td>
  </tr>
  <tr>
    <td>C</td>
    <td>Material tensor of the structure</td>
  </tr>
  <tr>
    <td>d, ḋ, ḋ̃</td>
    <td>Displacement, velocity and acceleration fields of the structure</td>
  </tr>
  <tr>
    <td>t<sub>s</sub>, ρ<sub>s</sub>, E, ν</td>
    <td>Thickness, density, Young's modulus and Poisson's ratio of the structure</td>
  </tr>
  <tr>
    <td>N<sub>i</sub></td>
    <td>Standard scalar-valued finite element basis functions</td>
  </tr>
  <tr>
    <td>N<sup>pα</sup><sub>i</sub>, R<sup>pα</sup><sub>i</sub></td>
    <td>Scalar-valued B-Spline and NURBS basis functions in θ<sub>α</sub> parametric direction</td>
  </tr>
  <tr>
    <td>φ<sup>h</sup><sub>r</sub>, φ<sup>h</sup><sub>r</sub></td>
    <td>Vector-valued standard finite element and NURBS basis functions</td>
  </tr>
  <tr>
    <td>Θ<sub>α</sub></td>
    <td>Knot vector in θ<sub>α</sub> parametric direction</td>
  </tr>
  <tr>
    <td>m<sub>α</sub>, n<sub>α</sub>, p<sub>α</sub></td>
    <td>Number of knots, number of control points and polynomial order in θ<sub>α</sub> parametric direction</td>
  </tr>
  <tr>
    <td>w<sub>i,j</sub></td>
    <td>(i, j)-th control point weight of a surface</td>
  </tr>
  <tr>
    <td>ê̂, ḋ̃, ḋ̃̃</td>
    <td>Vector of discrete control point displacement, velocity and acceleration DOFs</td>
  </tr>
  <tr>
    <td>M, D, K, K̄</td>
    <td>Mass, damping, steady-state tangent stiffness matrix and dynamic stiffness matrix of the structure</td>
  </tr>
  <tr>
    <td>R, R̄</td>
    <td>Steady-state and dynamic residual vector of the structure</td>
  </tr>
  <tr>
    <td>F̂<sub>t</sub></td>
    <td>External discrete force vector at time t</td>
  </tr>
  <tr>
    <td>Δt</td>
    <td>Time-step size</td>
  </tr>
  <tr>
    <td>β<sub>n</sub>, γ<sub>n</sub></td>
    <td>Newmark parameters</td>
  </tr>
  <tr>
    <td>α<sub>t</sub>, β<sub>t</sub></td>
    <td>Rayleigh damping parameters</td>
  </tr>
  <tr>
    <td>n<sub>h</sub></td>
    <td>Number of nodes in the finite element discretization</td>
  </tr>
  <tr>
    <td>C<sup>nn</sup>, C<sup>rr</sup></td>
    <td>Mortar mapping matrices when transferring a field from IGA to FEM and from FEM to IGA</td>
  </tr>
  <tr>
    <td>g</td>
    <td>Gravitational acceleration</td>
  </tr>
  <tr>
    <td>u, v, w</td>
    <td>Cartesian components of the flow velocity field</td>
  </tr>
  <tr>
    <td>U<sub>∞</sub></td>
    <td>Free-stream velocity</td>
  </tr>
  <tr>
    <td>μ</td>
    <td>Dynamic viscosity of the fluid</td>
  </tr>
  <tr>
    <td>D</td>
    <td>Diameter of the flexible hemisphere</td>
  </tr>
  <tr>
    <td>Re</td>
    <td>Reynolds number of the flow</td>
  </tr>
  <tr>
    <td>δ</td>
    <td>Thickness of the boundary layer</td>
  </tr>
  <tr>
    <td>Δp<sub>FSI</sub></td>
    <td>Pressure difference between the inner and the outer domain in the FSI hemisphere case</td>
  </tr>
  <tr>
    <td>R<sub>a</sub>, R<sub>q</sub></td>
    <td>Arithmetic average and root-mean-square roughness</td>
  </tr>
  <tr>
    <td>A<sub>α</sub> · A<sub>β</sub></td>
    <td>Scalar product of two vectors A<sub>α</sub> and A<sub>β</sub></td>
  </tr>
  <tr>
    <td>A<sub>α</sub> × A<sub>β</sub></td>
    <td>Cross product of two vectors A<sub>α</sub> and A<sub>β</sub></td>
  </tr>
  <tr>
    <td>A<sub>α</sub> ⊗ A<sub>β</sub></td>
    <td>Dyadic product of vectors A<sub>α</sub> and A<sub>β</sub></td>
  </tr>
  <tr>
    <td>C : ε</td>
    <td>Double contraction of tensors C and ε</td>
  </tr>
  <tr>
    <td>||·||<sub>2</sub></td>
    <td>2-norm in Euclidean space</td>
  </tr>
  <tr>
    <td>•*</td>
    <td>Dimensionless quantity</td>
  </tr>
</table>
# 1. Introduction

*Fluid-Structure Interaction* (FSI) simulations play an important role in modern engineering for the accurate prediction of phenomena which govern the mutual interaction between a fluid flow and a flexible structure. Their importance is especially highlighted in that this kind of events are governed by a highly nonlinear coupled behavior which cannot be predicted using solely a structural or a fluid simulation code but rather both together in a coupled environment. Within this study, the iterative Gauss-Seidel partitioned approach (Sicklinger, 2014; Sicklinger et al., 2014) is used for the FSI simulation, which is also known as serial staggered procedure with an iterative correction loop (Felippa and Park, 1980). It allows for a highly modular framework within which various fluid or structural solvers can be employed, whereas the fulfillment of the interface constraints is obtained in an iterative manner.

At the so-called *wet surface* between the fluid and the structural domains, the interface constraints have to be fulfilled typically at distinct descriptions of the FSI interface from each subdomain. This is typically due to the fact that diverse solution approaches are employed for solving the underlying fluid and structural *Boundary Value Problems* (BVPs) which are characterized by non-matching interface discretizations. For instance, the structural BVP is typically solved numerically using the *Finite Element Method* (FEM), whereas the fluid BVP is often numerically solved using the *Finite Volume Method* (FVM). Additionally, highly diverse mesh resolutions are encountered: Typically, fluid simulations need highly resolved boundary layers around the FSI interface, especially when applying the large-eddy simulation technique for turbulent flow fields, whereas a coarser discretization is typically sufficient for obtaining a satisfactory structural solution. Therefore, schemes have to be developed which allow the transformation of interface fields, for instance structural displacements and fluid tractions, across the common interface with minimal error. These methods comprise the *nearest neighbor* method, the *nearest element* interpolation, the *barycentric* interpolation and the *mortar* method amongst others (see Park et al. (2002), de Boer et al. (2007, 2008) and Wang (2016) for more details). It is especially important that the structural displacements are smoothly transferred onto the fluid FSI interface. Kinks on the surface due to a coarse finite element resolution on the structural side may lead to artificial shifts in separation and reattachment points which substantially influence the flow development and which in turn can cause the evolution of inaccurate deformation patterns due to the coupling.

Herein, standard finite element and isogeometric discretizations of geometrically nonlinear pre-stressed membrane structures are considered concerning the structural BVP. The widely used FVM method is employed for the discretization of the incompressible Navier-Stokes equations. Isogeometric analysis (IGA) is a recent numerical method for the approximation of BVPs which uses the exact geometric parametrization, in this study the NURBS basis functions, for the approximation of the problem's unknown field(s). Among others, smoothness and high-order field approximations are some of the main benefits of IGA (see Cottrell et al. (2009) for more information). For the transfer of the structural displacement and fluid traction fields across the FSI interface, the mortar method is chosen in the present study. This method was initially developed for enforcing the continuity of the solution to variational problems on non-overlapping subdomains, which attain non-conforming discretizations along their common interfaces, see Bernardi et al. (1994). It was then used in partitioned FSI simulations with standard discretizations of the underlying BVPs in de Boer et al. (2008) and Klöppel et al. (2011) for coupling the underlying physics between the fluid and structural subdomains. Moreover, a special mortar formulation was presented in De Lorenzis et al. (2012) in the context of contact mechanics. The extension of the mortar method for FSI simulations involving isogeometric structural discretizations and the *Variational Multiscale* (VMS) stabilized FEM (Codina, 2001; Bazilevs et al., 2007) for the fluid discretization was presented in Bazilevs et al. (2012) and denoted as $L^2$-projection. The methodology of the present work bases its fundamentals on this aforementioned $L^2$-projection. Moreover, this study details a framework for computing the mortar transformation matrices in an optimal way. Particularly, the presented isogeometric mortar-based mapping method deals with highly diverse, non-matching NURBS-based parametrizations stemming from Computer-Aided Design (CAD) and low-order discretizations of the common interface. Additionally, the latter approach is also used in the context of the *Exact Coupling Layer* (ECL) method, which can be used for smoothing fields transformed between two low-order discretizations of the common interface through the NURBS-parametrized model of the exact interface description in terms of CAD. As detailed in Wang et al. (2016) the mortar method is superior to the other aforementioned mapping methods in terms of accuracy and consistency for the transformation of fields across non-matching discretizations with highly diverse mesh densities. Additionally, since the geometries of the computational models are generated in CAD software, it is highly efficient to use IGA, avoiding the time-consuming mesh generation step which is required in standard finite element methods. It is concluded that the most suitable method for transferring fields between a low-order discretized and a NURBS parametrized surface is the mortar-based method. Besides the aforementioned advantages, another main reason is that the *Degrees of Freedom* (DOFs) of the isogeometrically discretized field(s) lie on the *Control Points* (CPs). Thus, they do not represent physically relevant values since the CPs do not necessarily interpolate the geometry (Cottrell et al., 2009). This renders nearest neighbor or nearest element methods inappropriate. The isogeometric mortar-based mapping method naturally allows a two-fold strategy: It can be used in a coupled FSI scenario when the structural BVP is discretized using IGA (Breitenberger et al., 2015; Philipp et al., 2016; Wüchner et al., 2018) or it can be used to smooth fields when transformed from one low-order discretization (e.g., FEM) to another (e.g., FVM) taking advantage of the high-order nature and smoothness that the NURBS basis functions offer. Therefore, in the latter case the interface can be exactly represented in terms of NURBS, while the transformation of the interface fields can be done through this surface, denoted as the *Exact Coupling Layer* (ECL).

The paper is organized as follows: Firstly, the numerical methodology is detailed. The membrane BVP is briefly introduced along with a short discussion on the NURBS basis functions and the NURBS-based IGA spatial discretization
<figure><img src="image_4.png" /><figcaption>Fig. 1. Continuum mechanics of thin-walled structures: Mid-surface reduction and membrane stress state (Apostolatos, 2019).</figcaption></figure>

of the underlying variational form. This section is complemented with the time discretization of the structural semi-discretized weak form. Then, a short description of the solution of the fluid problem is given. The Navier–Stokes equations for an incompressible fluid are written in an Arbitrary Lagrangian Eulerian (ALE) frame of reference and discretized using finite volumes. In case of turbulent flows, the LES technique is applied. Subsequently, the isogeometric mortar-based mapping method is introduced along with the employed assumptions and implementation aspects. The theoretical part is concluded by introducing the FSI solution approach including the residual energy computation associated with the coupling interface. For the assessment of the proposed methodology, three examples are chosen. The first one is the lid-driven cavity FSI benchmark proposed by Wall (1999). It is used for validation purposes, since the flow is laminar and the geometry simple. Then, an air-inflated hemisphere exposed to a turbulent boundary layer is investigated in order to demonstrate the applicability of the enhanced methodology for problems of physical relevance. Experimental and numerical reference data allow a direct comparison. Furthermore, in order to evaluate the influence of the FSI interface representation on the fluid solution, a significantly larger and time-invariant deformation is imposed to the hemispherical membrane. Finite element discretizations of different resolution for the deformed shape are subsequently considered as well as a NURBS parametrization. The standard mortar-based mapping method is used for the finite element discretizations, whereas the isogeometric mortar-based mapping method is employed for the NURBS-parametrized deformed shape. The results highlight the importance of the accurate interface representation in FSI simulations in terms of accuracy. Lastly, the computational costs related to the proposed methodology, particularly for the description of the FSI interface, are evaluated.

## 2. Numerical methodology

The numerical simulation methodology is based on a strong and efficient coupling scheme especially developed for fluid-structure interaction (FSI) problems including thin-walled structures within laminar or turbulent flows. The method relies on a partitioned coupling procedure. Thus, the structural and the fluid problems are solved separately. Since the details concerning the fluid part and the FSI coupling are included in Breuer et al. (2012) and corresponding validation studies can be found in De Nayer et al. (2014) and De Nayer and Breuer (2014), these parts are kept short. Herein, the focus is put on the isogeometric description of the membrane structure and on the isogeometric FSI mapping strategy.

### 2.1. Computational structural dynamics of membranes

In the following the basic equations governing the computational structural dynamics of membranes are detailed. Accordingly, this section starts with the continuum mechanics of membranes, continues with the isogeometric spatial semi-discretization and complemented with the time discretization using the Newmark scheme.

#### 2.1.1. Continuum mechanics of membranes

Within this study, the interest concerning Computational Structural Dynamics (CSD) is focused on so-called thin-walled structures whose mechanics are governed by the membrane theory (Apostolatos et al., 2019). Subsequently, let $S: \hat{\Omega} \to \Omega$ be a surface in 3D, where $\hat{\Omega} \subset \mathbb{R}^2$ and $\Omega \subset \mathbb{R}^3$ are its parametric and physical images, respectively. The domain $\hat{\Omega}$ is also known as the parametric image of the physical domain $\Omega$. The parametric coordinates $\theta_1$ and $\theta_2$ are defined along the surface, whereas $\theta_3$ is chosen in the direction normal to the surface such that $\theta_3 \in [-t_s/2, t_s/2]$. Here, $t_s$ represents the thickness of the actual structure and is considered constant. Therefore, the domain $\Omega$ consists of all material points $\mathbf{X}$ which have a parametric image in $\hat{\Omega}$ through $S$, that is, $\mathbf{X} = S(\theta_1, \theta_2)$. The CSD problem is then characterized by the unknown displacement field $\mathbf{d}$ in a total Lagrangian description of the motion, meaning that the underlying stresses and strains are defined with respect to the initial undeformed configuration. In this way, the curvilinear space is bent
on the physical space following the deformation of the structure. The curvilinear base vectors on the surface $\Omega$ may be constructed (Basar and Krätzig, 1985) as follows,

$$ \mathbf{A}_{\alpha} = \frac{\partial \mathbf{S}}{\partial \theta_{\alpha}}, \qquad (1a) $$

$$ \mathbf{A}_3 = \frac{1}{\|\mathbf{A}_1 \times \mathbf{A}_2\|_2} \mathbf{A}_1 \times \mathbf{A}_2 . \qquad (1b) $$

The displacement field can be defined in either the Cartesian or the aforementioned curvilinear space:

$$ \mathbf{d} = d_i^0 \mathbf{e}_i = d^\alpha \mathbf{A}_\alpha + d_3 \mathbf{A}_3 . \qquad (2) $$

In Eq. (2), $d_i^0$, $\mathbf{e}_i$ and $(d^\alpha, d_3)$ are the Cartesian components of the displacement field, the Cartesian base vectors and the contravariant components of the displacement field, respectively. The base vectors of the current configuration $\mathbf{a}_\alpha$ and $\mathbf{a}_3$ are defined in a similar manner as Eq. (1) using the surface in the current configuration, that is, $\mathbf{s} : \hat{\Omega} \to \Omega_t$. Herein $\Omega_t$ represents the current configuration consisting of all material points $\mathbf{x} = \mathbf{X} + \mathbf{d}$ for which a parametric image through $\mathbf{s}$ exists, namely $\mathbf{x} = \mathbf{s}(\theta_1, \theta_2)$ at a given time instance $t$. Concerning the kinematics, given the contravariant base vectors $\mathbf{A}^\alpha$ on surface $\Omega$, the nonlinear Green-Lagrange second-order strain tensor $\boldsymbol{\varepsilon} = \varepsilon_{\alpha\beta} \mathbf{A}^\alpha \otimes \mathbf{A}^\beta$ is employed which is defined over the contravariant basis (Holzapfel, 2000), and whose covariant components read,

$$ \varepsilon_{\alpha\beta} = \frac{1}{2} (\mathbf{a}_{\alpha\beta} - \mathbf{A}_{\alpha\beta}), \qquad (3) $$

where $\mathbf{A}_{\alpha\beta}$ and $\mathbf{a}_{\alpha\beta}$ stand for the components of the metric tensors of the reference and the current configuration, respectively,

$$ \mathbf{A}_{\alpha\beta} = \mathbf{A}_{\alpha} \cdot \mathbf{A}_{\beta}, \qquad (4a) $$

$$ \mathbf{a}_{\alpha\beta} = \mathbf{a}_{\alpha} \cdot \mathbf{a}_{\beta}. \qquad (4b) $$

It is important to note that the strains in Eq. (3) are assumed to be constant along the thickness thus allowing only in-plane strains, as typical for membrane structures. Within the thin structure assumption $r/t_s \ll 20$, $r$ being the radius of curvature, the out-of-plane strain and stress components in the continuum, i.e., in $\theta_3$-direction, can be neglected. The stress-strain relations are established using the Saint-Venant-Kirchhoff material law. For a given fourth-order material tensor $\mathbf{C} = C^{\alpha\beta\gamma\delta} \mathbf{A}_{\alpha\beta} \otimes \mathbf{A}_{\beta\gamma} \otimes \mathbf{A}_{\delta\alpha}$ this relation reads:

$$ \mathbf{n} = \mathbf{C} : \boldsymbol{\varepsilon}. \qquad (5) $$

Herein, $\mathbf{n} = n^{\alpha\beta} \mathbf{A}_{\alpha} \otimes \mathbf{A}_{\beta}$ stands for the 2nd Piola-Kirchhoff second-order stress resultant force tensor. The components of the stress resultant force tensor $n^{\alpha\beta}$ are given as pre-integrated components of the actual stress tensor components $S^{\alpha\beta}$ along the thickness (see Fig. 1), that is, $n^{\alpha\beta} = t_s S^{\alpha\beta}$, since the stresses are assumed to be constant across the thickness (Bletzinger et al., 2005). The name stress resultant force tensor is herein used, since after the pre-integration $\mathbf{n}$ represents neither a stress nor a force for this kind of structures (Kiendl, 2011).

Membrane structures in 3D typically require a tensile pre-stress field, herein represented by the second-order tensor $\mathbf{n}_0$, in order to withstand general external forces which have also a transversal component (Bletzinger and Ramm, 1999). In this way, the weak form of the dynamic equilibrium for membrane structures can be formulated as follows: Find the displacement field $\mathbf{d}$ such that,

$$ \int_{\Omega} \delta \mathbf{d} \cdot \rho_s t_s \ddot{\mathbf{d}} \, d\Omega + \int_{\Omega} \delta \mathbf{d} \cdot c_s t_s \dot{\mathbf{d}} \, d\Omega + \int_{\Omega} \delta \boldsymbol{\varepsilon} : (\mathbf{n} + \mathbf{n}_0) \, d\Omega = \int_{\Omega} \delta \mathbf{d} \cdot \mathbf{b} \, d\Omega, \qquad (6) $$

for all test displacement fields $\delta \mathbf{d}$ which fulfill the Dirichlet boundary conditions along a portion of the boundary $\Gamma_d \subset \partial \Omega$. In Eq. (6), $\dot{\mathbf{d}}$ is the velocity, $\ddot{\mathbf{d}}$ the acceleration, $\rho_s$ the density and $c_s$ the damping coefficient of the structure. Pre-integrated body forces and surface tractions are contained in vector $\mathbf{b}$. However, in the present work the damping term is approximated using the Rayleigh damping approach in the discrete form of the problem as described in Section 2.1.3.

### 2.1.2. Spatial discretization using the isogeometric structural analysis

Isogeometric Analysis (IGA) is a recent numerical method for solving Partial Differential Equations (PDEs) using the exact geometric information for the computational model (Hughes et al., 2005; Cottrell et al., 2009) as opposed to “standard” Finite Element Methods (FEMs), where the geometry of the computational model is typically approximated by lower-order polynomials. Herein, the Non-Uniform Rational B-Spline (NURBS)-based IGA is employed for the spatial discretization of the PDE describing the structural mechanics (Apostolatos et al., 2014, 2015). Thus, the geometric parametrization of the reference and the current configuration, $\mathbf{S}$ and $\mathbf{s}$, respectively, is done using the NURBS basis functions.

In the following, Pieglin and Tiller (2012) is closely followed for defining the NURBS and the B-Spline basis functions as well as the corresponding parametrically defined surfaces. Given two knot vectors $\Theta_\alpha = [\theta_1^\alpha \dots \theta_{m_\alpha}^\alpha]$, $\theta_i^\alpha$ and $m_\alpha$ being
the set of knots and the total number of knots within each knot vector, the two-dimensional NURBS basis functions are defined by means of the two-dimensional B-Spline basis functions $N_{i,j}^{p_1 \cdot p_2}$ and a set of weights $w_{i,j}$ as,

$$R_{i,j}^{p_1 \cdot p_2} (\theta_1, \theta_2) = \frac{w_{i,j} N_{i,j}^{p_1 \cdot p_2} (\theta_1, \theta_2)}{\sum_{k=1}^{n_1} \sum_{l=1}^{n_2} w_{k,l} N_{k,l}^{p_1 \cdot p_2} (\theta_1, \theta_2)}, \quad \forall (\theta_1, \theta_2) \in \hat{\Omega}, \quad (7)$$

where $\hat{\Omega} = [\theta_1^1, \theta_1^1] \times [\theta_2^1, \theta_2^1]$. The two-dimensional B-Spline basis functions $N_{i,j}^{p_1 \cdot p_2}$ are then defined by the one-dimensional B-Spline basis functions $N_i^{p_\alpha}$ in $\Theta_\alpha$ as follows,

$$N_{i,j}^{p_1 \cdot p_2} (\theta_1, \theta_2) = N_i^{p_1} (\theta_1) N_j^{p_2} (\theta_2), \quad \forall (\theta_1, \theta_2) \in \hat{\Omega}, \quad (8)$$

where $p_\alpha$ stands for the polynomial order of B-Spline basis functions $N_i^{p_\alpha}$ in $\theta_\alpha$-direction. The one-dimensional B-Spline basis functions $N_i^{p_\alpha}$ are defined by means of the Cox-De Boor recursion formula (Piegl and Tiller, 2012), that is,

$$N_i^{\bar{p}_\alpha}(\theta_\alpha) = \frac{\theta_\alpha - \theta_i^\alpha}{\theta_{i+\bar{p}_\alpha}^\alpha - \theta_i^\alpha} N_i^{\bar{p}_\alpha-1}(\theta_\alpha) + \frac{\theta_{i+\bar{p}_\alpha+1}^\alpha - \theta_\alpha}{\theta_{i+\bar{p}_\alpha+1}^\alpha - \theta_{i+1}^\alpha} N_{i+1}^{\bar{p}_\alpha-1}(\theta_\alpha), \quad \forall \theta_\alpha \in \Theta_\alpha \quad (9)$$

where $\bar{p}_\alpha = 0, \dots, p_\alpha$. Thus, B-Spline basis functions $N_i^{\bar{p}_\alpha}$ are obtained by a recursive construction $\bar{p}_\alpha = 0, \dots, p_\alpha$. Additionally, for the piecewise constant basis functions it holds $N_i^0(\theta_\alpha) = 1$ for $\theta_\alpha \in [\theta_i^\alpha, \theta_{i+1}^\alpha)$, while $N_i^{p_\alpha}(\theta_\alpha) = 0$ elsewhere and the definition $0/0 = 0$ is assumed in Eq. (9). The B-Spline and subsequently the NURBS basis functions attain $C^\infty$-continuity within each knot span $[\theta_\alpha^i, \theta_\alpha^{i+1}]$ and $C^{p_\alpha-k_\alpha}$-continuity across knots $\theta_\alpha^i$, whereas $k_\alpha^i$ represents the multiplicity of knots $\theta_\alpha^i$ in $\Theta_\alpha$. Given a set of Control Points (CPs) with the associated position vector $\hat{\mathbf{X}}_{i,j}$ in the physical space, the NURBS surface $\mathbf{S}: \hat{\Omega} \to \Omega$ can be constructed as

$$\mathbf{S} = \sum_{(i,j) \in \mathcal{I}} R_{i,j}^{p_1 \cdot p_2} \hat{\mathbf{X}}_{i,j}, \quad (10)$$

where $\mathcal{I}$ contains all pairs $(i, j)$ with $i = 1, \dots, n_1, j = 1, \dots, n_2$ and $n_\alpha = m_\alpha - p_\alpha - 1$ being the number of basis functions along $\Theta_\alpha$. Within NURBS-based IGA, the primal variable, herein the displacement field $\mathbf{d}$, is approximated using the aforementioned NURBS basis functions. For simplicity, let the aforementioned pairs $(i, j) \in \mathcal{I}$ be sequentially sorted through a map $(i, j) \to k$ for all $i = 1, \dots, n_1$ and $j = 1, \dots, n_2$, where $k = 1, \dots, n$ and $n = n_1 n_2$. In this way, the approximation of the displacement field using NURBS reads,

$$\mathbf{d} \approx \sum_{r=1}^{3n} \phi_r \hat{d}_r, \quad (11)$$

where $\phi_r$ stands for the vector-valued NURBS basis functions, namely, $\phi_r = R_k^{p_1 \cdot p_2} \mathbf{e}_i$ with $(k, i) = (\lceil r/3 \rceil, r - 3\lceil r/3 \rceil + 3)$ for all $r = 1, \dots, 3n$. In this way, each pairing of a two-dimensional scalar-valued NURBS basis function $R_k^{p_1 \cdot p_2}$ with a Cartesian base vector $\mathbf{e}_i$ is associated to a two-dimensional vector-valued NURBS basis function $\phi_r$ in the fashion $r \to (k, i)$, that is, $1 \to (1, 1)$, $2 \to (1, 2)$, $3 \to (1, 3)$, $4 \to (2, 1)$, etc. Additionally, $\hat{d}_r$ stands for the displacement of the $\lceil r/3 \rceil$-th CP along $r - (3\lceil r/3 \rceil - 3)$-Cartesian direction. Since the CPs do not interpolate the geometry, displacement components $\hat{d}_r$ have no direct physical meaning. By bringing the approximation given in Eq. (11) into the weak form of dynamic equilibrium defined by Eq. (6), the semi-discretized problem in terms of the dynamic residual vector $\mathbf{R}$ reads,

$$\bar{\mathbf{R}}(\hat{\mathbf{d}}) = \mathbf{M}\ddot{\hat{\mathbf{d}}} + \mathbf{D}\dot{\hat{\mathbf{d}}} + \mathbf{R}(\hat{\mathbf{d}}) - \hat{\mathbf{F}}_t, \quad (12)$$

where $\hat{\mathbf{d}} = [\hat{d}_1 \cdots \hat{d}_{3n}]^T$ stands for displacements of the CPs, $\dot{\hat{\mathbf{d}}}$ and $\ddot{\hat{\mathbf{d}}}$ being the corresponding discrete vectors of the velocity and acceleration Degrees of Freedom (DOFs). Moreover, the components of the mass matrix $\mathbf{M}$ and the steady-state residual vector $\mathbf{R}$ are given by

$$M_{ij} = \int_{\Omega} \rho_s t_s \phi_i \cdot \phi_j \, d\Omega, \quad (13a)$$

$$R_i = \int_{\Omega} \frac{\partial \boldsymbol{\varepsilon}}{\partial \hat{d}_i} : (\mathbf{n} + \mathbf{n}_0) \, d\Omega. \quad (13b)$$

The damping matrix $\mathbf{D}$ is approximated using the Rayleigh approach and its corresponding definition is postponed to Section 2.1.3. Additionally, $\hat{\mathbf{F}}_t$ in Eq. (12) stands for the external load vector at each time instance $t$. In this way the components of the steady-state stiffness matrix are as follows:

$$K_{ij} = \frac{\partial R_i}{\partial \hat{d}_j} = \int_{\Omega} \frac{\partial \boldsymbol{\varepsilon}}{\partial \hat{d}_i} : \frac{\partial \mathbf{n}}{\partial \hat{d}_j} + \frac{\partial^2 \boldsymbol{\varepsilon}}{\partial \hat{d}_i \partial \hat{d}_j} : (\mathbf{n} + \mathbf{n}_0) \, d\Omega. \quad (14)$$
### 2.1.3. Time-discretization

The Newmark method is applied for the time discretization of Eq. (12), as it is unconditionally stable and provides second-order accuracy (proven for linear problems). In this way, the time domain $T$ is discretized into a finite number of time steps $t_n \in T$. Furthermore, if $\hat{\mathbf{d}}(t_n) = \hat{\mathbf{d}}_n$ is the discrete vector of DOFs at time step $t_n$, given the Newmark parameters $\beta_n, \gamma_n \in \mathbb{R}_+^*$ the following approximations are chosen (Newmark, 1959):

$$ \hat{\mathbf{d}}_n = \hat{\mathbf{d}}_{n-1} + \Delta t \hat{\mathbf{d}}_{n-1} + (\Delta t)^2 \left( \frac{1}{2} - \beta_n \right) \ddot{\hat{\mathbf{d}}}_{n-1} + (\Delta t)^2 \beta_n \ddot{\hat{\mathbf{d}}}_n, \quad (15a) $$

$$ \dot{\hat{\mathbf{d}}}_n = \dot{\hat{\mathbf{d}}}_{n-1} + \Delta t (1 - \gamma_n) \ddot{\hat{\mathbf{d}}}_{n-1} + \Delta t \gamma_n \ddot{\hat{\mathbf{d}}}_n, \quad (15b) $$

$$ \bar{\mathbf{R}}(\hat{\mathbf{d}}_n) = \mathbf{M} \ddot{\hat{\mathbf{d}}}_n + \mathbf{D} \dot{\hat{\mathbf{d}}}_n + \mathbf{R}(\hat{\mathbf{d}}_n) - \hat{\mathbf{F}}_n. \quad (15c) $$

The Newton-Raphson linearization of the nonlinear system of equations is employed resulting in a sequence of iterations,

$$ \bar{\mathbf{K}}(\hat{\mathbf{d}}_{n,i}) \Delta_i \hat{\mathbf{d}}_n = -\bar{\mathbf{R}}(\hat{\mathbf{d}}_{n,i}), \quad (16) $$

where the solution increment is defined as $\Delta_i \hat{\mathbf{d}}_n = \hat{\mathbf{d}}_{n,i+1} - \hat{\mathbf{d}}_{n,i}$, in the $i$-th iteration. $\bar{\mathbf{K}}$ stands for the dynamic tangent stiffness matrix whose entries are given in terms of the steady-state tangent stiffness matrix $\mathbf{K}$, the mass and damping matrix components via

$$ \bar{\mathbf{K}}_{ij}(\hat{\mathbf{d}}_{n,i}) = \frac{\partial \bar{\mathbf{R}}_i}{\partial \hat{\mathbf{d}}_j} = \left( \frac{1}{\beta_n (\Delta t)^2} M_{ij} + \frac{\gamma_n}{\beta_n \Delta t} D_{ij} \right) + K_{ij}(\hat{\mathbf{d}}_{n,i}). \quad (17) $$

As aforementioned, the computation of the damping matrix $\mathbf{D}$ is herein based on the Rayleigh damping approach (Clough and Penzien, 1993; Thomsen et al., 2000): Given the Rayleigh damping parameters $\alpha_r, \beta_r \in \mathbb{R}_+^*$, $\mathbf{D}$ is approximated as follows

$$ \mathbf{D} = \alpha_r \mathbf{M} + \beta_r \mathbf{K}(\hat{\mathbf{d}}_{0,0}), \quad (18) $$

where $\hat{\mathbf{d}}_{0,0}$ stands for the initial condition of the displacement DOFs. Thus, the damping matrix is computed using the initial linear stiffness matrix of the structure.

### 2.2. Computational fluid dynamics

In order to predict the flow field, the three-dimensional Navier-Stokes equations for an incompressible fluid are solved based on a finite volume discretization, which is second-order accurate in space and time. Several methods are available to solve the pressure-velocity coupling problem: In case of a laminar flow a standard implicit SIMPLE scheme can be applied. However, when the flow becomes turbulent, a semi-implicit predictor-corrector scheme (projection method), which is more appropriate for the large-eddy simulation (LES) technique, is preferred (Breuer et al., 2012). The large scales of the turbulent flow are resolved directly, whereas the small scales are modeled by a subgrid-scale (SGS) model such as the standard Smagorinsky model (1963). In order to reasonably predict the flow applying LES, an appropriate representation of the inflow boundary conditions is required (Wood et al., 2016; Breuer, 2018). In the present LES framework turbulent perturbations mimicking the reality are synthetically generated by an inflow generator based on the digital filter concept of Klein et al. (2003). These data are then injected as source terms inside the domain (Schmidt and Breuer (2017), De Nayer et al. (2018b) and Breuer (2018)).

In this study, both the pressure and the shear stresses of the fluid are taken into account along the FSI interface for the computation of the fluid tractions. Based on the pressure $p$, the velocity $\mathbf{u}$ and the dynamic viscosity $\mu$ of the fluid, the fluid tractions are given by

$$ \mathbf{t} = -p \mathbf{e}_n + \mu \mathbf{e}_n \cdot \nabla_s \mathbf{u}. \quad (19) $$

In Eq. (19), $\mathbf{e}_n$ and $\nabla_s = 1/2(\nabla + \nabla^T)$ stand for the unit vector normal to the FSI surface and the symmetric gradient tensor, respectively.

In case of FSI applications the temporally varying domain is taken into account by the Arbitrary Lagrangian-Eulerian (ALE) formulation. Since the grid is body-fitted, a mesh adaption method is necessary: A very fast algebraic technique relying on a combination of linear and transfinite interpolations (TFI) is employed for FSI problems undergoing small displacements and deformations (Thompson et al., 1985). In the context of larger displacements a hybrid method especially developed for FSI-LES problems and block-structured grids (Sen et al., 2017) is preferred. For further details concerning the CFD technique used, the reader is referred to the cited references.

### 2.3. Isogeometric mortar-based mapping

In this section the isogeometric mortar-based mapping (Bazilevs et al., 2012) for field transformations between non-matching discretizations which share a common interface is introduced. The proposed methodology is an extension of
<figure><img src="image_5.png" /><figcaption>Fig. 2. Isogeometric mortar-based mapping: Problem definition.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Isogeometric mortar-based mapping: Projection of low-order elements onto NURBS-based CAD model.</figcaption></figure>

the mortar-based mapping method suggested in Wang et al. (2016) involving two low-order faceted discretizations to account also for isogeometric discretizations. Let the surface $\Omega$ be modeled using the NURBS basis functions as introduced in Section 2.1.2 with the parametrization $S$ (see Fig. 2(a) for example). Furthermore, let $\Omega_h$ be a low-order faceted approximation of $\Omega$ with parametrization $S_h$ (see Fig. 2(b)). In the present work, $\Omega$ is linked to the wet surface of the membrane and $\Omega_h$ to the fluid FSI interface. Since in the current study the fluid domain is discretized using finite volumes on block-structured grids (see Section 4), the surface $\Omega_h$ is represented by quadrilaterals. These can be split into triangles in order to describe $\Omega_h$ by a low-order representation based on linear triangles. The underlying parametrization of the surface $\Omega_h$ is done using the standard FE basis functions meaning that

$$S_h = \sum_{i=1}^{n_h} N_i \mathbf{X}_i, \qquad (20)$$

where $N_i$, $\mathbf{X}_i$ and $n_h$ stand for the low-order basis functions (linear for triangles or bilinear for quadrilaterals), the corresponding set of nodal coordinates and the total number of nodes in the low-order discretization $\Omega_h$, respectively. In this way, given an isogeometrically discretized vector field $\mathbf{q}$ defined in $\Omega$, which may represent either the interface displacement or traction field, a piecewise linear (or bilinear) vector field $\mathbf{q}^h$ on $\Omega_h$ is sought such that,

$$\mathbf{q}^h = \arg \min_{\bar{\mathbf{q}}^h \in \mathcal{L}^2(\Omega)} \|\mathbf{q} - \bar{\mathbf{q}}^h\|_{0,\Omega}, \qquad (21)$$

where $\mathcal{L}^2(\Omega)$ stands for the space of all square integrable vector functions in $\Omega$ and $\|\cdot\|_{0,\Omega}^1$ is its corresponding norm defined by

$$\|\mathbf{q}\|_{0,\Omega} = \left( \int_{\Omega} \mathbf{q} \cdot \mathbf{q} \, d\Omega \right)^{\frac{1}{2}}, \quad \text{for all } \mathbf{q} \in \mathcal{L}^2(\Omega). \qquad (22)$$

<footnote n="1">¹ Symbol • stands for an arbitrary variable.</footnote>
A common interface between the domains $\Omega$ and $\Omega_h$ has to be established for the problem defined in Eq. (21) to be meaningful. Therefore, the exact representation of surface $\Omega$ by means of NURBS is herein employed as the common integration area and subsequently the low-order surface discretization $\Omega_h$ is mapped onto $\Omega$. Consider the triangle $X_i-X_j-X_k$ in $\Omega_h$ (see Fig. 2(b)). Each of its nodes is projected onto $\Omega$ using a Newton-Raphson scheme, finding the corresponding local parameters of its vertices in $\hat{\Omega}$, namely, $\hat{\theta}_i = (\hat{\theta}_1^i, \hat{\theta}_2^i)$, $\hat{\theta}_j = (\hat{\theta}_1^j, \hat{\theta}_2^j)$ and $\hat{\theta}_k = (\hat{\theta}_1^k, \hat{\theta}_2^k)$, respectively (see Fig. 3(a)). The nodes are assumed to be linearly connected in $\hat{\Omega}$ forming the projection of the specific triangle onto $\Omega$. Subsequently, the image of the triangle transferred back to the physical space is depicted in Fig. 3(b) through its projected vertices $X_i^p, X_j^p$ and $X_k^p$, respectively, where $X_i^p = S(\hat{\theta}_i)$ stands for the corresponding geometric map. Having defined the common interface between the domains of the fields $\mathbf{q}$ and $\mathbf{q}^h$, the weak form of the problem in Eq. (21) states: Given a field $\mathbf{q} \in \mathcal{L}^2(\Omega_h)$, find a field $\mathbf{q}^h \in \mathcal{L}^2(\Omega_h)$ such that,

$$ \int_{\Omega} \delta \mathbf{q}^h \cdot \mathbf{q}^h \, d\Omega = \int_{\Omega} \delta \mathbf{q}^h \cdot \mathbf{q} \, d\Omega, \quad \text{for all } \delta \mathbf{q}^h \in \mathcal{L}^2(\Omega_h). \qquad (23) $$

Since the fields $\mathbf{q}$ and $\mathbf{q}^h$ are discretized by using isogeometric and “standard” finite element basis functions in $\Omega$ and $\Omega_h$, respectively, they can be expressed as

$$ \mathbf{q} = \sum_{r=1}^{3n} \phi_r \hat{q}_r, \qquad (24a) $$

$$ \mathbf{q}^h = \sum_{r=1}^{3n_h} \phi_r^h \hat{q}_r^h, \qquad (24b) $$

where $\phi_r^h = N_k e_i$ (with $k = 1, \dots, n_h$, $i = 1, \dots, 3$ and $r$ is defined as in Section 2.1.2) are the standard vector-values linear or bilinear basis functions. The definition of the DOFs is similar to the isogeometric discretization of the displacement field in Section 2.1.2. Herein, $\hat{d}_i$ and $\hat{d}_j$ are the discrete values of the fields defined on each CP and each node, respectively. This leads to the matrix form of the isogeometric mortar-based mapping, namely,

$$ \hat{\mathbf{q}}^h = (\mathbf{C}^{nn})^{-1} \mathbf{C}^{nr} \hat{\mathbf{q}}, \qquad (25) $$

where $\hat{\mathbf{q}}^h = [\hat{q}_1^h \dots \hat{q}_{3n}^h]^T$ and $\hat{\mathbf{q}} = [\hat{q}_1 \dots \hat{q}_{3n}]^T$ stand for the discrete vectors of DOFs for the finite element and the isogeometric discretization, respectively. The entries of the mortar matrices $\mathbf{C}^{nn}$ and $\mathbf{C}^{nr}$ in Eq. (25) are then given by

$$ C_{ij}^{nn} = \int_{\Omega} \phi_i^h \cdot \phi_j^h \, d\Omega, \qquad (26a) $$

$$ C_{ij}^{nr} = \int_{\Omega} \phi_i^h \cdot \phi_j^h \, d\Omega. \qquad (26b) $$

On the other hand, if the field $\mathbf{q}^h$ is given and the isogeometric field $\mathbf{q}$ is sought, then the problem in Eq. (23) reverses and the weak form accordingly becomes: Given a field $\mathbf{q}^h \in \mathcal{L}^2(\Omega_h)$, find a field $\mathbf{q} \in \mathcal{L}^2(\Omega) \cup \mathcal{L}^2(\Gamma_d)$ such that,

$$ \int_{\Omega} \delta \mathbf{q} \cdot \mathbf{q} \, d\Omega + \int_{\Gamma_d} \bar{\alpha} \delta \mathbf{q} \cdot \mathbf{q} \, d\Gamma = \int_{\Omega} \delta \mathbf{q} \cdot \mathbf{q}^h \, d\Omega, \quad \text{for all } \delta \mathbf{q} \in \mathcal{L}^2(\Omega) \cup \mathcal{L}^2(\Gamma_d). \qquad (27) $$

In Eq. (27) the field $\mathbf{q}$ is weakly set to zero along the boundary $\Gamma_d$ by means of the penalty method using a penalty parameter $\bar{\alpha} : \Gamma_d \to \mathbb{R}$. The latter is optional and can be selectively used in case a fixed value of the mapped field along the boundary is desired. Since the considered patches are untrimmed, it is possible herein to enforce the constraint along $\Gamma_d$ strongly. However, the weak enforcement is herein chosen, since it is more general and can be later on extended to account also for trimmed patches. In addition, the corresponding admissible space does not have to comply with the Dirichlet boundary conditions thus enabling a better approximation of the mapped field in the domain $\Omega$ (Babuska, 1973). Accordingly, the discrete system of equations becomes,

$$ \hat{\mathbf{q}} = (\mathbf{C}^{rr})^{-1} \mathbf{C}^{rn} \hat{\mathbf{q}}^h, \qquad (28) $$

where the entries of the underlying matrices accordingly read

$$ C_{ij}^{rr} = \int_{\Omega} \phi_i \cdot \phi_j \, d\Omega + \int_{\Gamma_d} \bar{\alpha} \phi_i \cdot \phi_j \, d\Gamma, \qquad (29a) $$

$$ C_{ij}^{rn} = C_{ji}^{nr} = \int_{\Omega} \phi_i^h \cdot \phi_j^h \, d\Omega. \qquad (29b) $$

Notably, $\mathbf{C}^{nn}$ and $\mathbf{C}^{rr}$ are square, symmetric and positive definite matrices rendering the mortar-based mapping well-posed. In addition, these matrices are sparse and thus the employed linear solvers benefit in terms of storage and efficiency. Moreover, since no coupling between the three Cartesian directions appears in these matrices, they can be reduced to
<figure><img src="image_7.png" /><figcaption>Fig. 4. Sequential isogeometric mortar-based mapping via the Exact Coupling Layer (ECL).</figcaption></figure>

account for each Cartesian direction separately, making the procedure even more efficient. The evaluation of the integrals in Eqs. (26) and (29) over $\Omega$ is performed at the subdomains $\Omega_i \subset \Omega$ with $i = 1, \dots, n_t \in \mathbb{N}$ (Fig. 3(a)) which emerge when clipping the projection of the low-order elements (into the NURBS parameter space) with the knots of the NURBS parametrization (as depicted in Fig. 3(a)). For this purpose, the Gauss integration method is used on a triangulation of each $\Omega_i$, where given a corresponding set of Gauss points $(\theta_{1,i}^{g,j}, \theta_{2,i}^{g,j})$ in the NURBS parametric space $\hat{\Omega}_i$ (Fig. 3(b)) of each $\Omega_i$ and Gauss weights $w_i^{g,j}$ with $j = 1, \dots, n_i^g$ one has,

$$ \int_{\Omega} \cdot d\Omega = \sum_{i=1}^{n_t} \int_{\Omega_i} \cdot d\Omega = \sum_{i=1}^{n_t} \int_{\hat{\Omega}_i} \cdot \|\mathbf{A}_3\|_2 d\hat{\Omega} = \sum_{i=1}^{n_t} \sum_{j=1}^{n_i^g} (\|\mathbf{A}_3\|_2)_{(\theta_{1,i}^{g,j}, \theta_{2,i}^{g,j})} w_i^{g,j}. \quad (30) $$

This ensures that the basis function products reach $C^\infty$-continuity within the integration subdomains. The low-order finite element basis functions are then computed at each Gauss point with respect to the shape of the low-order element in the NURBS parametric space. This means that a one-to-one correspondence between the original and the projected low-order elements is established.

As already mentioned in the introduction, the isogeometric mortar-based mapping method can also be used for smoothing the transformed fields taking advantage of the arbitrarily smooth high-order NURBS basis functions introduced in Eq. (7). This methodology is especially useful when fields from a coarsely discretized surface are mapped onto a finely discretized surface: To be more specific, let field $q_{\text{coarse}}^h$ be defined by a coarse low-order discretization and that its transformation onto a field $q_{\text{fine}}^h$ defined on a fine low-order discretization is sought. This can be realized by establishing the transformation via the smooth NURBS-based CAD model of the underlying geometry, the so-called Exact Coupling Layer (ECL): The isogeometric mortar-based mapping can be applied from the coarse low-order discretization onto the CAD surface and from the CAD surface onto the fine low-order discretization sequentially, that is,

$$ q_{\text{fine}}^{h} = ( (\mathbf{c}^{rr})^{-1} \mathbf{c}^{rn} ) ( (\mathbf{c}^{nn})^{-1} \mathbf{c}^{nr} ) q_{\text{coarse}}^{h} . \quad (31) $$

In this way, the artificial geometric roughness defining the coarse low-order discretization is smoothed through the ECL, thus enabling a smooth field representation on the fine low-order discretization (see Fig. 4). Typically, the coarse low-order discretization is associated with the structural discretization and the fine low-order discretization with the fluid FSI interface mesh. Thus, naturally the displacement field is transformed from the structural discretization to the fluid FSI interface mesh where the smoothing through the ECL proves highly beneficial, since kinks which may naturally occur on the displaced structural mesh are not transferred to the fluid FSI interface. Within this work, the transformation of the traction field from the fluid FSI interface to the structural mesh is also filtered through the ECL within the isogeometric mortar-based mapping method so that the methodology is kept consistent in both mapping directions.

## 2.4. FSI coupling

Based on the aforementioned mapping methodology, the structural displacement can be transformed to the fluid side and vice versa the fluid tractions can be transformed from the fluid FSI interface onto the structural side. Since the pressure and the velocity fields are defined at the cell center in the CFD problem, the fluid tractions (see Eq. (19)) are also computed at the center of the cell face which is associated with the fluid FSI interface. Since the traction is considered constant in each of these faces, it is discontinuous across neighboring faces along the fluid FSI interface. Thus, the traction vector at each node of the fluid FSI interface is linearly averaged from the neighboring elements. Then, a bilinear finite element interpolant at each face of the fluid FSI interface is applied to construct a continuous traction field. This is used within the proposed mortar-based mapping method in Section 2.3 to transform the tractions from the fluid to the structural FSI interface. Lastly, the consistent load vector of the corresponding fluid tractions on the structural FSI interface (which is a part of vector b in Eq. (6)) is constructed by integrating the transformed fluid tractions along the structural FSI interface.

In order to couple both the structural and fluid problems in time and thus to solve the coupled FSI problem, the present numerical methodology relies on a standard partitioned approach. If the implicit SIMPLE scheme (in case of laminar flows) is applied on the fluid side, the FSI coupling is carried out by an implicit iterative coupling strategy similar to the one described by Glück et al. (2003). With the semi-implicit predictor-corrector method used for the large-eddy simulations in case of turbulent flows, the adapted semi-implicit iterative FSI coupling published in Breuer et al. (2012) is preferred, since it significantly reduces the CPU-time.
Both strategies are based on one or several field exchange cycles (displacements and tractions) between the structural and the fluid problems during one time step which are called FSI-subiterations. In case of applications with low added-mass effect, where the ratio of the structural density to the fluid density is large, only one FSI-subiteration is necessary: The structural and the fluid problems are successively solved once per time step. In that case the procedure is called *loose or explicit coupling scheme* (Causin et al., 2005; Farhat et al., 2006; Wüchner et al., 2007). On the contrary, i.e., in case of a non-negligible added-mass effect, several FSI-subiterations are required to ensure the equilibrium along the FSI interface. The FSI-coupling is then denoted as *strong or implicit coupling scheme*. During this FSI-subcycle an underrelaxation of one of the fields (displacements or tractions) is often necessary to guarantee the FSI convergence in case of an unfavorable added-mass. The relaxation factor can be constant or dynamically computed by the Aitken algorithm (Aitken, 1926). The dynamically updated relaxation factor can considerably speed-up the FSI convergence (Mok, 2001). The coupled solution is considered converged if the following FSI criterion is reached: The 2-norm of the displacement difference<fnref n="2" /> between two FSI-subiterations (index $\hat{k}$) normalized by the 2-norm of the change in the displacements between the current and the last time step (index $\hat{n}$) drops below a predefined limit $\epsilon_{FSI}$, i.e.,

$$ \frac{\| \hat{\mathbf{d}}_{\hat{n},\hat{k}} |_{\Omega_{\text{FSI}}} - \hat{\mathbf{d}}_{\hat{n},\hat{k}-1} |_{\Omega_{\text{FSI}}} \|_2}{\| \hat{\mathbf{d}}_{\hat{n},\hat{k}} |_{\Omega_{\text{FSI}}} - \hat{\mathbf{d}}_{\hat{n}-1} |_{\Omega_{\text{FSI}}} \|_2} < \epsilon_{\text{FSI}}. \quad (32) $$

where $\hat{\mathbf{d}}_{\hat{n},\hat{k}} |_{\Omega_{\text{FSI}}}$ and $\hat{\mathbf{d}}_{\hat{n}-1} |_{\Omega_{\text{FSI}}}$ stand for the FSI interface displacement DOFs at the $\hat{k}$-th coupling iteration at time $t_{\hat{n}}$ and the converged displacement DOFs along the FSI interface $\Omega_{\text{FSI}}$ at time step $t_{\hat{n}-1}$, respectively. Moreover, the index $\hat{i}$ of the Newton-Raphson iteration is herein omitted since the aforementioned displacement DOFs of the converged nonlinear structural solution are considered. Typically, $\epsilon_{\text{FSI}}$ is set to $10^{-6}$ and the procedure needs about 5 FSI-subiterations to converge up to this accuracy.

The important role of the energy transfer within partitioned FSI simulations has been shown in Piperno and Farhat (2001) and Dervieux et al. (2010) among others. Therefore, this matter is also investigated herein. Let $W_{\Omega_{\text{FSI}}}$ be the work done by the fluid forces on the moving FSI interface $\Omega_{\text{FSI}}$ at a given time. This is defined as follows,

$$ W_{\Omega_{\text{FSI}}} = \int_{\Omega_{\text{FSI}}} \mathbf{t} \cdot \mathbf{d} \, d\Omega, \quad (33) $$

where $\mathbf{t}$ and $\mathbf{d}$ stand for the traction and displacement fields along the FSI interface $\Omega_{\text{FSI}}$ at a given time defined on either the fluid or the structural subdomain. The superscript $\bullet^{(f)}$ or $\bullet^{(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 Section 2.3), their discrete representations (see Eq. (24)) can be substituted in Eq. (33). 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:

$$ W_{\Omega_{\text{FSI}}}^{(f)} = (\hat{\mathbf{t}}^{(f)})^{\mathrm{T}} \mathbf{c}^{\mathrm{nn}} \hat{\mathbf{d}}^{(f)}, \quad (34a) $$

$$ W_{\Omega_{\text{FSI}}}^{(s)} = (\hat{\mathbf{t}}^{(s)})^{\mathrm{T}} \mathbf{c}^{\mathrm{rr}} \hat{\mathbf{d}}^{(s)}. \quad (34b) $$

If the ECL is applied, the structure is discretized by FEM and thus $W_{\Omega_{\text{FSI}}}^{(s)}$ is computed using Eq. (34a) using the corresponding structural $\mathbf{c}^{\mathrm{nn}}$ matrix. The residual interface energy $\mathcal{E}_{\Omega_{\text{FSI}}}$ is defined by the difference between the structural and the fluid interface work:

$$ \mathcal{E}_{\Omega_{\text{FSI}}} = W_{\Omega_{\text{FSI}}}^{(s)} - W_{\Omega_{\text{FSI}}}^{(f)}. \quad (35) $$

### 3. Lid-driven cavity with a flexible membrane attached at its bottom

#### 3.1. Case setup

The FSI cavity benchmark is employed for the evaluation of the proposed methodology as it offers a simple 2D geometry while relatively large deformations are allowed. This case was initially proposed by Wall (1999) and Mok (2001) and then used by many other groups as benchmark (see, e.g., Gerbeau and Vidrascu, 2003; Kassiotis et al., 2011). It involves a square domain where the top wall is impermeable and moving with a given time-dependent cosinusoidal velocity $\mathbf{u}_i$. The two impermeable side-walls are fixed, while a flexible impermeable membrane is attached at the bottom (see Fig. 5(a)). A small area on the left and the right side-walls is chosen as inlet and outlet, respectively. The previously described numerical methodology is applied combining a strong FSI coupling strategy with an FSI convergence criterion of $\epsilon_{\text{FSI}} = 10^{-6}$. Since the Reynolds number in the fluid domain does not exceed 200, the flow is considered as laminar and no turbulence model is applied. The kinematic viscosity of the fluid is set to $10^{-2} \text{ m}^2/\text{s}$. A membrane model is used for the flexible structure with Young's modulus $E = 250 \text{ Pa}$, Poisson's ratio $\nu = 0$, density $\rho_s = 500 \text{ kg/m}^3$ and thickness $t_s = 2 \times 10^{-3} \text{ m}$.

<footnote n="2">² The 2-norm of a discrete vector of DOFs is equal to the square root of the sum of its squared components.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 5. FSI cavity benchmark: Setup.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. FSI cavity benchmark: Refinement study using the time displacement curves of a node in the middle of the flexible membrane (FSI simulation with 100 structural FEM elements).</figcaption></figure>

Based on a preliminary study taking different time steps into account (see Fig. 6(a)), the finally applied time step size is $\Delta t = 10^{-2}$ s for both the fluid and structural problems. To ensure the spatial grid independency of the fluid solution, two equidistant Cartesian fluid grid ($30 \times 30$ and $50 \times 50$ control volumes (CVs)) were tested (see Fig. 6(b)). The $30 \times 30$ fluid mesh (see Fig. 5(b)) is found to be sufficient. For the structure multiple discretizations of the flexible membrane are considered. These discretizations involve standard quadrilateral finite element and isogeometric discretizations. Due to the simplicity of the geometry of the membrane, a B-Spline patch is chosen for the isogeometric discretization. The B-Spline surface is subsequently also used as the *Exact Coupling Layer*: The tractions mapped from the fluid to the structural domain as well as the displacements mapped from the structural to the fluid domain are filtered through the B-Spline surface, thus smoothing the underlying vector fields. The forthcoming simulations involve standard $C^0$-continuous FEM with 6 elements with and without the ECL and IGA with $10 C^1$-continuous elements. The latter NURBS parametrization, that is the $10 C^1$-continuous elements, is also employed for modeling the ECL. As a reference solution, the time-dependent solution of the problem using 100 quadrilateral equidistant finite elements for the discretization of the flexible membrane is employed using the standard mortar-based mapping method.

## 3.2. Results

A qualitative comparison of the streamlines and the deformation of the membrane at $t = 19$ s simulation time, where the maximum of the displacement in the periodic FSI cycle is reached, is shown in Fig. 7. Accordingly, the time history of the displacements of the membrane's mid-point is shown in Fig. 8(a) for the reference data and the solution obtained by using $10$ biquadratic $C^1$-continuous isogeometric elements. Both qualitatively and quantitatively the comparison shows a very good agreement between both simulations, i.e., when using IGA and the reference finite element mesh for the structural discretization. Furthermore, the simulation results coincide with the results of Kassiotis et al. (2011). Then, the
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** FSI cavity benchmark: Streamlines corresponding to the FSI simulation with the reference mesh (100 quadrilateral membrane finite elements) and with IGA single patch mesh (10 biquadratic C¹-continuous isogeometric membrane elements) at time $t = 19$ s.</figcaption></figure>

Exact Coupling Layer (ECL) concept is exploited by using a very coarse structural finite element mesh consisting of only 6 elements, while the previously employed isogeometric discretization is used as the *Exact Coupling Layer*. Note that the ECL is used for both displacement and traction field transformations in order to keep the framework consistent, i.e., the same methodology in both transformation directions. Fig. 9 provides a qualitative comparison between the streamlines and the corresponding deformation of the fluid interface mesh without and with the application of the *Exact Coupling Layer*. Note that the faceted interface visible in Fig. 9(a) is defined on the fluid FSI interface and is caused by the low resolution of the FEM structural mesh generating an artificial roughness at the fluid FSI interface. The time history of the displacements at the membrane's mid-point and that of the fluid interface without and with the *Exact Coupling Layer* are shown in Fig. 8(b). Subsequently, they are compared against the reference solution. Since the structure is discretized with only 6 elements, both solutions (FEM 6 elements without and with ECL) are deviating from the reference solution, as expected. However, they produce similar results. The latter demonstrates that the *Exact Coupling Layer* concept can be used for appropriately smoothing the solution when transferred from one low-order discretized interface to another, thus eliminating artificial interface roughness triggering spurious flow field solutions (see Section 5.1). Additionally, a quantitative convergence evaluation of the mapping results is presented. The displacement field $d$ defined on the B-Spline surface at $t = 19$ s of the FSI simulation is mapped using Eq. (25) onto the fluid FSI interface, while the fluid FSI interface is gradually refined using an equidistant finite element mesh (5, 10, 20, 40, 80, 160 finite elements). Similarly, the traction field $t$ defined over the fluid interface at the same time of the FSI simulation is mapped using Eq. (28) onto the NURBS patch, while this NURBS patch is also uniformly refined (24, 48, 192, 384, 768 and 1536 equidistant knots) using piecewise C¹-continuous quadratic elements in the length direction, whereas one linear element in the width direction is used, since the problem is two-dimensional and effects in the width direction need to be suppressed. The results are shown in Fig. 10, where $h$ stands for either the minimum finite element edge size in the considered finite element or the square root of the minimum element area size found in the isogeometric discretization, which is the area of the smallest knot span in the Cartesian space. A uniform convergence to the target field can be observed with a convergence of second-order for both mapping transformations (displacements and tractions). The target displacement and traction fields defined over the NURBS patch and the fluid interface, respectively, at $t = 19$ s and the displacement and traction fields mapped onto the fluid FSI interface and the NURBS patch, respectively, are depicted in Fig. 11 for a given refinement level. The original mesh has 30 equidistant elements, so that the mapping is comparable with the one used for the FSI cavity benchmark (the CFD finite volume mesh has $30 \times 30$ CVs). An excellent agreement can be observed which is also supported by the corresponding convergence graphs.

The last step is to evaluate the residual energy at the interface $\varepsilon_{\Omega_{FSI}}$ with the help of Eq. (35) (see Fig. 12). For the four different setups the absolute value of the residual energy at each time step starts from zero and increases during the initialization phase. Then, in the periodic phase, it stays below $10^{-5}$ N m for FEM-6 and even below $10^{-9}$ N m for FEM-100, IGA and FEM-6 with ECL. These errors are negligible compared to the values of the work of the fluid forces at the FSI interface $W_{\Omega_{FSI}}$ (between $3 \times 10^{-3}$ N m and $5 \times 10^{-3}$ N m). This demonstrates that a better representation of the FSI interface with the presented isogeometric mortar mapping leads to lower residual energies and thus lower errors at this coupling interface. Although $\varepsilon_{\Omega_{FSI}}$ is not equal to zero as it would be in an exactly energy-conservative scheme, the corresponding error remains small during the simulation.

Relying on the simple lid-driven cavity case by Mok (2001), the numerical FSI methodology presented in Section 2 is therefore successfully validated.
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** FSI cavity benchmark: Time history of the displacements.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** FSI cavity benchmark: Streamlines corresponding to the FSI simulations without and with the *Exact Coupling Layer* (ECL) at time $t = 19$ s.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** FSI cavity benchmark: Relative error in the mapping of displacements from IGA to FEM and tractions from FEM to IGA.</figcaption></figure>

## 4. Air-inflated flexible hemisphere in turbulent flow

In order to test the numerical FSI methodology including the isogeometric analysis described in Section 2, the wall-mounted air-inflated membrane immersed in a turbulent boundary layer (Wood et al., 2018; De Nayer et al., 2018a,c) is studied. In comparison with the rather simple lid-driven cavity test case the geometry is three-dimensional, the flow
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** FSI cavity benchmark: The 2-norm contours of the target fields (top) compared with their corresponding mapped fields (bottom).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** FSI cavity benchmark: Time evolution of the interface work in the fluid and solid subdomains and time evolution of the residual interface energy $\mathcal{E}_{\Omega_{FSI}}$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Air-inflated flexible hemisphere: Problem definition.</figcaption></figure>

regime is turbulent at a high Reynolds number and experimental reference data are available by Wood et al. (2018). Fig. 13 depicts the geometry in detail and Table 1 summarizes the operating conditions. Among the three cases considered in Wood et al. (2018) the case with the highest Reynolds number $Re = DU_\infty \rho_{air}/\mu_{air} = 100,000$ is selected, since the largest displacements and the strongest fluctuations of the membrane are observed for this setup. For more details concerning the choice of the material, the determination of the material parameters and the choice of the characteristics of the incoming boundary layer and the inner pressure, refer to Wood et al. (2018). The designation of the Cartesian axes as $X_1-X_2-X_3$ used in the prequel for the governing equations in continuum mechanics is shifted to x-y-z in the sequel since that is more common in computational fluid dynamics.

## 4.1. Setup of the test case

The current section describes the computational setup to solve the FSI problem using an isogeometric membrane structure. The corresponding simulations herein are denoted as FSI-LES-IGA for the FSI simulation using IGA structural
<table><caption>**Table 1**
Air-inflated flexible hemisphere: Operating conditions.</caption><thead><tr><th>Re</th><th>D<br>[m]</th><th>ρ<sub>air</sub><br>[kg/m<sup>3</sup>]</th><th>μ<sub>air</sub><br>[kg/(m s)]</th><th>U<sub>∞</sub><br>[m/s]</th><th>ρ<sub>s</sub><br>[kg/m<sup>3</sup>]</th><th>t<sub>s</sub><br>[m]</th><th>E<sub>silicone</sub><br>[Pa]</th><th>ν<sub>silicone</sub><br>[-]</th><th>Δp<sub>FSI</sub><br>[Pa]</th></tr></thead><tbody><tr><td>100,000</td><td>0.15</td><td>1.225</td><td>18.27 × 10<sup>-6</sup></td><td>10.24</td><td>1050</td><td>1.65 × 10<sup>-4</sup></td><td>7 × 10<sup>5</sup></td><td>0.45</td><td>43</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Air-inflated flexible hemisphere: Single patch modeling of the hemisphere.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Air-inflated flexible hemisphere: Single patch modeling of the hemisphere: Form-found geometries and error distribution for each refinement level.</figcaption></figure>

analysis for the CSD problem and FSI-LES-FEM with ECL concerning the FSI simulation using a standard FEM discretization for the CSD problem while the field exchange takes place through the *Exact Coupling Layer* (ECL) which is modeled via NURBS. Numerical results of the same FSI problem but computed with a “standard” finite element membrane discretization are used in the present investigation for comparison and the corresponding simulation is abbreviated with FSI-LES-FEM. The setup of the FSI-LES-FEM simulation with 1926 constant strain triangle elements can be found in De Nayer et al. (2018a).

### 4.1.1. Single patch isogeometric modeling of the flexible hemisphere

Although the hemispherical shape of the membrane belongs to the category of generic shapes, which is also a synclastic surface, its modeling using NURBS results in parametrizations which attain singularities either along its boundary or at its center. Furthermore, there may also exist parametric directions along which the continuity drops to $C^0$. In order to obtain a hemispherical shape which is at least $C^1$-continuous in its interior, the *Updated Reference Strategy* (URS) form-finding method is employed (see Bletzinger and Ramm (1999), Wüchner and Bletzinger (2005) and Apostolatos et al. (2018) for more information). Given the pressure difference of magnitude $\Delta p_{cal} = 43$ Pa between the inner and the outer domain (based on this value and without a flow field the structure reaches an almost hemispherical shape in the calibration phase of the experiments), the membrane thickness $t_s$ and the diameter of the hemisphere $D$, the calibration pre-stress of the membrane for an equilibrium in an almost hemispherical configuration is given by Barlow’s formula

$$n_{cal} = \frac{\Delta p_{cal} D}{4 t_s} = 9772.3 \text{ Pa.} \qquad (36)$$

The setting of the form-finding problem can be seen in Fig. 14(a). The parametrization with four singularities at the Dirichlet boundary $\Gamma_d$ is chosen. The corresponding determinant of the geometrical Jacobian is shown in Fig. 14(b). Starting from a flat reference configuration for the form-finding problem, a series of form-finding steps according to the URS algorithm are carried out applying the constant pressure difference and the pre-stress $n_{cal}$ as well as the zero displacement field along $\Gamma_d$.

For the form-finding analysis, three different isogeometric discretizations of the flat circular plate are chosen, one with 144 biquadratic elements, one with 625 biquadratic elements and one with 625 bicubic elements. The results of the form-finding analysis for the three aforementioned configurations are represented in Fig. 15, respectively. The difference
of the form-found shape to the exact hemispherical shape is quantified using the following error measure,

$$e = (\mathbf{X}_h - \mathbf{X}) \cdot \mathbf{A}_3, \qquad (37)$$

where $\mathbf{X}_h$, $\mathbf{X}$ and $\mathbf{A}_3$ stand for all material points of the form-found geometry, all material points of the exact hemisphere and the outward normal to the exact hemisphere, respectively. The isogeometric discretizations with 625 elements<fnref n="3" /> are one order of magnitude more accurate than the coarser parametrization.

In Appendix the three different isogeometric discretizations are evaluated based on a steady and a dynamic pure structural case. In these pure CSD cases and in the last FSI application described below, the initial pre-stress has to be reduced to the value $n_{membrane} = 7794.5$ Pa as assigned in the FEM model of De Nayer et al. (2018a), since the gravitational forces are taken into account by means of a dead load. Furthermore, the silicone rubber used in the experiments is characterized by the St. Venant–Kirchhoff material model. Its two parameters, the Young's modulus and the Poisson's ratio are set to the experimentally evaluated values $E_{silicone} = 7 \times 10^5$ Pa and $\nu_{silicone} \approx 0.45$, respectively. Based on this setup, the steady and dynamic pure CSD cases in Appendix are used to determine the influence of the grid resolution. It is found that the grid with 625 biquadratic elements is the best compromise between accuracy and CPU consumption. Therefore, it is used in the following FSI study.

For the FSI case described in Wood et al. (2018), the inner pressure is maintained constant in the experiment. This is modeled in the CSD solver by means of a follower pressure load acting on each element surface (Schweizerhof and Ramm, 1984). The pressure difference between the inner and outer domain is set to $\Delta p_{FSI} = \Delta p_{cal}$.

The CSD solver applies the standard second-order Newmark scheme with a dimensionless time step $\Delta t^* = \Delta t U_{\infty} / D = 3.317 \times 10^{-5}$. This small value is driven by the fluid solver, since both solvers must have the same time step size within the employed solution approach. In addition, the Rayleigh damping, necessary for modeling the structural damping of the hemisphere, is used. The mass- and stiffness-proportional parameters are set to $\alpha_t = 17.47$ and $\beta_t = 1.89 \times 10^{-4}$, as determined in De Nayer et al. (2018a).

### 4.1.2. Large-eddy simulation of the flow

The turbulent flow is predicted by the large-eddy simulation technique applying the standard Smagorinsky subgrid-scale model with $C_s = 0.1$ and Van Driest damping as used in the FSI-LES-FEM simulation of De Nayer et al. (2018a). The flow is only solved in the outer domain, since the constant inner pressure is modeled by a follower pressure load as mentioned above. The external CFD domain is a large hemispherical expansion with its origin at the center of the hemisphere and an outer diameter of 20 D. In De Nayer et al. (2018b) CFD results of the flow around the rigid hemisphere obtained on a fine grid with about $30 \times 10^6$ CVs and on a medium grid with about $4.3 \times 10^6$ CVs were compared. It shows that the medium grid provides sufficient accuracy and saves a huge amount of computational time. Therefore, this medium grid is selected for the present FSI simulation.

The boundary conditions are exactly the same as applied for the FSI-LES-FEM simulation. The bottom wall is a no-slip fixed wall and the flexible membrane is also considered as a no-slip smooth wall. The correct characteristics of the turbulent boundary layer (distributions, length and time scales) are imposed by turbulence data synthetically generated by an inflow generator. These inflow data are injected inside the domain at $x/D = -1.5$ by a source-term formulation (De Nayer et al., 2018b).

Concerning the spatial discretization, the convective fluxes are approximated by the flux blending method to stabilize the simulation: 5% first-order accurate upwind scheme is combined with 95% second-order accurate central scheme. The CFD grid and the explicit time integration technique used yield a small dimensionless time step size of $\Delta t^* = 3.317 \times 10^{-5}$ to ensure stable simulations with CFL numbers below unity.

### 4.1.3. FSI setup

In order to couple the structural and fluid solvers and exchange data at the common interface, the FSI coupling described in Section 2.4 is selected. Since the ratio of the structural density to the fluid density is large ($\rho_s / \rho_{air} \approx 857$), the added-mass effect is not overly pronounced in this case and therefore a loose FSI coupling scheme can be applied (Causin et al., 2005; Farhat et al., 2006).

Analogous to the FSI-LES-FEM simulation the fluid domain is initialized by a pure LES before starting to solve the FSI problem. Moreover, at the beginning of the FSI simulation the structure is stabilized by a higher Young's modulus ($7 \times 10^6$ Pa) to ensure a stable coupled simulation. Then, $E_{silicone}$ is decreased to its correct value of $7 \times 10^5$ Pa in a time interval of 5000 time steps.

In the following two additional simulations are considered and compared with FSI-LES-FEM: An FSI simulation with an isogeometric discretization of the hemisphere (FSI-LES-IGA) and an FSI simulation with a “standard” triangular FEM computational model (FEM 1926 used in De Nayer et al. (2018a)) for the hemisphere combined with ECL (FSI-LES-FEM with ECL). The isogeometric discretization with 625 biquadratic $C^1$-continuous elements (see Fig. 15(b)) is applied for the hemisphere and for the ECL. Concerning FSI-LES-IGA, the fluid mesh at the FSI interface and its projection on the NURBS hemisphere are shown in Fig. 16(a) and Fig. 16(b), respectively. Concerning FSI-LES-FEM with ECL, the structural mesh and its projection on the NURBS hemisphere are shown in Fig. 16(c) and Fig. 16(d), respectively. For the transformation of

<footnote n="3">³ An element in NURBS-based IGA is the knot span of the underlying NURBS parametrization.</footnote>
<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Air-inflated flexible hemisphere: Low-order surface discretizations and their projections on the NURBS-parametrized hemisphere (Fig. 15(b)) along with the generated integration subdomains (triangles) in x-y view.</figcaption></figure>

the interface fields (displacements and tractions) between the fluid interface mesh and the FEM computational structural model the mortar-based mapping method described in Section 2.3 is applied. Both simulations are carried out on 40 processors for the CFD part, 1 processor for the CSD part and 1 processor for the coupling part.

## 4.2. Results

The coupled FSI problem of the turbulent flow around the flexible air-inflated membrane is investigated in detail, both experimentally by Wood et al. (2018) and numerically by De Nayer et al. (2018b). Fig. 17 depicts a snapshot of the flow field (pressure and vorticity magnitude) and the deformed membrane predicted by the FSI-LES-IGA simulation. The approaching turbulent boundary layer excites the membranous structure. The flow separates close to the apex of the hemisphere. The oscillations of the structure in this region are the result of the separation and roll-up of the shear layer arising from the Kelvin-Helmholtz instability. Different vortex shedding processes are found in the wake of the hemisphere: An asymmetric von Kármán vortex street irregularly alternating with an arc-type symmetric vortex shedding, both with different characteristic frequencies. Small vortices merging in the shear layer form medium-sized vortices which move downstream. The succession of the corresponding pressure fluctuations induces a series of indentations and elevations on the membrane. A wave of deformations is generated on the rear side due to the downward movement of the vortices. Furthermore, several eigenfrequencies of the hemispherical membrane are found in the experimental and numerical results.

In order to evaluate the turbulent flow around the flexible hemisphere predicted by the FSI-LES-IGA simulation, the results are averaged over a time interval of $\Delta T^* = 125$. Fig. 18 compares the time-averaged deformation of the flexible hemisphere obtained by the FSI-LES-IGA simulation with the FSI-LES-FEM results (with and without ECL) and the experimental data. The comparison takes place at the symmetry plane $y/D = 0$ in the wake region, since this location is of major importance for the generation of the recirculation area. The flow deflected by the bluff body accelerates in the
<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Air-inflated flexible hemisphere: Snapshot of the instantaneous flow predicted by the FSI-LES-IGA simulation (top: pressure; bottom: vorticity magnitude) in the symmetry plane. A magnification factor of 3 is applied to the displacement and added to the reference configuration to highlight the differences in shape.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Air-inflated flexible hemisphere: Comparison of the time-averaged deformation of the flexible hemisphere in the symmetry plane of the wake region. The experimental data and the numerical results of the FSI-LES-FEM simulation can be found in Wood et al. (2018) and De Nayer et al. (2018a), respectively. A magnification factor of a = 3 is applied to the average displacements $\Delta x/D$ and $\Delta z/D$ and added to the reference configuration to highlight the differences in shape.</figcaption></figure>

region above the apex. It induces a low pressure zone on the top (see Fig. 17), i.e., a vertical lift force, leading to a vertical displacement of the membrane. The time-averaged deformation at the apex predicted by FSI-LES-FEM with and without ECL is in very good agreement with the experimental data, but the FSI-LES-IGA results exhibit an even better agreement with the experimental data. Moreover, the shape of the FSI-LES-IGA deformations on the rear side tends towards the experimental data.

For the evaluation of the present FSI case the fluctuations of the displacements $\Delta x/D$ and $\Delta z/D$ are also relevant. They are compared based on the dimensionless standard deviations $S_x^* = S_x/D$ and $S_z^* = S_z/D$ and displayed as a function of the angle $\varphi$ in Fig. 19. The shape of $S_x^*$ and $S_z^*$ obtained by both FSI-LES-FEM and FSI-LES-IGA simulations
<figure><img src="image_22.png" /><figcaption>Fig. 19. Air-inflated flexible hemisphere: Standard deviation $S_x^*$ and $S_z^*$ of the corresponding displacements.</figcaption></figure>

<table><caption>**Table 2**

Time-averaged transition and separation locations in the midplane.</caption><thead><tr><th rowspan="2">Case</th><th colspan="2">Transition</th><th colspan="2">Separation</th></tr><tr><th>x<sub>tra</sub>/D</th><th>z<sub>tra</sub>/D</th><th>x<sub>sep</sub>/D</th><th>z<sub>sep</sub>/D</th></tr></thead><tbody><tr><td>FSI-LES-FEM</td><td>-1.26 &times; 10<sup>-2</sup></td><td>5.19 &times; 10<sup>-1</sup></td><td>4.96 &times; 10<sup>-2</sup></td><td>5.08 &times; 10<sup>-1</sup></td></tr><tr><td>FSI-LES-IGA</td><td>-5.72 &times; 10<sup>-3</sup></td><td>5.19 &times; 10<sup>-1</sup></td><td>5.57 &times; 10<sup>-2</sup></td><td>5.20 &times; 10<sup>-1</sup></td></tr></tbody></table>

are qualitatively similar but with deviations of the amplitude as explained below: Near the apex ($\varphi \approx 90^{\circ}$) the standard deviations are low and nearly constant. Then, they increase reaching a maximum and abruptly decrease to zero at the fixed boundary ($\varphi = 0^{\circ}$). Whereas the locations of the minimum and maximum of the predicted shape of $S_x^*$ are similar to the experimental curve, a discrepancy appears for $S_z^*$: A local peak present in the experimental data around $\varphi = 70^{\circ}$ does not exist in the simulations, neither in the previous FSI-LES-FEM nor in the new FSI-LES-IGA results. The most visible difference between both simulation results is found for the amplitude of the standard deviations: The FSI-LES-IGA setup predicts lower values fitting better to the experiments in the lower part of the hemisphere. The simulation FSI-LES-FEM with ECL shows almost the same behavior as the FSI-LES-FEM simulation. This is somehow expected since the smoothing induced by ECL has no significant effect on an FSI interface under the current condition of small deformations. Nevertheless, the ECL produces the expected results, thus successfully demonstrating its application to real-world problems. Conclusively, all performed simulations produce satisfactory results compared to the experiment, thus demonstrating that the structural analysis based on IGA and the ECL method can be reliably used in the context of FSI.

In order to compare the time-averaged flow, the statistics in the midplane $y/D = 0$ are once more selected. Fig. 20 provides a comparison of the first- and second-order moments predicted by FSI-LES-IGA with the previous FSI-LES-FEM results and the experimental data. As aforementioned the ECL method has no significant effect on the deformation of the FSI interface, which is only slightly deformed in this case. Therefore, the results obtained by FSI-LES-FEM with ECL are very similar to those predicted by FSI-LES-FEM and are not shown here. Additionally, the time-averaged transition<fnref n="4" /> and separation<fnref n="5" /> locations are summarized in Table 2. In order to determine the position of transition, a simple method relying on the Reynolds shear stress $\overline{u'w'}/U_{\infty}^2$ is applied, which can be considered as a quantity defining the exchange of momentum into the boundary layer. The threshold for $\overline{u'w'}/U_{\infty}^2$ is set to -0.001 as used in Schmidt and Breuer (2014).

The flow characteristics predicted by both simulations are globally similar. The approaching flow in front of the obstacle remains nearly identical, since the flexible structure has a minor impact on it. The flow accelerates over the front part of the hemisphere, becomes turbulent and detaches near the apex. The transition location of the present FSI-LES-IGA case moves slightly downstream in comparison to FSI-LES-FEM. This observation can be attributed to the better representation of the FSI interface on the fluid side due to IGA. A smoother surface leads to a later transition from the laminar to the turbulent flow regime. Similar to the transition, the separation location shifts downstream, but it also slightly moves in vertical direction. This change is rather small, implying no major difference in the size of the recirculation area. However, as depicted in Fig. 20(a) the global shape of the recirculation area predicted by FSI-LES-IGA fits better to the experimental data. The same is true for the reattachment point. The wake does not strongly differ, except concerning the quasi horizontal dividing streamline which is found at a higher position in case of FSI-LES-IGA tending towards the experiments. Moreover, the region of high Reynolds stresses of the shear layer (visible in Figs. 20(b)–20(d)) tends to decrease for the FSI-LES-IGA simulation in agreement with the experiments. It correlates with the lower standard deviations of the displacements observed with FSI-LES-IGA in Fig. 19. The comparison of both simulation results with the experimental data shows a global improvement of the flow field predicted by the FSI-LES-IGA setup.

<footnote n="4">⁴ The transition point is the location where the laminar flow becomes turbulent.</footnote>

<footnote n="5">⁵ The separation point or line is defined where the flow velocity gradient in wall-normal direction is equal to zero.</footnote>
<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Air-inflated flexible hemisphere: Comparison of the time-averaged flow characteristics at $Re = 100,000$ in the symmetry x-z-plane at $y/D = 0$ based on the experiments of Wood et al. (2018) (top), the results of the FSI-LES-FEM simulation of De Nayer et al. (2018a) (middle) and the current FSI-LES-IGA simulation (bottom).</figcaption></figure>

The residual energy $\varepsilon_{\Omega_{\text{FSI}}}$ at the FSI interface defined in Eq. (35) against the time is shown in Fig. 21. Additionally, the time evolution of the work of the fluid forces at this interface $W_{\Omega_{\text{FSI}}}$ is plotted. These results are gathered for a given dimensionless time interval of $\Delta T^* = 4$, but not at the same simulation time spans. Moreover, the flow around the air-inflated flexible hemisphere is turbulent, adding a chaotic contribution to the instantaneous results. Therefore, a direct quantitative comparison of the curves is not possible. For the setups FSI-LES-FEM without and with ECL the residual energy curves are almost identical, which is due to the small displacements observed in this case. The absolute value of $\varepsilon_{\Omega_{\text{FSI}}}$ stays around $10^{-6}$ N m, which is small compared to the work of the fluid forces. Concerning the FSI-LES-IGA setup the associated residual energy $\varepsilon_{\Omega_{\text{FSI}}}$ is three orders of magnitude lower (around $10^{-9}$ N m). Similar to the FSI cavity
<figure><img src="image_24.png" /><figcaption>Fig. 21. Air-inflated flexible hemisphere: Time evolution of the interface work in the fluid and solid subdomains and time evolution of the residual interface energy $E_{ΩFSI}$.</figcaption></figure>

benchmark, a better representation of the FSI interface based on IGA and the proposed isogeometric mortar mapping leads to lower errors at the coupling interface.

Summarizing, this FSI application shows the relevance of the applied numerical methodology for 3D flexible membrane structures. The smooth representation of the FSI interface provided by NURBS allows for an improvement of the static and dynamic behavior compared to the case where standard FEM discretizations are used.

## 5. Influence of the artificial roughness of the deformed interface

In an FSI simulation based on non-matching grids, a mapping procedure is required to ensure a correct exchange of the displacement and traction fields. The mapping between both meshes implies that the quality of the FSI interface is driven by the coarser grid. Contrary to the previous investigation on the flexible hemisphere using IGA, in the case where the discretization is based on FEM as in De Nayer et al. (2018a), the coarser grid is usually the structural one. In such a case, when the structure is discretized using a coarse finite element mesh but only subject to small displacements (as observed in the flexible hemisphere application), the interface quality is not overly degraded and has generally only a moderate influence on the results. However, if the structure undergoes large displacements, a low resolution of the FSI interface due to the FEM mesh is particularly critical on the fluid solution. Facets, corresponding to the structural elements, break through (see Fig. 22(c)). The surface is not smooth anymore and an artificial roughness appears on the fluid side, which is of particular importance for eddy-resolving simulation methodologies such as LES. In order to investigate the influence of this artificial roughness on the predicted flow field, a test case derived from the wall-mounted flexible hemisphere is studied.

### 5.1. Case setup

Initially, the geometry of the problem is exactly a rigid hemisphere (De Nayer et al., 2018b) (see Fig. 22(a)), while the flow field is initially predicted by a normal large-eddy simulation. Then, the prediction is continued by a one-way coupled simulation<fnref n="6" /> based on the partitioned approach described in Section 2. A significantly large displacement is progressively imposed to the hemispherical shape on the structural side within a time interval of 1000 time steps. The final form represents the 8th eigenmode (see Fig. 22(b)) of the linear modal analysis carried out in De Nayer et al. (2018a). After reaching its final shape, the large-eddy simulation continues as a pure CFD prediction. Prior to any analysis the CFD simulation is performed for a dimensionless time interval of $\Delta T^* = 20$ in order to reach a new statistically steady state.

Four different set-ups are taken into account: The first one relies on IGA for the structural discretization (see Section 4.1.1) and is considered as the reference. In the three other cases, standard finite element discretizations are used for the structure. Three of the unstructured meshes used in De Nayer et al. (2018a) are tested: A coarse one with 908 FE, an intermediate with 1926 FE and a fine one with 4080 FE.

The artificial roughness seen by the fluid solver can be evaluated using the FSI interface from the simulation with the isogeometric structural discretization as the reference. The arithmetic average roughness $R_a$ and the root-mean-square roughness $R_q$ are computed as defined in DIN EN ISO 4287:2010-07:

$$R_a = \frac{1}{n} \sum_{i=1}^{n} |r_i^{\text{FEM}} - r_i^{\text{IGA}}| \quad \text{and} \quad R_q = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (r_i^{\text{FEM}} - r_i^{\text{IGA}})^2}, \qquad (38)$$

<footnote n="6">⁶ One-way coupled herein means that the displacement is prescribed on the hemisphere without interaction with the fluid flow.</footnote>
<figure><img src="image_25.png" /><figcaption>**Fig. 22.** Rigid deformed hemisphere: Initial and deformed final shape of the FSI interface.</figcaption></figure>

<caption>**Table 3**</caption>

Rigid deformed hemisphere: Artificial roughness seen on the fluid side generated by the mapping of the structure side. The reference is the IGA FSI interface.

<table><caption>**Table 4**</caption><thead><tr><th></th><th>FEM 908</th><th>FEM 1926</th><th>FEM 4080</th></tr></thead><tbody><tr><th>R<sub>a</sub> [µm]</th><td>161</td><td>87</td><td>45</td></tr><tr><th>R<sub>q</sub> [µm]</th><td>210</td><td>111</td><td>58</td></tr></tbody></table>

Rigid deformed hemisphere: Location of the reattachment point in the wake. The reference is the IGA case.

<table><thead><tr><th></th><th>FEM 908</th><th>FEM 1926</th><th>FEM 4080</th><th>IGA</th></tr></thead><tbody><tr><th>x<sub>reattach</sub>/D</th><td>9.54 × 10<sup>-1</sup></td><td>8.21 × 10<sup>-1</sup></td><td>8.21 × 10<sup>-1</sup></td><td>7.92 × 10<sup>-1</sup></td></tr><tr><th>Relative error [%]</th><td>20.5</td><td>3.7</td><td>3.7</td><td>0</td></tr></tbody></table>

where $r_i$ is the position of the i-th node of the fluid grid at the FSI interface and n the total number of fluid nodes at the interface. Both formulas evaluate the roughness based on the grooves and the peaks of a rough surface (here FEM) in comparison with a smooth one (herein IGA). The definition of $R_q$ is similar to $R_a$. However, $R_q$ is more sensitive to the grooves and peaks than $R_a$. The grid applied in the fluid solver is the same for all cases considered. The medium grid (4.3 × 10⁶ CVs) is used as in Section 4. The number of nodes on the interface n is equal to 35,280. Table 3 summarizes these artificial roughness values found on the surface of the deformed hemisphere for the three simulations utilizing different finite element meshes.

As expected, both roughnesses $R_a$ and $R_q$ decrease by increasing the grid resolution. Due to its definition, $R_q$ generates higher values emphasizing even more the deviations of the deformed FSI interface from the IGA reference shape.

In order to quantify the impact of this artificial roughness on the flow field, data are gathered and time-averaged during a dimensionless time interval of $\Delta T^* \ge 500$. Two planes are selected to compare the results: The x-z-plane, i.e., the symmetry plane at $y/D = 0$ and the x-y-plane at $z/D = 0.35$, where the deformation of the shape is the most prominent. In these planes some relevant quantities are compared and the results obtained by IGA are considered as the reference.

## 5.2. Results and discussion

Fig. 23 represents the first-order moments, i.e., the time-averaged streamwise $\bar{u}/U_\infty$ and vertical $\bar{w}/U_\infty$ velocities. Significant differences between the reference solution and the solution obtained using FEM 908 case are observed. The artificial roughness on the top strongly influences the location of the separation and thus the recirculation area and the wake. Table 4 provides a comparison of the location of the reattachment point<fnref n="7" /> for each case using the IGA case as reference. The FEM 908 case leads to an error of more than 20% on this flow characteristic. The first-order moments in the wake obtained for the FEM 1926 and 4080 cases nicely converge to the IGA reference data. These observations are supported by the second-order moments in the same plane in Fig. 24 and 25, particularly for the FEM 4080 case. As expected, the changes in the Reynolds stresses show that an increase of the artificial roughness on the interface leads to a stronger mixing and thus a strong overestimation of the Reynolds stresses in the separated shear layer.

The flow results in the horizontal plane at $z/D = 0.35$ visible in Fig. 27 show even more discrepancies between the low resolution case (FEM 908) and the high resolution ones (FEM 4080, IGA): For the former the flow detaches earlier and the

<footnote n="7">⁷ The reattachment point is located in the wake where the dividing streamline of the recirculation area attaches to the wall again (saddle point visible in Fig. 26).</footnote>
<figure><img src="image_26.png" /><figcaption>**Fig. 23.** Rigid deformed hemisphere: Time-averaged first-order moments in the symmetry plane $y/D = 0$.</figcaption></figure>

recirculation areas are consequently larger (see Table 4). Moreover, an asymmetry in the results appears. This asymmetry is also visible in the time-averaged streamlines on the wall and on the deformed structure depicted in Fig. 26. The initial hemispherical shape and the displacements applied are symmetric with respect to $y = 0$. However, the mapping combined with an unstructured and therefore not exactly symmetric finite element mesh results in a small asymmetry in the final deformed shape. This leads to slightly different local values of the artificial roughness, which generate a perceptible flow asymmetry in the wake. It has to be noted that this asymmetry is not due to a too short time-averaging process, since the results are time-averaged over a dimensionless time interval of $\Delta T^* \ge 500$.

To conclude, the mapping between non-matching grids forms an FSI interface, whose resolution corresponds to the coarser grid, which is typically used on the structure side. In case a FEM mesh is used for the discretization of the structure, the interface is represented by the structural elements as facets implying an artificial roughness. This roughness increases with the decrease of the grid resolution and with the increase of the structural displacement. In case of large displacements the effects of this artificial roughness on the flow field are not negligible: Separation and reattachment locations are significantly influenced and expected symmetries of the flow field are lost. Obviously, IGA is of invaluable benefit here since it automatically guarantees a smooth interface.

## 6. Computational costs related to the interface description

The FSI methodology in the present study (i.e., Sections 2–4) is based on a standard partitioned approach. Therefore, the total CPU consumption for an FSI simulation is the sum of the computational time of the CFD problem, the computational time of the CSD problem and the computational time of the mapping and transfer between the CFD and CSD solvers. For each case the CFD problem is solved using a parallelized and well-optimized solver. The CSD problem is solved serially. The mapping and data transfer are also done serially, since the sizes of the different problems under consideration are small enough. The computational costs related to the FSI cavity benchmark are not provided here because the size of the corresponding coupled problem is very small and the CPU-time needed for each problem and for one time-step is under 10 ms. Thus, the cavity FSI benchmark cannot be used for the evaluation of the proposed methodology in terms of CPU resources and computational time. Therefore, the computational cost of the presented methodology is evaluated based on the flexible hemisphere case. The related simulations are carried out on a cluster with Intel Xeon Ivy Bridge E5-2650 v2 processors operating at 2.60 GHz coupled with an Infiniband QDR network.
<figure><img src="image_27.png" /><figcaption>Fig. 24. Rigid deformed hemisphere: Time-averaged second-order moments in the symmetry plane $y/D = 0$.</figcaption></figure>

<caption>**Table 5**</caption>

Comparison of the required CPU-times for the different setups of the flexible hemisphere case of Section 4. The values are given in seconds and in percent of the total CPU-time of one time-step.

<table><thead><tr><th rowspan="3">Simulation</th><th>Initialization</th><th colspan="6">Averaged CPU-time for 1 time-step</th></tr><tr><th rowspan="2">Mapping</th><th colspan="2">CFD</th><th colspan="2">CSD</th><th colspan="2">Coupling (Transfer)</th></tr><tr><th>[s]</th><th>[s]</th><th>[%]</th><th>[s]</th><th>[%]</th><th>[s]</th><th>[%]</th></tr></thead><tbody><tr><td>FSI-LES-FEM</td><td>8.0</td><td>34.4</td><td>99.3</td><td>0.15</td><td>0.4</td><td>&lt;0.10</td><td>&lt;0.3</td></tr><tr><td>FSI-LES-IGA</td><td>56.0</td><td>34.2</td><td>92.8</td><td>2.40</td><td>6.5</td><td>&lt;0.25</td><td>&lt;0.7</td></tr><tr><td>FSI-LES-FEM with ECL</td><td>59.0</td><td>34.2</td><td>98.1</td><td>0.15</td><td>0.4</td><td>&lt;0.50</td><td>&lt;1.5</td></tr></tbody></table>

Table 5 summarizes the CPU-time required for the initialization phase of the mapping and the CPU-time, which is consumed by each sub-problem during a typical time-step for each setup. Since the flexible hemisphere case is simulated with a loose FSI coupling algorithm, only one FSI-iteration is carried out at each time-step, leading to a quasi constant computational resources for the CFD problem.

As expected for an FSI case, the CPU-time consumed by the fluid problem is significantly larger than the CPU-time consumed by the structural problem for all simulations considered. However, the computational resources required by the IGA structural solver are not negligible compared to the FEM structural solver. This can be attributed to the use of high-order NURBS basis functions (see Section 2.1.2).

The initialization phase of the mapping comprises MPI data exchanges between the CFD and CSD solvers and the steps of the algorithmic procedure described in Section 2.3 except Eq. (28) or Eq. (31). The transfer step for the coupling part comprises the transformation of the field $q$ (displacement or traction field) in both directions (see Eq. (28) in case of FSI simulations where the structural and the fluid problems interact directly or Eq. (31) in case of FSI simulations using the ECL). The given computational time of this transfer step includes also the MPI data exchange between the solvers.

The computational cost of the initial mapping phase is negligible compared to the total time of the full FSI simulation, which requires more than $10^6$ time steps. The computational part related to the transfer of the fields (displacements and tractions) between the domains depends on the problem size and the chosen discretization. Within FSI-LES-IGA, the amount of data is typically larger leading to larger coupling matrices ($C^{\text{nr}}$, $C^{\text{nr}}$ and $C^{\text{rr}}$) and therefore higher CPU time consumption. As it can be seen in Eq. (31), the application of the ECL adds supplementary costs at each transfer:
<figure><img src="image_28.png" /><figcaption>Fig. 25. Rigid deformed hemisphere: Time-averaged second moments in the symmetry plane $y/D = 0$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 26. Rigid deformed hemisphere: Time-averaged streamlines near the bottom wall and on the surface of the rigid deformed hemisphere.</figcaption></figure>

Indeed, two successive transformations are required for one field leading to doubled computational costs compared to the standard FSI-LES-IGA.

Although the use of the ECL introduces additional computational costs at each transfer, the solution of the whole coupled problem with FEM using ECL is generally much faster than solving the same FSI problem with IGA. Therefore, combining an available FEM structural discretization with the ECL is an attractive concept for engineers in term of efficiency. It is worth noting that the FSI simulations using a structural solver based on IGA have the big advantage that no mesh generation is needed as the model can be directly imported from CAD (Teschemacher et al., 2018). On the contrary for the structural solvers based on FEM the costly mesh generation step is unavoidable.

The influence of the FSI interface description on the computational costs of the fluid solver is worth noting. The average CPU-time values for one time-step of the CFD problem presented in Section 5 are summarized in Table 6. It can be deduced
FEM
908

<figure><img src="image_30.png" /></figure>

FEM
1926

<figure><img src="image_31.png" /></figure>

FEM
4080

<figure><img src="image_32.png" /></figure>

IGA

<figure><img src="image_33.png" /></figure>

(a) $\bar{u}/U_\infty$

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

(b) $\overline{u'u'}/U_\infty^2$

**Fig. 27.** Rigid deformed hemisphere: Time-averaged streamwise velocity and Reynolds stress in the horizontal plane $z/D = 0.35$ (same contours as in Figs. 23 and 24).

that the computational costs are reduced when a better interface description is used: The fluid solver converges faster and the CPU-time is reduced by 20% when a smooth interface description (IGA) is used contrary to simulation FEM-908, where the interface description has a higher artificial roughness due to the faceted discretization.
<table><caption>**Table 6**
Comparison of the required CPU-time between the different setups of the rigid deformed hemisphere of Section 5. The values are given in seconds.</caption><thead><tr><td>Simulation</td><td>Averaged CPU-time for 1 time-step [s]</td></tr></thead><tbody><tr><td>FEM-908</td><td>30.3</td></tr><tr><td>FEM-1926</td><td>29.0</td></tr><tr><td>FEM-4080</td><td>28.2</td></tr><tr><td>IGA</td><td>24.1</td></tr></tbody></table>

## 7. Conclusions and outlook

In the context of FSI simulations involving complex geometries an exact representation of the FSI interface is of great importance. Recently, with the introduction of the isogeometric analysis in computational structural mechanics, the exact description of the geometry from a CAD software can be directly used as basis for the structural solution. To allow the coupling between structures discretized using IGA and fluid flows discretized with the finite volume scheme, an isogeometric mortar-based mapping method is elaborated and evaluated in the present study. It enables FSI simulations with isogeometrically discretized structures and hence a smooth representation of the FSI interface.

At first the complete numerical methodology for solving fluid-structure interaction problems involving thin-walled membranous structures is shown based on a partitioned solution strategy relying on the ALE approach and the LES technique for turbulent flows. An introduction to membrane structural analysis is provided, while subsequently the application of IGA for the discretization of this kind of problems is shown. In order to accurately solve the FSI problem, the displacement and traction fields need to be properly exchanged at the interface between both domains, whereas the FSI interface has to be accurately represented. This is done by the isogeometric mortar-based mapping method, which is the core of the present study and enables a rigorous evaluation of the influence of the interface description. The variational form of the isogeometric mortar-based mapping method is described along with the corresponding discrete system of equations and implementation aspects. Firstly, its application to FSI simulations involving isogeometrically discretized structures, where the displacement and traction fields need to be transformed between a low-order discretized surface and a surface which is parametrized using high-order NURBS, is demonstrated. Subsequently, the *Exact Coupling Layer* concept is presented: The ECL method allows for a smooth transformation of the displacement and traction fields across the FSI interface between two low-order discretizations by taking advantage of the smoothing properties of the isogeometric mortar-based mapping technique. This is especially beneficial when the structure is discretized with significantly less elements than the fluid interface for, i.e., FSI simulations involving turbulent flows. The ECL concept prevents kinks resulting from the structural deformation, which is highly advantageous for the CFD solution.

The aforementioned applications are firstly shown for the simple cavity FSI benchmark, where the consistency of the method is validated. A single $C^1$-continuous membrane patch is employed. The results are compared against the standard mortar-based method, when a standard finite element discretization of the membrane is used demonstrating an excellent agreement. Subsequently, the effects of smoothing the transferred fields with the ECL concept is highlighted based on the same FSI benchmark. A very coarse finite element discretization for the membrane is used while keeping the fluid discretization unchanged. Although the discretization of the membrane is very coarse, the application of ECL guarantees a smooth representation of the FSI interface and the artificial kinks arising from the structural deformation on a coarse low-order mesh are not transferred into the fluid domain. Then, the FSI case of the turbulent flow around an air-inflated flexible hemisphere is carried out. A smooth $C^1$-continuous isogeometric discretization of the hemisphere is used for the structural boundary value problem, whereas the flow is described by a wall-resolved LES. The numerical results are compared with the experimental data and the numerical results obtained with the low-order FEM discretization of the membrane. The results suggest an improvement of the overall FSI behavior predicted by FSI-LES-IGA. This can be attributed to the smoothness and the better representation of the FSI interface by NURBS. Lastly, a significantly scaled deformation pattern of the hemisphere taken from one of the structural eigenmode shapes is imposed on the fluid FSI interface. Three gradually refined finite element meshes and the smooth isogeometric parametrization of the hemisphere used within the aforementioned FSI case are considered. A pure CFD analysis is performed for these four cases. The low-order FEM discretization combined with the large displacements imply a faceted representation of the surface, which is generating artificial surface roughness. The comparison of the time-averaged first- and second-order moments suggest a significant influence of this artificial interface roughness on the flow results. It highlights the necessity of interface smoothness for CFD and subsequently FSI simulations, showing the advantage of the presented isogeometric mortar-based mapping method as well as the isogeometric structural analysis in case of large deformation FSI simulations.

Although the proposed methodology does not ensure the conservation of the interface energy by construction, the results obtained for the FSI cases tested show that the corresponding error is negligible when the mortar-based mapping method is employed. Moreover, a significant improvement of the interface energy conservation is observed when using the isogeometric mortar-based mapping method in comparison to the standard mortar-based mapping method.

A smooth representation of the FSI interface obtained during simulations based on high-order discretizations globally leads to higher computational costs at each time step than low-order discretizations due to the higher computational
<figure><img src="image_38.png" /><figcaption>Fig. A.28. Steady-state structural analysis: Load-displacement curves for the numerical models against the experimental values.</figcaption></figure>

effort for evaluating the basis functions. Nevertheless, they are advantageous due to the following reasons: Firstly, the NURBS-based high-order parametrized geometries can be directly imported from CAD avoiding the costly mesh generation step, which is necessary for low-order FEM discretizations. Secondly, a smooth FSI interface leads to a faster convergence of the fluid solver reducing the required CPU-time. The best compromise between low computational costs and a smooth representation of the FSI interface is found with the proposed *Exact Coupling Layer* (ECL) concept.

Future work can be focused on extending the herein proposed isogeometric mortar-based method to complex real-world CAD geometries involving trimmed multipatches, thus enabling FSI on geometric models stemming directly from industrial scale CAD. Additionally, FSI problems with very large deformations can be studied in order to further investigate the influence of the FSI interface description.

## Acknowledgments

The authors thank F. Péan (Computer Vision Laboratory, Swiss Federal Institute of Technology in Zurich, Switzerland) and A. Ghantasala (Chair of Structural Analysis, Technical University of Munich, Germany) for their contributions to this work. The project is financially supported by the German Research Foundation (Deutsche Forschungsgemeinschaft) under the contract numbers BR 1847/12-2 and BL 306/26-2. The computations were partially carried out on the German Federal Top-Level Computer SuperMUC at LRZ Munich under the project number pr53ne.

# Appendix. Evaluation of isogeometric model for the flexible hemisphere

## A.1. Steady-state structural analysis

In order to evaluate the isogeometric model of the flexible hemisphere case, the steady calibration case described in Wood et al. (2018) is considered for the three isogeometric discretizations: One using 144 piecewise biquadratic $C^1$-continuous elements, one using 625 piecewise biquadratic (low-order) $C^1$-continuous elements and one using 625 piecewise bicubic (high-order) $C^2$-continuous elements. The corresponding geometries are obtained using the URS form-finding method as explained in Section 4.1.1. The steady calibration case works as follows: The pressure difference between the interior and the exterior of the membrane is steadily increased, while measuring the membrane displacement at the apex. The flexible structure is subjected to its self-weight. Its homogeneous and isotropic pre-stress is chosen equal to $n_{\text{membrane}} = 7794.5$ Pa as computed in De Nayer et al. (2018a).

The results of the three different isogeometric discretizations are shown in Fig. A.28 and compared with the experimental data of Wood et al. (2018) and the numerical FEM results of De Nayer et al. (2018a). It can be observed that all isogeometric discretizations deliver satisfactory results, whereas the apex displacement provided by the discretizations involving 625 biquadratic and bicubic elements are nearly identical. Therefore, it can be deduced that a sufficiently fine isogeometric discretization for obtaining a highly accurate discretization-independent solution is achieved for the steady validation case.
<figure><img src="image_39.png" /><figcaption>Fig. A.29. Dynamic structural analysis: Comparison of the simulation results with the experimental data.</figcaption></figure>

## A.2. Dynamic structural analysis

The three isogeometric discretizations used in A.1 are herein employed within the unsteady validation case introduced in Wood et al. (2018). An impulse load $F(t)$ of the form

$$F(t) = \begin{cases} 0.018 \text{ N}, & t < 0.001 \text{ s}, \\ 0, & \text{else}, \end{cases} \qquad (A.1)$$

is applied at an angle $\phi = \pi/4$ in the meridian direction (see De Nayer et al. (2018a) for more details). The standard second-order Newmark scheme with a time step of magnitude $\Delta t = 10^{-3}$ s is chosen for the time discretization as explained in Section 2.1.3. The simple Rayleigh damping model is applied with the mass- and the stiffness-proportional parameters set to $\alpha_t = 17.47$ and $\beta_t = 1.89 \times 10^{-4}$, respectively.

For the three isogeometric parametrizations the vertical displacements $d_z$ depicted in Fig. A.29(a) are compared with the results obtained by the finite element mesh with 1926 elements (De Nayer et al., 2018a) and the experimental data (Wood et al., 2018). In addition, the corresponding spectra are shown in Fig. A.29(b). It can be observed that a very good agreement of both isogeometric discretizations using 625 elements with the finite element solution and the experimental data is achieved, whereas the isogeometric discretization using 144 elements deviates significantly. Thus, the latter is not considered for the FSI simulation. Since both isogeometric discretizations with 625 elements deliver almost the same solution, the one with biquadratic elements is chosen for the FSI setup for efficiency reasons.

## References

* Aitken, A.C., 1926. On Bernoulli's numerical solution of algebraic equations. Proc. Roy. Soc. Edinburgh 46, 289–305.

* Apostolatos, Andreas, 2019. Isogeometric Analysis of Thin-Walled Structures on Multipatch Surfaces in Fluid-Structure Interaction. Chair of Structural Analysis, Technical University of Munich.

* Apostolatos, Andreas, Bletzinger, Kai-Uwe, Wüchner, Roland, 2018. Nitsche's method for form-finding of multipatch isogeometric membrane analysis. PAMM 18 (1), e201800106. http://dx.doi.org/10.1002/pamm.201800106, https://onlinelibrary.wiley.com/doi/abs/10.1002/pamm.201800106.

* Apostolatos, Andreas, Bletzinger, Kai-Uwe, Wüchner, Roland, 2019. Weak imposition of constraints for structural membranes in transient geometrically nonlinear isogeometric analysis on multipatch surfaces. Comput. Methods Appl. Mech. Engrg. http://dx.doi.org/10.1016/j.cma.2019.01.023, http://www.sciencedirect.com/science/article/pii/S0045782519300489.

* Apostolatos, A., Breitenberger, M., Wüchner, R., Bletzinger, K.-U., 2015. Domain decomposition methods and Kirchhoff-Love shell multipatch coupling in isogeometric analysis. In: Isogeometric Analysis and Applications.

* Apostolatos, A., Schmidt, R., Wüchner, R., Bletzinger, K.-U., 2014. A Nitsche-type formulation and comparison of the most common domain decomposition methods in isogeometric analysis. Internat. J. Numer. Methods Engrg. 97 (7), 473–504.

* Babuska, I., 1973. The finite element method with penalty. Math. Comp. 27 (122), 221–228.

* Basar, Y., Krätzig, W., 1985. Mechanik der Flächentragwerke. Vieweg, Braunschweig.

* Bazilevs, Y., Calo, V., Cottrell, J., Hughes, T., Reali, A., Scovazzi, G., 2007. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Comput. Methods Appl. Mech. Engrg. 197 (1), 173–201.

* Bazilevs, Y., Hsu, M.-C., Scott, M., 2012. Isogeometric fluid-structure interaction analysis with emphasis on non-matching discretizations, and with application to wind turbines. Comput. Methods Appl. Mech. Engrg. 249–252, 28–41.

* Bernardi, C., Maday, Y., Patera, A.T., 1994. A new nonconforming approach to domain decomposition: The mortar element method. In: Brezis, H., Lions, J.L. (Eds.), Nonlinear Partial Differential Equations and their Applications. Pitman Wiley, London/New York.

* Bletzinger, K.-U., Ramm, E., 1999. A general finite element approach to the form finding of tensile structures by the updated reference strategy. Int. J. Space Struct. 14 (2), 131–145.

* Bletzinger, K.-U., Wüchner, R., Daoud, F., Camprubí, N., 2005. Computational methods for form finding and optimization of shells and membranes. Comput. Methods Appl. Mech. Engrg. 194 (30), 3438–3452.

* Breitenberger, M., Apostolatos, A., Philipp, B., Wüchner, R., Bletzinger, K.-U., 2015. Analysis in computer aided design: Nonlinear isogeometric b-rep analysis of shell structures. Comput. Methods Appl. Mech. Engrg. 284, 401–457.
Breuer, M., 2018. Effect of inflow turbulence on an airfoil flow with laminar separation bubble: An LES study. Flow, Turbul. Combust. 101 (2), 433–456.

Breuer, M., De Nayer, G., Münsch, M., Gallinger, T., Wüchner, R., 2012. Fluid-structure interaction using a partitioned semi-implicit predictor-corrector coupling scheme for the application of large-eddy simulation. J. Fluids and Structures 29, 107–130.

Causin, P., Gerbeau, J.-F., Nobile, F., 2005. Added-mass effect in the design of partitioned algorithms for fluid-structure problems. Comput. Methods Appl. Mech. Engrg. 194, 4506–4527.

Clough, R.W., Penzien, J., 1993. Dynamics of Structures. McGraw-Hill, New York.

Codina, R., 2001. A stabilized finite element method for generalized stationary incompressible flows. Comput. Methods Appl. Mech. Engrg. 190 (20–21), 2681–2706.

Cottrell, J., Hughes, T., Bazilevs, Y., 2009. Isogeometric Analysis: Toward Integration of CAD and FEA. John Wiley & Sons Ltd, Chichester, United Kingdom.

de Boer, A., van Zuijlen, A., Bijl, H., 2007. Review of coupling methods for non-matching meshes. Comput. Methods Appl. Mech. Engrg. 196 (8), 1515–1525.

de Boer, A., van Zuijlen, A., Bijl, H., 2008. Comparison of conservative and consistent approaches for the coupling of non-matching meshes. Comput. Methods Appl. Mech. Engrg. 197 (49), 4284–4297.

De Lorenzis, L., Wriggers, P., Zavarise, G., 2012. A mortar formulation for 3D large deformation contact using NURBS-based isogeometric analysis and the augmented Lagrangian method. Comput. Mech. 49, 1–20.

De Nayer, G., Apostolatos, A., Wood, J.N., Bletzinger, K.U., Wüchner, R., Breuer, M., 2018a. Numerical studies on the instantaneous fluid-structure interaction of an air-inflated flexible membrane in turbulent flow. J. Fluids and Structures 82, 577–609. http://dx.doi.org/10.1016/j.jfluidstructs. 2018.08.005.

De Nayer, G., Breuer, M., 2014. Numerical FSI investigation based on LES: Flow past a cylinder with a flexible splitter plate involving large deformations (FSI-Pfs-2a). Int. J. Heat Fluid Flow 50, 300–315.

De Nayer, G., Kalmbach, A., Breuer, M., Sicklinger, S., Wüchner, R., 2014. Flow past a cylinder with a flexible splitter plate: A complementary experimental-numerical investigation and a new FSI test case (FSI-Pfs-1a). Comput. & Fluids 99, 18–43.

De Nayer, G., Schmidt, S., Wood, J.N., Breuer, M., 2018b. Enhanced injection method for synthetically generated turbulence within the flow domain of eddy-resolving simulations. Comput. Math. Appl. 75 (7), 2338–2355.

De Nayer, Guillaume, Wood, Jens Nikolas, Breuer, Michael, Apostolatos, Andreas, Wüchner, Roland, 2018c. Coupled simulations involving light-weight structures within turbulent flows: a complementary experimental and numerical application. PAMM 18 (1), e201800030. http://dx.doi.org/10. 1002/pamm.201800030, https://onlinelibrary.wiley.com/doi/abs/10.1002/pamm.201800030.

Dervieux, A., Farhat, C., Koobus, B., Vázquez, M., 2010. Total energy conservation in ALE schemes for compressible flows. Eur. J. Computat. Mech./Rev. Eur. Méc. Numér. 19 (4), 337–363.

Farhat, C., van der Zee, K.G., Geuzaine, P., 2006. Provably second-order time-accurate loosely-coupled solution algorithms for transient non-linear computational aeroelasticity. Comput. Methods Appl. Mech. Engrg. 195 (17–18), 1973–2001.

Felippa, C.A., Park, K.C., 1980. Staggered transient analysis procedures for coupled mechanical systems: Formulation. Comput. Methods Appl. Mech. Engrg. 24 (1), 61–111.

Gerbeau, J.-F., Vidrascu, M., 2003. A quasi-Newton algorithm based on a reduced model for fluid-structure interaction problems in blood flows. ESAIM Math. Model. Numer. Anal. 37 (4), 631–647.

Glück, M., Breuer, M., Durst, F., Halfmann, A., Rank, E., 2003. Computation of wind-induced vibrations of flexible shells and membranous structures. J. Fluids and Structures 17 (5), 739–765.

Holzapfel, G.A., 2000. Nonlinear Solid Mechanics: A Continuum Approach for Engineering. John Wiley & Sons Ltd., Chichester, United Kingdom.

Hughes, T., Cottrell, J., Bazilevs, Y., 2005. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement. Comput. Methods Appl. Mech. Engrg. 194, 4135–4195.

Kassiotis, C., Ibrahimbegovic, A., Niekamp, R., Matthies, H.G., 2011. Nonlinear fluid-structure interaction problem. Part I: Implicit partitioned algorithm, nonlinear stability proof and validation examples. Comput. Mech. 47 (3), 305–323.

Kiendl, J.M., 2011. Isogeometric Analysis and Shape Optimal Design of Shell Structures (Ph.D. thesis). Lehrstuhl für Statik, Technische Universität München, Germany.

Klein, M., Sadiki, A., Janicka, J., 2003. A digital filter based generation of inflow data for spatially-developing direct numerical or large-eddy simulations. J. Comput. Phys. 186, 652–665.

Klöppel, T., Popp, A., Küttler, U., Wall, W.A., 2011. Fluid-structure interaction for non-conforming interfaces based on a dual mortar formulation. Comput. Methods Appl. Mech. Engrg. 200 (45), 3111–3126.

Mok, D.P., 2001. Partitionierte Lösungsansätze in der Strukturdynamik und der Fluid-Struktur-Interaktion (Ph.D. thesis). Institut für Baustatik, Universität Stuttgart, Germany.

Newmark, N.M., 1959. A method of computation for structural dynamics. J. Eng. Mech. Div. 85 (3), 67–94.

Park, K.C., Felippa, C.A., Rebel, G., 2002. A simple algorithm for localized construction of non-matching structural interfaces. Internat. J. Numer. Methods Engrg. 53 (9), 2117–2142.

Philipp, B., Breitenberger, M., D'Auria, I., Wüchner, R., Bletzinger, K.-U., 2016. Integrated design and analysis of structural membranes using the isogeometric b-rep analysis. Comput. Methods Appl. Mech. Engrg. 303, 312–340.

Piegl, L., Tiller, W., 2012. The NURBS Book. Springer Science & Business Media, Berlin Heidelberg.

Piperno, S., Farhat, C., 2001. Partitioned procedures for the transient solution of coupled aeroelastic problems–Part II: Energy transfer analysis and three-dimensional applications. Comput. Methods Appl. Mech. Engrg. 190 (24–25), 3147–3170.

Schmidt, S., Breuer, M., 2014. Hybrid LES-URANS methodology for the prediction of non-equilibrium wall-bounded internal and external flows. Comput. & Fluids 96, 226–252.

Schmidt, S., Breuer, M., 2017. Source term based synthetic turbulence inflow generator for eddy-resolving predictions of an airfoil flow including a laminar separation bubble. Comput. & Fluids 146, 1–22.

Schweizerhof, K., Ramm, E., 1984. Displacement dependent pressure loads in nonlinear finite element analyses. Comput. Struct. 18 (6), 1099–1114.

Sen, S., De Nayer, G., Breuer, M., 2017. A fast and robust hybrid method for block-structured mesh deformation with emphasis on FSI-LES applications. Internat. J. Numer. Methods Engrg. 111 (3), 273–300.

Sicklinger, S., 2014. Stabilized Co-Simulation of Coupled Problems Including Fields and Signals (Ph.D. thesis). Lehrstuhl Für Statik, Technische Universität München, Germany.

Sicklinger, S., Belsky, V., Engelmann, B., Elmqvist, H., Olsson, H., Wüchner, R., Bletzinger, K.-U., 2014. Interface Jacobian-based co-simulation. Internat. J. Numer. Methods Engrg. 98 (6), 418–444.

Smagorinsky, J., 1963. General circulation experiments with the primitive equations I: The basic experiment. Mon. Weather Rev. 91 (3), 99–165.

Teschemacher, T., Bauer, A.M., Oberbichler, T., Breitenberger, M., Rossi, R., Wüchner, R., Bletzinger, K.-U., 2018. Realization of cad-integrated shell simulation based on isogeometric b-rep analysis. Adv. Model. Simul. Eng. Sci. 5 (1), 276.

Thompson, J.F., Warsi, Z.U.A., Wayne Mastin, C., 1985. Numerical Grid Generation: Foundations and Applications, vol. 45. North-Holland.
Thomsen, K., Petersen, J.T., Nim, E., Øye, S., Petersen, B., 2000. A method for determination of damping for edgewise blade vibrations. Wind Energy 3 (4), 233–246.

Wall, W.A., 1999. Fluid-Struktur-Interaktion Mit Stabilisierten Finiten Elementen (Ph.D. thesis). Institut für Baustatik, Universität Stuttgart, Germany.

Wang, T., 2016. Development of Co-Simulation Environment and Mapping Algorithms (Ph.D. thesis). Lehrstuhl für Statik, Technische Universität München, Germany.

Wang, T., Wüchner, R., Sicklinger, S., Bletzinger, K.-U., 2016. Assessment and improvement of mapping algorithms for non-matching meshes and geometries in computational FSI. Comput. Mech. 57 (5), 793–816.

Wood, J.N., Breuer, M., De Nayer, G., 2018. Experimental studies on the instantaneous fluid-structure interaction of an air-inflated flexible membrane in turbulent flow. J. Fluids and Structures 80, 405–440. http://dx.doi.org/10.1016/j.jfluidstructs.2018.02.006.

Wood, J.N., De Nayer, G., Schmidt, S., Breuer, M., 2016. Experimental investigation and large-eddy simulation of the turbulent flow past a smooth and rigid hemisphere. Flow, Turbul. Combust. 97 (1), 79–119.

Wüchner, Roland, Apostolatos, Andreas, De Nayer, Guillaume, Breuer, Michael, Bletzinger, Kai-Uwe, 2018. Coupled simulations involving light-weight structures within turbulent flows: FSI strategy and non-matching interface treatment for isogeometric b-rep analysis. PAMM 18 (1), e201800107. http://dx.doi.org/10.1002/pamm.201800107, https://onlinelibrary.wiley.com/doi/abs/10.1002/pamm.201800107.

Wüchner, R., Bletzinger, K.-U., 2005. Stress-adapted numerical form finding of pre-stressed surfaces by the updated reference strategy. Internat. J. Numer. Methods Engrg. 64 (2), 143–166.

Wüchner, R., Kupzok, A., Bletzinger, K.-U., 2007. A framework for stabilized partitioned analysis of thin membrane-wind interaction. Internat. J. Numer. Methods Engrg. 54 (6–8), 945–963.