Fluid-structure interaction with NURBS-based coupling

TL;DR

This paper presents a novel fluid-structure interaction method using a shared NURBS-based spline interface, enabling precise geometry communication and flexible refinement in coupled analysis.

Contribution

It introduces a coupled approach employing NURBS for both structural and fluid domains, facilitating direct transfer and independent refinement of interface quantities.

Findings

Enables exact geometry representation at the interface

Allows independent refinement of fluid and structural models

Facilitates direct transfer of interface quantities

Abstract

Engineering design via CAD software relies on Non-Uniform Rational B-Splines (NURBS) as a means for representing and communicating geometry. Therefore, in general, a NURBS description of a given design can be considered the exact description. The development of isogeometric methods has made the geometry available to analysis methods Hughes et al. (2005). Isogeometric analysis has been particularly successful in structural analysis; one reason being the wide-spread use of two-dimensional finite elements in this field. For fluid dynamics, where three-dimensional analysis is usually indispensable, isogeometric methods are more complicated, yet of course not impossible, to apply in a general fashion. This paper describes a method that enables the solution of fluid-structure interaction with a matching spline description of the interface. On the structural side, the spline is used in an isogeometric setting. On the fluid side, the same spline is used in the framework of a NURBS-enhanced finite element method (extension of Sevilla et al. (2011)). The coupling of the structural and the fluid solution is greatly facilitated by the common spline interface. The use of the identical spline representation for both sides permits a direct transfer of the necessary quantities, all the while still allowing an adjusted, individual refinement level for both sides.

Tables

Table 1. Table 1: Rigid featflow benchmark parameters.

parameter	identifier	value
inflow velocity	U	0.2*1.5 $m/s$
dynamic viscosity	$\mu$	0.001 $kg/ms$
density	$\rho$	1.0 $kg/m^3$
cylinder diameter	D	0.1 $m$
Reynolds number	Re	20

Table 2. Table 3: Flexible cylinder parameters.

parameter	identifier	value
inflow velocity	U	0.10282 $m/s$
dynamic viscosity	$\mu$	0.0000182 $kg/ms$
density (fluid)	${\rho }^f$	1.18 $kg/m^3$
cylinder diameter	D	0.006 $m$
Reynolds number	Re	40
Young’s modulus	$E$	100 Pa
Poisson number	${\nu }^s$	0.3
density (structure)	${\rho }^s$	1000 $kg/m^3$

Table 3. Table 5: Moving flag featflow benchmark parameters.

parameter	identifier	value
inflow velocity	U	1.0*1.5 $m/s$
dynamic viscosity	$\mu$	1.0 $kg/ms$
density	${\rho }^f$	1000.0 $kg/m^3$
cylinder diameter	D	0.1 $m$
Reynolds number	Re	100
Young’s modulus	$E$	1.4 10e6 $Pa$
Poisson number	${\nu }^s$	0.4
density (structure)	${\rho }^s$	10000.0 $kg/m^3$

Equations

$${\rho }^s\frac{d^2{\mathbf{d}}^s}{d{t}^2}=\nabla \cdot {\mathbf{\sigma }}^s+{\mathbf{b}}^s\text{mbox}on({\Omega }_t^s)\forall t\in [0,T],$$

$${\mathbf{S}}^s(\mathbf{E})={\lambda }^{\mathbf{s}}\mathbf{tr}\left({\mathbf{E}}^{\mathbf{s}}\right)+2{\mu }^{\mathbf{s}}{\mathbf{E}}^{\mathbf{s}}$$

$$\mathbf{E}=\frac{1}{2}\left({\mathbf{F}}^T\mathbf{F}-\mathbf{I}\right),$$

$$\mathbf{F}=\frac{\partial \mathbf{x}}{\partial {\mathbf{x}}_0}.$$

$${\lambda }^s=\frac{{\nu }^s{E}^s}{(1+{\nu }^s)(1-2{\nu }^s)},{\mu }^s=\frac{E^s}{2(1+{\nu }^s)}.$$

$${\rho }^s\frac{d^2{\mathbf{d}}^s}{d{t}^2}={\nabla }_0\cdot \left({\mathbf{S}}^s{\mathbf{F}}^T\right)+{\mathbf{b}}^s\text{mbox}on({\Omega }_0^s)\forall t\in [0,T].$$

$${\mathbf{d}}^s={\mathbf{g}}^s\mathrm{on}{\left({\Gamma }_t^s\right)}_g,{\mathbf{n}}^s\cdot {\mathbf{\sigma }}^s={\mathbf{h}}^s\mathrm{on}{\left({\Gamma }_t^s\right)}_h.$$

$${\mathbf{d}}^s\left(\mathbf{x},t=0\right)={\mathbf{d}}^0.$$

$${\rho }^f\left(\frac{\partial {\mathbf{u}}^f}{\partial t}+{\mathbf{u}}^f\cdot \nabla {\mathbf{u}}^f-{\mathbf{f}}^f\right)-\nabla \cdot {\mathbf{\sigma }}^f$$

$${\nabla \cdot \mathbf{u}}^f$$

$${\mathbf{\sigma }}^f({\mathbf{u}}^f,p^f)$$

$${\mathbf{\epsilon }}^f({\mathbf{u}}^f)$$

$${\mathbf{u}}^f$$

$${\mathbf{n}}^f\cdot {\mathbf{\sigma }}^f$$

$${\mathbf{u}}^f(\mathbf{x},0)={\mathbf{u}}^0(\mathbf{x})\text{mbox}in{\Omega }_t^f\text{mbox}att=0.$$

$$\begin{array}{rl}{\mathbf{d}}^f\left(\mathbf{x},t\right)& ={\mathbf{d}}^s\left(\mathbf{x},t\right)\mathrm{on}{\Gamma }_{FS},\\ {\mathbf{u}}^f\left(\mathbf{x},t\right)& ={\mathbf{u}}^s\left(\mathbf{x},t\right)\mathrm{on}{\Gamma }_{FS}.\end{array}$$

$${\mathbf{\sigma }}^f\left(\mathbf{x},t\right)\cdot {\mathbf{n}}^f=-{\mathbf{\sigma }}^s\left(\mathbf{x},t\right)\cdot {\mathbf{n}}^s\mathrm{on}{\Gamma }_{FS}.$$

$${\tilde{\mathcal{S}}}_d^{s,h}$$

$${\tilde{\mathcal{V}}}_d^h$$

$$\delta W\left({\mathbf{w}}^{s,h},{\mathbf{d}}^{s,h}\right)={\int }_{{\Omega }^s}{\mathbf{w}}^{s,h}\cdot \frac{d^2{\mathbf{d}}^{s,h}}{d{t}^2}+\nabla {\mathbf{w}}^{s,h}:{\mathbf{\sigma }}^{s}d\Omega -{\int }_{{\Gamma }_g^s}{\mathbf{w}}^{s,h}\cdot {\mathbf{h}}^{s,h}d\Gamma =0$$

$$\mathbf{C}(\Theta )=\sum _{i=1}^n{R}_{i,p}(\Theta ){\mathbf{P}}_i,$$

$${\mathbf{d}}^{s,h}(\Theta )=\sum _{i=1}^n{R}_{i,p}(\Theta ){\mathbf{d}}_i^s.$$

$$({\mathcal{S}}_{\mathbf{u}}^h{)}_n$$

$$({\mathcal{V}}_{\mathbf{u}}^h{)}_n$$

$$({\mathcal{S}}_p^h{)}_n$$

$$\displaystyle\begin{split}\int_{Q_{n}}{\bf w}^{f,h}\cdot\rho^{f}\left(\frac{\partial{\bf u}^{f,h}}{\partial t}+{{\bf u}^{f,h}\cdot{\bf\nabla u}^{f,h}}-{\bf f}^{f,h}\right)\;dQ+\int_{Q_{n}}\mbox{\boldmath$\varepsilon$}^{f}({\bf w}^{f,h}):\mbox{\boldmath$\sigma$}^{f,h}(p^{f,h},{\bf u}^{f,h})\;dQ\\ +\int_{Q_{n}}q^{f,h}\nabla\cdot{\bf u}^{f,h}\;dQ+\int_{\Omega_{n}}({\bf w}^{f,h})^{+}_{n}\cdot\rho^{f}\left(({\bf u}^{f,h})^{+}_{n}-({\bf u}^{f,h})^{-}_{n}\right)\;d\Omega\\ +\sum_{e=1}^{({n_{el})}_{n}}\int_{Q^{e}_{n}}\tau_{\mbox{\tiny{MOM}}}\frac{1}{\rho^{f}}\left[\rho^{f}\left(\frac{\partial{\bf w}^{f,h}}{\partial t}+{{\bf u}^{f,h}\cdot{\bf\nabla w}^{f,h}}\right)-\nabla\cdot\mbox{\boldmath$\sigma$}^{f,h}(q^{f,h},{\bf w}^{f,h})\right]\\ \cdot\left[\rho^{f}\left(\frac{\partial{\bf u}^{f,h}}{\partial t}+{{\bf u}^{f,h}\cdot{\bf\nabla u}^{f,h}}-{\bf f}^{f,h}\right)-\nabla\cdot{\bm{\sigma}}^{f,h}(p^{f,h},{\bf u}^{f,h})\right]\;dQ+\sum_{e=1}^{({n_{el})}_{n}}\int_{Q^{e}_{n}}\tau_{\mbox{\tiny{CONT}}}\nabla\cdot{\bf w}^{f,h}\rho^{f}\nabla\cdot{\bf u}^{f,h}\;dQ\\ =\int_{P_{n}}{\bf w}^{f,h}\cdot{\bf h}^{f}\;dP.\end{split}$$

$$({\mathbf{u}}^h{)}_n^{\pm }$$

$${\int }_{Q_n}\dots dQ$$

$${\int }_{P_n}\dots dP$$

$$({\mathbf{u}}^h{)}_0^{+}={\mathbf{u}}_0.$$

$$\mathbf{\Phi }(s,r)=(1-s-r){\mathbf{x}}_2+(s+r)\mathbf{C}\left(\frac{{\Theta }_{1}s+{\Theta }_{3}r}{s+r}\right)$$

$$\mathbf{\sigma }\left(r\right)=\sum _i{L}_i\left(r\right){\mathbf{\sigma }}_i^f,$$

$${\mathbf{F}}_i^f={\int }_{{\Gamma }_{FS}}{L}_i\left(r\right)\sum _j\left[L_j\left(r\right){\mathbf{\sigma }}_j^f\cdot \mathbf{n}\right]dr,$$

$${\mathbf{F}}_i^s=\sum _j^{n{n}_{{\Gamma }_{FS}}^f}{R}_i\left({\Theta }_j\right){\mathbf{F}}_j^f.$$

Full text

††journal: Computer Methods in Applied Mechanics and Engineering

1 Introduction

Geometries for engineering applications are generated using Computer-Aided-Design (CAD) systems: the CAD model is what we assume to be the exact geometry. Nowadays, all major CAD systems share one common basis for geometry representation: Non-Uniform Rational B-Splines (NURBS). NURBS provide a common standard to exchange geometry information. In the classic discretization methods, the exact NURBS geometry is lost when a finite element mesh — which is in general only an approximation of the geometry — is generated. Spline-based methods, such as Isogeometric Analysis [1, 3] or NURBS-Enhanced Finite Elements [2, 4] make use of the NURBS format in order to integrate the exact geometry into the finite element method, thus avoiding the approximative character of the finite element mesh. Note, the change of computational domain due to numerical analysis with isogeometric analysis, e.g., for structural deformation, is only an approximation, even though the NURBS representation of the initial geometry can be considered as exact. This paper presents a fluid-structure-interaction (FSI) method with spline-based, smooth geometry representation of the interface between fluid and structure.

To give a general overview, we will categorize different approaches to FSI according to (1) the solution process, (2) the temporal coupling, and (3) the spatial coupling.

(1) With respect to the solution method, one differentiates between monolithic approaches (e.g., [5, 6, 7]) or partitioned approaches (e.g., [8, 9, 10, 11, 12]). The monolithic approach includes the full solution process in one single solver: The main advantage gained is the robustness (in the sense of stability and the capability for large time steps); the price to pay is the significant implementation effort. Partitioned approaches rely on individual solvers for all physical systems, which generate individual contributions to the coupled system. The necessary communication of solution data is restricted to the interfaces between the single field solvers. This facilitates the use of existing, already very evolved solvers, but leads to stability issues. Only in the case of partitioned solvers, there is a need for temporal and spatial coupling.

(2) Methods for temporal coupling can be categorized into strong and weak coupling schemes. When using a strong coupling method, the idea is to iterate between the individual solvers until a converged solution for the current time step is obtained [11]. Within a weak coupling, each individual solver is only called once per time step; this usually requires the use of predictor-corrector methods [9, 13]. Weak coupling schemes profit from reduced computational cost. However, they may suffer from a lack of stability in cases where the density ratio between fluid and solid is close to one [14]. This limit may very well be reached for incompressible flows in combination with flexible structures.

(3) With regard to the spatial coupling, the use of individual solvers leads to – in general – non-conforming discretizations for fluid and structure. Consequently, the transfer of interfacial quantities (e.g., forces, displacements) from one computational mesh to the other requires projection methods [15]. Again, there are two fundamental approaches available: mesh-based and mesh-free methods. Mesh-based methods, such as finite-interpolation-elements (FIE), rely on the direct approaches of the individual solvers [16, 17]. An example for a mesh-free method is the spline-based moving-least-squares-method (MLS) [18], which interpolates the transfer information between the different interface discretizations. Both methods can also be adapted to individual solvers with spline-based solution methods [19, 20, 21].

In this work, we use a partitioned approach with strong coupling and rather specialized solvers. The key idea of our FSI approach is to generate a smooth and conforming spline description of the interface between fluid and structure, that might be also benefical for monolithic schemes. On the fluid side, the spline is supported through the NURBS-Enhanced Finite Element Method. On the structural side, this is done by means of IGA. The use of the identical spline representation for both sides permits a direct transfer of the necessary quantities, all the while still allowing an adjusted, individual refinement level for both sides. Although there are several approaches that use IGA on the structural side, to our knowledge this is the first time where a conforming interface for both sides was established using splines.

The paper has five main sections. Section 2 is devoted to the underlying governing equations. Sections 3–5 outline the numerical methods for fluid, structure, and coupling. Numerical examples are discussed in Section 6.

2 Modeling: Governing Equations

The test cases in this work feature a nonlinear elastic structure enclosed in an incompressible Newtonian fluid. This section introduces the relevant governing equations for the structure and the fluid.

2.1 Structural Deformation

The response of a structure to an external load is governed by classical elastodynamics. The primal variable is the displacement field, here denoted as ${\bf{d}}^{s}\left({\bf{x}},t\right)$. It describes the change between the current structural configuration ${\bf{x}}$ and the initial configuration ${\bf{x}}_{0}$. These alterations can be either due to deformations of the structure or due to rigid body motions.

Consider a deformable computational domain, which is, at each instant in time, denoted by $\Omega_{t}^{s}$. $\Omega_{t}^{s}$ is a subset of $\mathbb{R}^{nsd^{s}}$, with $nsd^{s}$ as the number of space dimensions for the structure. Its boundary is denoted as $\Gamma_{t}^{s}$. At each point in time $t\in[0,T]$, the displacement ${\bf{d}}\left({\bf{x}},t\right)$ of the structure is governed by an equation based on Newton’s second law:

[TABLE]

where $\rho^{s}$ indicates the density of the structure, ${\bf{b}}^{s}$ the prescribed body force per unit volume and $\bm{\sigma^{s}}$ the Cauchy stress tensor.

As constitutive relation, we employ a geometrically nonlinear approach for a hyperelastic material model; a model which enables us to account for large deformations. More specifically, we employ the St. Venant Kirchhoff model, which describes the constitutive equation for the stress tensor based on a stress-strain relation using the Green-Lagrange strain tensor ${\bf{E}}^{s}$ and the 2nd Piola-Kirchhoff stress tensor $\bf{S}^{s}(\bf{E})$ [22]. They are defined as:

[TABLE]

and

[TABLE]

where the deformation gradient ${\bf{F}}$ is defined as

[TABLE]

Here, $\lambda^{s}$ and $\mu^{s}$ are the Lamé parameters. They can be expressed in terms of the Young’s modulus $E^{s}$ and the Poisson ratio $\nu^{s}$ as

[TABLE]

Connecting all of this information, we can express the equation of motion in the reference configuration as

[TABLE]

In order to obtain a well-posed system, boundary conditions have to be imposed on $\Gamma^{s}$. Here, we distinguish between Dirichlet $({\bf{g}}^{s})$ and Neumann $({\bf{h}}^{s})$ boundary conditions given by:

[TABLE]

$\left(\Gamma_{t}^{s}\right)_{g}$ and $\left(\Gamma_{t}^{s}\right)_{h}$ denote the Dirichlet and Neumann part of the boundary, forming a complementary subset of $\Gamma^{s}_{t}$, i.e., $\left(\Gamma_{t}^{s}\right)_{g}\cup\left(\Gamma_{s}^{f}\right)_{h}=\Gamma^{s}_{t}$ and $\left(\Gamma_{t}^{s}\right)_{g}\cap\left(\Gamma_{t}^{s}\right)_{h}=\emptyset$. Here, ${\bf n}^{s}$ refers to the outer normal vector on $\Gamma_{t}^{s}$.

Furthermore, the initial displacement is prescribed as initial condition:

[TABLE]

2.2 Fluid Flow

Consider a deformable fluid domain, which is, at each instant in time, denoted by $\Omega_{t}^{f}$. Domain $\Omega_{t}^{f}$ is a subset of $\mathbb{R}^{nsd^{f}}$, with $nsd^{f}$ as the number of space dimensions of the fluid domain. Then at each point in time $t\in[0,T]$, the velocity, ${\bf u}^{f}({\bf x},t)$, and the pressure, $p^{f}({\bf x},t)$, of the fluid are governed by the unsteady, incompressible Navier-Stokes equations:

[TABLE]

with $\rho^{f}$ as the fluid density. In the Newtonian case, the stress tensor $\bm{\sigma}^{f}$ is defined as

[TABLE]

with

[TABLE]

where $\mu^{f}$ denotes the dynamic viscosity. ${\bf f}^{f}$ includes all external body forces per unit mass of fluid.

In order to obtain a well-posed system, boundary conditions have to be imposed on the external boundary of $\Omega^{f}_{t}$, denoted as $\Gamma^{f}_{t}$. Here, we distinguish between Dirichlet and Neumann boundary conditions given by:

[TABLE]

where ${\bf g}^{f}$ and ${\bf h}^{f}$ are prescribed velocity and stress values. $\left(\Gamma_{t}^{f}\right)_{g}$ and $\left(\Gamma_{t}^{f}\right)_{h}$ denote the Dirichlet and Neumann part of the boundary, forming a complementary subset of $\Gamma^{f}_{t}$, i.e., $\left(\Gamma_{t}^{f}\right)_{g}\cup\left(\Gamma_{t}^{f}\right)_{h}=\Gamma^{f}_{t}$ and $\left(\Gamma_{t}^{f}\right)_{g}\cap\left(\Gamma_{t}^{f}\right)_{h}=\emptyset$. Here, ${\bf n}^{f}$ refers to the outer normal vector on $\Gamma_{t}^{f}$.

In the transient case, a divergence-free velocity field for the whole computational domain is needed as an initial condition:

[TABLE]

2.3 Coupling Conditions at the Fluid-Structure Interface

The governing equations of fluid (Equations (8)– (9)) and structure (Equation (5)) need to be connected in order to represent the interaction between the two components. This interaction takes place only through the common interface $\Gamma_{FS}=\Gamma^{f}_{t}\cup\Gamma^{s}_{t}$, not the full volume; thus distinguishing FSI problems from other multiphysics problems [23].

For a consistent coupling, the following physical requirements are essential: (1) geometric compatibility between the fields, (2) kinematic and dynamic conditions at the shared interface $\Gamma_{FS}$, and (3) conservation of mass, momentum and energy. This leads to the following coupling conditions on $\Gamma_{FS}$:

Kinematic continuity:

[TABLE]

These coupling conditions ensure the continuity of displacements and velocities across the interface.

Dynamic continuity:

[TABLE]

In agreement with Newton’s third law – $Actio\;and\;Reactio$ – this coupling condition enforces continuity of fluid ($\bm{\sigma}^{f}$) and structural stresses ($\bm{\sigma}^{s}$) at the interface $\Gamma_{FS}$ [24].

3 Numerical Methods: Structural Solution

The governing equation for the structural deformation, Equation (5), is discretized using an isogeometric finite element method. The time discretization is performed using a generalized $\alpha$ scheme [25, 26]. The resulting weak form is linearized using a Newton-Raphson framework.

3.1 Variational Form

The employed isogeometric finite element method is a standard Galerkin formulation; only the interpolation functions differ from standard Lagrange interpolation functions. We introduce a finite-dimensional space $\mathcal{I}_{0}\subset C^{0}(\overline{\Omega^{s}})$ based on NURBS interpolation as spatial discretization. To obtain the weak form, Equation (5) is multiplied by the displacement test function ${\bf w}^{s}$, integrated, and the stress term $\bm{\sigma}^{s}$ is integrated by parts. The following finite element interpolation and weighting function spaces for the displacement ${\bf d}^{s}$ can be defined:

[TABLE]

Based on these spaces, the discrete formulation of the structure is defined as: Find ${\bf d}^{s,h}\in\tilde{\mathcal{S}}_{d}^{h}$ such that:

[TABLE]

holds for all ${\bf w}^{s,h}\in\tilde{\mathcal{V}}_{d}^{h}$.

3.2 Isogeometric Analysis (IGA)

Both the discrete function spaces (Equation (17)) and the weak form (Equation (19)) assume the use of isogeometric interpolation functions. In our implementation, we have resorted to a standard IGA method, as it was introduced in [1, 3]. In view of the ample publications in the field of structural analysis with isogeometric finite elements, e.g., [27], we will restrict ourselves to describing the general idea of IGA. IGA invokes an isoparametric finite element concept. This means that the unknown solution and the geometry are interpolated utilizing the same type of interpolation functions. In the case of IGA, these interpolation functions are inspired by the standard CAD description employed in design engineering. CAD systems are usually based on either NURBS [28] or, in very rare cases, also on T-splines [29]. Both cases belong to the category of parametric geometry description, with basis functions defined in terms of a local parameter and a mapping to the physical space; a concept very similar to Lagrange interpolation functions. The key advantages of employing the NURBS or T-spline basis are (1) exact geometry description in the initial state and (2) user-controlled smoothness of the basis.

A NURBS curve is represented through a combination of control points that guide the curve — ${\bf P}_{i}$ indicates the control point coordinates — and the NURBS basis functions $R_{i,p}$ of degree $p$. The curve definition is as follows:

[TABLE]

with $\Theta$ as the local parameter of the NURBS curve. Due to the isoparametric concept, the unknown solution ${\bf d}^{s,h}(\Theta)$ is represented in the following form:

[TABLE]

4 Numerical Methods: Flow Solution

The Navier-Stokes equations (8)–(9) are discretized using P1P1 finite elements, i.e., linear interpolation for both the velocity and pressure degrees of freedom [30, 31]. Along the boundary, the linear finite elements are supplemented by geometry information using the NURBS-Enhanced Finite Element Method described in Section 4.2 [4]. P1P1 finite elements are known to violate the Ladyzhenskaya-Babuska-Brezzi (LBB) compatibility condition. Consequently, without appropriate stabilization, the pressure field is likely to present spurious and oscillatory results. The stabilization technique used here is Galerkin/Least-Squares (GLS) stabilization. In the GLS method, the stabilization term consists of an element-by-element weighted least-squares form of the original differential equation [32].

The required deformability of the computational domain over time is included through the Deformable-Spatial-Domain/Stabilized Space-Time (DSD/SST) method [30]. In contrast to the classic choice within the finite element community of finite element discretization in the spatial domain and finite differences in the time direction, this method relies on finite elements in both space and time (cf. Figure 3). This approach allows to formulate the variational form directly over the deforming domain. The necessity of modifying the fluid-flow equations based on the mesh deformation velocity is avoided. As long as the domain deformation remains within a certain range, the mesh can be equipped with the necessary flexibility to account for the mesh motion without the need for remeshing. For this purpose, the Elastic Mesh Update Method [33] is employed.

4.1 Variational Form

In order to construct the finite element function spaces for the space-time method, the time interval $[0,T]$ is divided into subintervals $I^{f}_{n}=(t_{n},t_{n+1})$, with $t_{n}$ and $t_{n+1}$ representing an ordered series of time levels $0=t_{0}<t_{1}<\cdots<t_{N}=T$. Now, if $\Omega_{n}^{f}=\Omega_{t_{n}}^{f}$, the space-time slab $Q_{n}^{f}$ is defined as the domain enclosed by the surfaces $\Omega_{n}^{f},\Omega_{n+1}^{f}$ as well as the surface desribed by $\partial\Omega_{t}^{f}$ as t traverses $I_{n}^{f}$, which shall be named $P_{n}^{f}$. The following finite element interpolation and weighting function spaces for velocity ${\bf u}$ and pressure $p$ can be defined based on linear, $C^{0}$-continuous interpolation in space and also linear, but discontinuous interpolation in time:

[TABLE]

The stabilized space-time formulation of the incompressible Navier-Stokes equations (8)–(9) for deforming domains can then be expressed as follows: Given $({\bf u}^{h})^{-}_{n}$ find ${\bf u}^{h}\in(\mathcal{S}^{h}_{\bf u})_{n}$ and $p^{h}\in(\mathcal{S}^{h}_{p})_{n}$ such that $\forall{\bf w}^{h}\in(\mathcal{V}^{h}_{\bf w})_{n}$, $\forall q^{h}\in(\mathcal{V}^{h}_{p})_{n}$:

[TABLE]

In the above equation, the following notation is used:

[TABLE]

The problem is solved sequentially for each space-time slab, starting with:

[TABLE]

Details on the method and its parameters $\tau_{\mbox{\tiny{MOM}}}$ and $\tau_{\mbox{\tiny{CONT}}}$ can be found in [34].

Within the stabilization terms, $\nabla\cdot{\bm{\sigma}}^{f}$ needs to be computed. This involves the computation of second order derivatives of the velocity field, which are zero in case of linear shape functions. Thus, a least-squares recovery technique is applied to improve consistency of the method [35]. By that, the fluid stresses are computed, which are needed to evaluate tractions on the surface later on.

4.2 NURBS-Enhanced Finite Element Method

Comparable to the IGA approach pursued for the structure (cf. Section 3), also the fluid side benefits from an integration of the CAD geometry. As an additional aspect, the use of the CAD geometry also in the flow solver ensures an excellent compatibility in the fluid-structure-coupling. In principle, it would be possible to also solve the fluid equation (25) using IGA. However, despite recent progress in a variety of directions [36, 37], the generation of closed volume splines describing complex geometries remains a challenge. As an alternative on middle ground between isogeometric analysis and standard finite elements, Sevilla, Fernandez-Mendes, and Huerta have proposed the NURBS-Enhanced Finite Element Method (NEFEM) [38, 39, 2]. Already in [4], the NEFEM has been extended to space-time finite elements. The FSI-coupling requires further modifications detailed in Section 4.2.2

4.2.1 The Original Formulation of NEFEM

In NEFEM, instead of the full geometry, only its boundary is represented using NURBS at the cost of maximally a 2D spline (the boundary of a three-dimensional object). In the interior of the geometry, a standard finite element mesh is utilized, which preserves all advantages of existing meshing algorithms (cf. Figure 3). Note that this approach leads to two different kinds of elements: (1) standard finite elements in the interior and (2) elements with a NURBS edge alongside the boundary. Usually, elements of category (1) will be in the vast majority, keeping the computation very efficient. Through the elements of category (2), the geometry is made available during the process of evaluating the integrals of the appropriate weak form (e.g., of the variational equation (25)): The integration domain $\left(\Omega^{f}\right)^{h}$ is no longer only an approximation of the real domain $\Omega^{f}$. Translated to the numerical implementation, this factor is incorporated through the position of the integration points needed in the finite element method during numerical quadrature: These are determined from the curved NURBS geometry and not from the approximated geometry (cf. Figure 4 for a schematic illustration). If we consider boundary integrals — e.g., for the evaluation of the forces the fluid exerts on the structure — the integration points are distributed on the curved geometry.

The unknown function continues to be represented using Lagrange polynomials. These can be of any order, in particular independent of the order of the NURBS basis. In the NEFEM version proposed in [2, 38, 39] the choice was made to compute the shape functions in the global coordinates. This variant — termed by the authors as Cartesian FEM — has the advantage that one obtains undistorted polynomial functions in the global elements, an attribute of particular importance when used with polynomials of higher order. The price to pay is an oddity that arises from this representation: negative shape function values as well as values larger than one may occur along the boundary (cf. Figure 4).

In terms of implementation, there are two important modifications with respect to a standard finite element method: (1) the definition and evaluation of the shape functions and (2) the placement of the quadrature points. Regarding (1), Sevilla performs all evaluations in the global coordinate system. With respect to (2), the original method uses a bi-unit square reference element; even for triangular elements. This enables the use of tensorized 1D Gauss-Legendre quadrature rules. Furthermore, a clear distinction between the NURBS direction and the standard direction can be made; leading to straight interior edges.

4.2.2 The Modified Version of NEFEM Used in the FSI Context

The NEFEM implementation utilized in this work is based on the same idea as it was proposed by [38]. The differences lie in the fact that (1) we employ the space-time version, as derived in [4], and (2) we refrain from using Cartesian FEM, but compute both the shape functions and the position of the integration points on a — now triangular — reference element. What remains unchanged is the overall concept of adjusting of the position of the integration points to the curved NURBS shape as well as the representation of the unknown solution with Lagrange polynomials — even if restricted to linear polynomials in our case.

The reason for computing the shape functions on a reference element — and thereby accepting the complication of distorted polynomial functions (cf. Figure 6) — lies in the fact that for Cartesian FEM, also the shape functions of the interior node are non-zero along the boundary edge/face. As a consequence, interior nodes contribute to boundary integrals. In the context of Dirichlet boundaries and boundary integrals, this shape function definition, however, leads to a grave disadvantage: the values along the boundary edge are always affected by nodal values of interior nodes. Since both Dirichlet boundary conditions and the boundary integral connected to the load transfer in FSI are of utter importance to the scenarios considered in this work, it was imperative to modify the shape function definition such that shape functions of inner nodes are zero along the curved edge. The shape function definition — inspired by the so-called "p-FEM" in [38] — is performed on the reference element. It is transformed to the global elements using a non-linear mapping ${\bm{\Phi}}$, which includes the NURBS definition. The mapping — termed Triangle-Rectangle-Triangle (TRT) mapping — is given as:

[TABLE]

Here, $s$ and $r$ denote the coordinates in the reference triangle. ${\bf x}_{2}$ is the global coordinate of the iterior node and ${\Theta}_{1,3}$ are the NURBS coordinates of the FE boundary nodes.

The derivation of the mapping is illustrated in Figure 7. It is used for both the definition of the shape function and the placement of the quadrature points. Figure 8 compares the p-FEM mapping from [38] with the mapping in Equation (30). In principle, the two mappings lead to a similar — although still different — distribution of quadrature points. Note however, that the p-FEM mapping has a singularity at one of the boundary nodes; thus excluding this point in case of a boundary integral. The TRT mapping has its singularity at the interior node.

5 Numerical Methods: Coupling Approach

In this work, we adopt a partitioned coupling approach, where the fluid and the structure are treated as individual fields and solved separately. The coupling conditions are explicitly incorporated as a means for interchanging information between the structure and the fluid. Due to the use of individual solvers, one will usually encounter discrepancies between the discretizations of the structure and the fluid: starting from different levels of refinement, through different interpolation orders, to possibly even completely different discretization techniques. As a consequence, there is a demand for a strategy to interchange information/quantities between the individual fields. The most common approach is the "Neumann/Dirichlet" load transfer [24]: Forces resulting from the fluid boundary stresses are projected onto the structure as a Neumann boundary condition, while the structural deformations are transferred to the fluid as a Dirichlet boundary condition.

Two main aspects have to be considered:

temporal coupling: synchronization of individual fields,

• 2.

spatial coupling: transfer of loads and deformation.

5.1 Temporal Coupling

For the temporal coupling, we employ a strong coupling approach. Its advantage is that – in contrast to the often utilized weak coupling – the coupling conditions are fulfilled after each time step. This is achieved via fixed-point iterations between the structure and fluid within one time step until convergence.

A schematic description of the coupling strategy is given in the following (see Fig. 9):

compute an initial guess of the deformation $\tilde{\bf{d}}^{n+1}_{k}$ based on given quantities at $t_{n}$ and transfer it to the fluid, 2. 2.

deform the fluid domain based on $\tilde{\bf{d}}^{n+1}_{k}$ and calculate new fluid solution ${\bf{u}}^{n+1}_{k}$, $p^{n+1}_{k}$, 3. 3.

transfer stresses based on ${\bf{u}}^{n+1}_{k}$, $p^{n+1}_{k}$ onto structure, 4. 4.

use the stresses to compute ${\bf{d}}^{n+1}_{k+1}$, 5. 5.

test convergence $\lVert{\bf{d}}^{n+1}_{k+1}-{\bf{d}}^{n+1}_{k}\rVert\textless\epsilon$. If not converged go back to 2 and compute ${\bf{u}}^{n+1}_{k+1}$, $p^{n+1}_{k+1}$.

5.2 Spatial Coupling

For spatial coupling, we utilize a partitioned approach with individual solvers for both the structure and the fluid. As a consequence, information has to be transferred from one solver to the other, which then enters the respective other simulation as a boundary condition: From structure to fluid these are the stresses, from fluid to structure these are deformations. In this work, we have chosen two possible approaches for spatial coupling: (1) the finite interpolation method, which is also applicable to standard FEM, and (2) direct integration of the discrete forces.

5.2.1 Finite Interpolation Method

The general idea of this method is to transfer the discrete force at every fluid node onto its corresponding base point on the structure, and subsequently use the local non-zero basis function of the structural discretization to distribute the load to the corresponding nodes of the structure. In order to transfer the deformation back to the fluid domain, this procedure can be reversed. The individual steps are detailed in the following:

1) Calculation of discrete forces:

The method requires discrete forces at each fluid node at the interface $\Gamma_{FS}$. Therefore, the discrete stresses calculated by the fluid solution have to be integrated over the interface faces and distributed to each node. Given the discrete stress $\bm{\sigma}^{f}_{i}$ at each node on the boundary, the stress distribution is defined as follows:

[TABLE]

where $L_{i}$ denotes the Lagrangian shape functions and $r$ is the parametric coordinate along the interface. In order to regain the discrete force at each individual node, the stress has to be multiplied by the corresponding shape function and then integrated, resulting in the formula

[TABLE]

where $\bf{n}$ is the normal vector at the interface.

*2) Structural base point computation:

*In the next step, the base point on the structure of each fluid interface node has to be found. With our interface representation, this is straightforward: Since the spline representations are identical for both fluid and structure, the base point is identified through the local spline coordinate $\Theta_{j}$ associated with the specific NEFEM node. During preprocessing, it can, e.g., be calculated by using an orthogonal projection of a point onto a spline based on [40]. For multiple patches, the corresponding patch has to be identified before the parametric coordinate is computed.

3) Transfer of forces onto the structure:

We want to obtain the discrete force at every single control point of the spline. Using its NURBS basis function $R(\Theta)$, we can simply sum over all discrete forces multiplied by the basis function evaluated at the corresponding local parametric coordinate with:

[TABLE]

Here, ${nn^{f}_{\Gamma_{FS}}}$ stands for the number of fluid nodes on the interface. The NURBS basis functions fulfill the partition of unity property, meaning that

[TABLE]

holds. Therefore, the overall force acting on the interface is conserved during the projection.

*4) Transfer of deformation onto the fluid:

*The computed structural deformation has to be transferred back to the fluid in order to account for the mesh deformation. Once again, this procedure is straightforward due to the use of IGA. The deformation is computed at every spline control point. Knowing the parametric coordinates, we have to sum up the contributions of all control points that have non-zero NURBS basis function at that point. This is done in the following way:

[TABLE]

Here, ${\bf{d}}_{j}^{s}$ is the computed displacement on the structural side — a displacement applied to the control points —, whereas ${\bf{d}}_{i}^{f}$ is the displacement of fluid node $i$ — a displacement applied on FE node level.

5.2.2 Direct Transfer

If NEFEM is applied on the fluid side and IGA on the structural side, the geometry of the coupling interface is identical on both sides. Thus, a direct integration can be applied to obtain the fluid forces in sense of a weighted residual method [41].

1) Calculation of discrete forces:

Using the NURBS basis functions as test functions, the right-hand-side of the structural problem can be formulated as follows:

[TABLE]

The fluid stresses $\bm{\sigma}^{f}$ can be computed from the fluid solution and the face normals from the fluid grid. For NEFEM, the face normal can be determined exactly by evaluating the NURBS. The boundary integral can be compute piecewise by computing it on every fluid grid face:

[TABLE]

By applying Gaussian integration with $n_{GP}$ being the number of Gauss points we obtain

[TABLE]

If we insert now the Lagrangian representation of the fluid solution, the forces acting on the right-hand-side of the structural problem can be evaluated as

[TABLE]

Following this idea, the forces are already available on every control point and can be used as a right-hand-side for the system of equations of the structural problem. An additional projection method is not needed anymore.

Steps 2) and 4) remain the same as in section 5.2.1. Because the numerical integration is conducted on the same geometry, the resulting method is consistent and conservative [41].

6 Numerical Examples

In order to evaluate the performance of the new fluid-structure interaction approach, we have chosen three different test cases. The first was proposed by Schäfer et al.  [42], the second one by Wall [23] and the third one by Turek et al. [43]. The first two test cases involve steady laminar incompressible channel flow around a rigid and elastic cylinder respectively. The third test case is unsteady lamir flow around a structure that oscillates due to flow-induced vibrations. Thereby, the adapted NEFEM method with direct coupling can be tested for simple problems. What we expect to see is that the convergence rate will remain the same as for standard FEM, as we still utilize the same shape function definition, but the absolute value of the error will decrease. The latter is due to a decrease in geometrical error due to integration over a smooth and —- in the rigid case — even exact domain.

6.1 Cylinder

The first test case is used to compare the influence of using NEFEM in contrast to FEM and involves laminar flow around an attached, rigid cylinder. The configuration is proposed on the featflow benchmark site [42]. In particular this is the “DFG flow a around cylinder benchmark 2D-1”. A sketch of the computational domain is provided in Figure 10.

A parabolic inflow velocity is given at the inlet:

[TABLE]

No-slip conditions are applied to the upper and lower wall as well as to the cylinder and a free-flow boundary condition to the outlet. All parameters used for the simulation can be found in Table 2. The Reynolds number is defined as

[TABLE]

where $\overline{U}$ is the average inflow velocity.

For the given Reynolds number, the flow is steady. Therefore, the steady Navier-Stokes equations are solved.

An unstructured grid with 50 boundary elements on the cylinder and 5411 elements in total is used as a starting point. Four further grids are considered, which are obtained by a simple h-refinement of the initial grid. This is done by splitting every edge in half. An overview over the utilized grids is given in Table 2. In case of NEFEM, the geometry of the cylinder is represented by a second order NURBS line with 10 control points.

The reference value that serves as a basis for the comparison between the results for NEFEM and standard FEM is the drag coefficient. In both cases — NEFEM and FEM — the drag coefficient computed on the final grid (NEFEM: $c_{d}=5.570928$ and FEM: $c_{d}=5.569983$) is in good agreement with the reference value of $c_{d}=5.579535$ given by featflow. We therefore consider both implementations to be valid for this test case. In addition, Figure 11 gives a convergence plot for the drag obtained with the two methods. The featflow result is employed as a reference result. As expected, both methods feature the same convergence rate. However, one can observe that NEFEM slightly improves the error constant.

6.2 Flexible Cylinder

The second test case can be viewed as an extention of test case 6.1. Instead of a rigid cylinder, it now considers an elastic structure [23]. The dimensions of the domain are detailed in Fig. 12. The inflow velocity is constant; on the lower and upper wall a slip condition is applied. A no-slip boundary condition is used on the cylinder. The mass enclosed by the cylinder is neglected. The elastic structure is completely fixed on its rightmost point (on the horizontal symmetry axis). In order to avoid asymmetric results induced by numerical errors, the displacement of the leftmost cylinder point is suppressed in y-direction.

The given Reynolds number of 40 leads to a steady solution and a significant deformation of the cylinder. Pressure distributions around the rigid and deformed cylinder are given in Fig. 13.

On the fluid side, an unstructured grid with 40 boundary elements on the cylinder and 2352 elements in total is generated initially. Similar to the first test case, six additional grids are created by applying simple h-refinement to the initial grid. An overview of the grids used is given in Table 4. The geometry of the structure is defined by a second-order NURBS with 40 elements in circular direction and 5 elements in radial direction. Because the effect of NEFEM on FSI should be isolated, the discretization of the elastic structure remains unchanged within the presented study.

Two reference quantities are used to assess the convergence: (1) drag on the cylinder and (2) the displacement of the leftmost cylinder point. The results are presented as relative results to grid 6 in Fig. 14. As in the rigid case, the slope of the graphs are almost identical, but the relative error of NEFEM is clearly lower. In the presented case, the offset between both graphs is in the order of one refinement level.

6.3 Moving Flag

The third testcase is used in order to show the functioning of the method for transient FSI problems. The configuration is proposed in [43] and is referred to as "FSI benchmark - FSI2 ". A sketch of the computational domain is given in Fig. 15. The deformable elastic structure, or "flag", is fixed at the flag-cylinder interface.

Similar to Section 6.1 a parabolic velocity profile is set at the inlet. No-slip conditions are applied to the upper and lower wall as well as to the cylinder and the flag and a free-flow boundary conditions is applied to the outlet. All parameters relevant for the simulations are shown in Table 6. For the given Reynolds number of 100 the interaction of the flow and the elastic structure causes self-induced oscillations of the structure which results in a quasi-steady periodic wave-like deformation of the flag. The pressure distribution at time t=13.6 s is given in Fig. 16.

We discretize the fluid domain using an unstructured grid with 75 boundary elements on the flag surface and 3572 elements in total. Similar to the previous test cases we create two other grids by applying simple h-refinement to the intitial grid. As reference for the convergence study we use the solution computed on a fine grid with 584 elements on the flag surface. An overview over the grids used is given in Table 6. The structure is discretized using a quadratic NURBS with 60x15 elements. It is kept constant throughout the study. In all computations a time step of $\Delta t=0.002s$ was applied.

The quantity for comparison is the amplitude of the vertical displacement of point A, located at the tail of the flag. Fig. 17 shows the results as relative error to the reference solution. As in the previous test cases, the offset between both graphs is identical, but again the relative error of NEFEM is clearly lower.

7 Conclusion

The paper presents a novel coupling scheme for fluid-structure-interaction problems. Building on the recent advances in structural mechanics based on isogeometric analysis, the idea is to involve the splines also in the flow simulations. This was achieved by employing the NURBS-Enhanced Finite Element Method for the fluid. A new geometrical mapping from physical to reference space ensuring that Dirichlet boundary conditions are fulfilled on the spline-based surface was presented for NEFEM. With that, a direct transfer of the necessary coupling variables is possible, as the interface description is identical for both the structure and the fluid. In addition to the simplified implementation of the coupling, an increased accuracy of the flow solution — and with this the FSI solution — was to be expected due to the matching interface representation. This expected improvement was confirmed with steady computations on a rigid and deformable cylinder as well as transient computations of flow-induced vibrations, see Section 6. Even if the convergence rate remains the same, the resulting error of NEFEM compared to standard finite elements can be decreased. With the second and third example it was shown that the method is directly applicable to elastic configurations using isogeometric analysis on the structural side. Finally, the presented method can be used to generate fluid grids for more complex problems with existing tools and allows for improved accuracy with minor implementation effort.

8 Acknowledgements

The authors gratefully acknowledge the support of the German research foundation DFG grant “Geometrically Exact Methods for Fluid-Structure Interaction”. Furthermore, we thank Daniel Hilger for his support and ideas.