# A rheological interface model and its space-time finite element formulation for fluid-structure interaction

Antoine Legay, Andreas Zilian, Christian Janssen

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

## To cite this version:

Antoine Legay, Andreas Zilian, Christian Janssen. A rheological interface model and its space-time finite element formulation for fluid-structure interaction. International Journal for Numerical Methods in Engineering, 2010, 86, pp.667 - 687. 10.1002/nme.3060. hal-01422440

HAL Id: hal-01422440

https://hal.science/hal-01422440

Submitted on 25 Dec 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>

Distributed under a Creative Commons Attribution 4.0 International License
# A rheological interface model and its space–time finite element formulation for fluid–structure interaction

Antoine Legay¹, Andreas Zilian² and Christian Janssen³

¹Structural Mechanics and Coupled Systems Laboratory, Conservatoire National des Arts et Métiers, 2 Rue Conté, 75003 Paris, France

²Institut für Statik, Technische Universität Braunschweig, Beethovenstraße 51, 38106 Braunschweig, Germany

³Institute for Computational Modeling in Civil Engineering, Technische Universität Braunschweig, Mühlenpfordstr. 4-5, 38106 Braunschweig, Germany

This contribution discusses extended physical interface models for fluid-structure interaction problems and investigates their phenomenological effects on the behavior of coupled systems by numerical simulation. Besides the various types of friction at the fluid-structure interface the most interesting phenomena are related to effects due to additional interface stiffness and damping. The paper introduces extended models at the fluid-structure interface on the basis of rheological devices (Hooke, Newton, Kelvin, Maxwell, Zener). The interface is decomposed into a Lagrangian layer for the solid-like part and an Eulerian layer for the fluid-like part. The mechanical model for fluid-structure interaction is based on the equations of rigid body dynamics for the structural part and the incompressible Navier-Stokes equations for viscous flow. The resulting weighted residual form uses the interface velocity and interface tractions in both layers in addition to the field variables for fluid and structure. The weak formulation of the whole coupled system is discretized using space-time finite elements with a discontinuous Galerkin method for time-integration leading to a monolithic algebraic system. The deforming fluid domain is taken into account by deformable space-time finite elements and a pseudo-structure approach for mesh motion. The sensitivity of coupled systems to modification of the interface model and its parameters is investigated by numerical simulation of flow induced vibrations of a spring supported fluid-immersed cylinder. It is shown that the presented rheological interface model allows to influence flow-induced vibrations.

KEY WORDS: fluid-structure interaction; rheological interface model; space-time finite element method

## 1. INTRODUCTION

Real-life situations of fluid-structure interaction may involve complex physical effects at the coupling interface that separates the fluid domain from the structural part. Besides the various types of friction at the fluid-structure interface, the most interesting and effective phenomena are related to additional interface stiffness and damping. Interactions between fluid and structure occur in a wide range of problems of structural engineering. In the context of noise and vibration reduction, engineers have developed passive damping material. These materials may be used to
design effective damping devices and to shift the eigenfrequencies of the coupled system appropriately. Many approaches have been proposed in the literature to model the absorbing material. They are often based on poroelastic theory [1–4] and have to be discretized as three-dimensional bodies, which leads to a large number of degrees of freedom. In the design process, parametric studies are often needed and are accompanied by reduced models to provide efficient computational methods. This approach has been successfully used for vibro-acoustic problems where the fluid is an acoustic fluid described by the pressure in [5, 6] or described by both pressure and displacement in [7].

In this work, the rheological interface is divided into two layers. The first, namely the Lagrangian layer, models the solid part of the interface in which the fluid is stagnant. The second, namely the Eulerian layer, models the fluid part of the interface where the fluid behaves as an incompressible flow along the solid part. For instance, a foam or porous medium interface would be seen as two parts: (i) a part where the fluid is within closed solid cavities (Lagrangian layer) and (ii) a part where the fluid flows in and along open solid cavities (Eulerian layer). The developed rheological model allows stiffness and/or damping in both normal and tangential directions to the interface, which permits a large number of possibilities to model a given material. For instance, a damping in normal direction for the Eulerian layer could describe permeability of a tissue. The aim of this paper is to provide a numerical strategy for reduced-order interface description in the context of fluid-structure interaction problems. For this purpose, a generalized constitutive interface law is assumed without restriction to specific interface physics.

The general governing equations of the extended interface physics have to be integrated in a numerical solution scheme for fluid-structure interaction. *Partitioned approaches* solve the involved equations sequentially for each of the subproblems and are designed for highly modular application of specific discretization and solution methods. The interaction between the fluid and structural subsystems is accomplished by interface conditions. Weakly coupled solutions in the time domain are carried out if exchange of information between the subsystems takes place once per time step only [8–11], and therefore, the overall system will not be in a physically consistent state. While for weakly coupled systems such an approach is acceptable, it cannot be applied successfully for numerical investigations of strong couplings. Strongly coupled solutions are obtained from partitioned methods by iterating between the solutions of the subsystems within one time step [12, 13], although hereby a consistent state of the overall system can often only be achieved by under-relaxation in block Gauss-Seidel and Dirichlet-Neumann methodologies [14, 15]. Alternative partitioned schemes utilize exact or inexact Newton strategies [16–21]. *Monolithic approaches* do not partition the overall physical domain so that all model equations are considered in a single coupled system of equations. This monolithic method together with a consistent time discretization is an ideal basis for conservative coupling formulations and guarantees the transient-precise coupling. For physically strongly coupled systems such schemes may lead to improved convergence and numerical efficiency compared to partitioned approaches [22, 23].

In this work, a monolithic approach for fluid and structural domain is used, leading to an implicit coupling of the two subdomains including localized interface effects. The weighted residual method is applied to the set of governing equations, including the two-layer rheological interface model, and the resulting weak form is discretized using a time-discontinuous Galerkin method based on space-time finite elements. Space-time finite elements for the solution of the incompressible Navier-Stokes equations on moving meshes have been introduced by Tezduyar et al. [24]. For applications in elastodynamics, Hughes and Hulbert [25] introduced time-discontinuous, stabilized space-time (SST) finite elements. Tezduyar et al. developed the deforming-spatial-domain/stabilized space-time (DSD/SST) finite element method for fluid-structure interaction problems. The space-time finite element discretization of a velocity-based monolithic model for strongly coupled fluid-structure interaction problems has been shown to be eligible for a number of different applications [26–29].

The sensitivity of coupled systems to modification of the interface model and its parameters is investigated by numerical simulation of flow induced vibrations of a spring-supported cylinder immersed in laminar incompressible flow.
## 2.1. Generalized rheological model

The fluid-structure interaction problem is modeled by a fluid domain $\Omega_F$ with an immersed structural domain $\Omega_S$. At the common interface $\Sigma = \Omega_F \cup \Omega_S$ of the boundary-coupled system, a lower order manifold is defined with the normal vector $\mathbf{n}$ and tangential vectors $\mathbf{s}_i$ as shown in Figure 1. Both continua, fluid and structure, are assumed to be described in terms of velocity $\mathbf{v}_F$ and $\mathbf{v}_S$.

In order to be able to represent different types of rheological devices, a generalized three-parameter model is utilized. Each geometrical direction may have its own set of parameters denoted by subscript $n$ for normal and $s$ for tangential direction. The interface is assumed to be in between two bodies $\Omega_S$ and $\Omega_F$. By convention, the normal vector $\mathbf{n}$ is pointing from $\Omega_S$ to $\Omega_F$ (Figure 1). The stress, displacement, and velocity vectors are defined in the local coordinate interface basis ($\mathbf{n}, \mathbf{s}_1, \mathbf{s}_2$). We assume that the general constitutive relation at the interface is linear [30] and written as

$$t + \alpha \dot{t} - \beta \Delta u - \gamma \Delta v = 0. \quad (1)$$

This relation is a differential equation in time where $\Delta u$ and $\Delta v$ are, respectively, the displacement and velocity jumps at the interface. Herein the coefficient matrices $\alpha, \beta$, and $\gamma$ denote the structural material properties of the rheological interface which depend on the thickness. In the general case these are fully occupied, leading to a coupled behavior in the different basis directions. In the following we assume decoupled material behavior such that the coefficient matrices given by

$$\alpha = \begin{bmatrix} \alpha_n & 0 & 0 \\ 0 & \alpha_s & 0 \\ 0 & 0 & \alpha_s \end{bmatrix}, \quad \beta = \begin{bmatrix} \beta_n & 0 & 0 \\ 0 & \beta_s & 0 \\ 0 & 0 & \beta_s \end{bmatrix}, \quad \gamma = \begin{bmatrix} \gamma_n & 0 & 0 \\ 0 & \gamma_s & 0 \\ 0 & 0 & \gamma_s \end{bmatrix} \quad (2)$$

are diagonal and regular in the local coordinate interface basis.

The actual choice of $\alpha, \beta$, and $\gamma$ determines the specific active rheological model in each local interface direction. The switch into the model choice can be made according to Table I and Figure 2 depicts associated specializations of the generalized equilibrium relation. The spring and damper devices used in the rheological model have the inherent spatial extent of interface thickness $h$. The unit of the stiffness parameters $\beta_n$ and $\beta_s$ is defined as Pascal per unit length (physical interface thickness) while the damping parameters $\gamma_n$ and $\gamma_s$ are expressed in terms of Pascal·second per unit length (physical interface thickness). The unit of the parameters $\alpha_n$ and $\alpha_s$ is second.

Determination of the parameters for a given interface material is subjected to experimental or numerical investigations. For example, in Section 5.2.1 of [6], parameters of a Kelvin-Voigt model are identified from a full three-dimensional numerical analysis of a poro-elastic material. On the other hand, experiments on small samples can provide characteristic curves allowing a parameter identification as is done in [5] using results from [31].

<figure><img src="image_3.png" /><figcaption>Figure 1. Fluid-structure interface with local basis and stress vector.</figcaption></figure>
<table><caption>Table I. Rheological models in the interface of thickness h.</caption><thead><tr><th>Model</th><th>α [s]</th><th>h [m] · β [Pa·m<sup>-1</sup>]</th><th>h [m] · γ [Pa·s·m<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Hooke (a)</td><td>0</td><td>E</td><td>0</td></tr><tr><td>Newton (b)</td><td>0</td><td>0</td><td>η</td></tr><tr><td>Kelvin (c)</td><td>0</td><td>E</td><td>η</td></tr><tr><td>Maxwell (d)</td><td>η/E</td><td>0</td><td>η</td></tr><tr><td>Zener-M (e)</td><td>η/E<sub>1</sub></td><td>E<sub>2</sub></td><td>η(E<sub>1</sub>+E<sub>2</sub>)/E<sub>1</sub></td></tr><tr><td>Zener-K (f)</td><td>η/(E<sub>1</sub>+E<sub>2</sub>)</td><td>E<sub>1</sub>E<sub>2</sub>/(E<sub>1</sub>+E<sub>2</sub>)</td><td>ηE<sub>2</sub>/(E<sub>1</sub>+E<sub>2</sub>)</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 2. Different possible rheological interface models.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. Concept of a two-layer rheological interface.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4. Model of the two-layer interface.</figcaption></figure>

## 2.2. *Composition of interface by Lagrange–Euler layers*

In the proposed rheological interface model, the interface links a solid body in Lagrangian description to a fluid body in Eulerian description. Therefore two layers (Figure 3) are introduced in the interface domain: the Lagrangian layer (subscript *L*) and the Eulerian layer (subscript *E*). By this construction, an independent model can be constructed in each layer. Moreover, this enables the Lagrangian layer to represent a visco-elastic solid-like behavior as well as a mechanical behavior of a stagnant compressible fluid, for instance such as in foams or porous media. In contrast, the Eulerian layer is thought of containing a flowing media which is for a large number of applications a viscous incompressible flow or a creeping flow. In the present approach, introduction of rheologic interface models is accompanied by the interface velocity $v_i$ as additional unknown as well as the tractions in the Lagrangian layer $t_L$ and in the Eulerian layer $t_E$ (Figure 4).
## 2.3. Strong form interface momentum balance

The stress continuity at the boundary between the Lagrangian and Eulerian layers

$$ \mathbf{t}_E - \mathbf{t}_L = \mathbf{0} \qquad (3) $$

ensures local conservativity of interface forces.

## 2.4. Strong form governing equation in the Lagrangian layer

The kinematic consistency equations written in terms of relative velocity for the Lagrangian layer are obtained by multiplying Equation (1) by $\beta^{-1}$ and taking its time derivative:

$$ \beta_L^{-1} \dot{\mathbf{t}}_L + \alpha_L \beta_L^{-1} \ddot{\mathbf{t}}_L - \Delta \mathbf{v}_L - \gamma_L \beta_L^{-1} \Delta \dot{\mathbf{v}}_L = \mathbf{0} \qquad (4) $$

with

$$ \Delta \mathbf{v}_L = \mathbf{v}_i - \mathbf{v}_S. \qquad (5) $$

The rate formulation assumes non-zero and time-independent parameters $\beta_{nL}$ and $\beta_{sL}$.

## 2.5. Strong form governing equation in Eulerian layer

Since in the Eulerian layer adjacent material particles do not necessarily remain neighbors during deformation, the following relation applies, neglecting $\alpha$ and $\beta$ terms in Equation (1):

$$ \gamma_E^{-1} \mathbf{t}_E - \Delta \mathbf{v}_E = \mathbf{0} \qquad (6) $$

with

$$ \Delta \mathbf{v}_E = \mathbf{v}_F - \mathbf{v}_i. \qquad (7) $$

The matrix $\gamma_E$ contains viscous damping $\gamma_{sE}$ in tangential direction and $\gamma_{nE}$ in normal direction. The tangential component can be linked to sliding friction while the normal component can be used for interface permeability as it is, e.g., for membrane tissues.

# 3. SPACE-TIME WEAK FORM OF THE RHEOLOGICAL INTERFACE

A space-time discretization uniformly includes the temporal axis in addition to the spatial dimensions in a finite element approach. For numerical efficiency, the interface space-time domain $\mathcal{R}$ is divided into a sequence of $N$ time slabs $\mathcal{R}^n = \sum^n \times [t_n, t_{n+1}]$, as shown in Figure 5, which are solved successively. The spatial discretization at the fluid-structure interface is assumed to be such that a fitting mesh is realized (Figure 5(a)).

At time instant $t_n$ the energy of the discretized system at the end of the previous time slab $t_n^-$ has to be equal to the energy at the beginning of the next time step $t_n^+$ (Figure 5(b)). For time-discontinuous approximations of field unknowns this leads to additional jump terms in the weak form. Moreover, spatial discretizations from $t_n^-$ and $t_n^+$ do not need to be conforming. For first-order ordinary differential equations the resulting time integration scheme is A-stable and third-order accurate for linear temporal interpolation.

## 3.1. Weak form of interface momentum balance

The interface momentum balance is weakly enforced by multiplying Equation (3) by the interface velocity weighted function $\delta \mathbf{v}_i$ supposed to be sufficiently smooth. By integrating over the space time slab $\mathcal{R}^n$, this leads to

$$ \int_{\mathcal{R}^n} \delta \mathbf{v}_i \cdot (\mathbf{t}_E - \mathbf{t}_L) d\mathcal{R} = 0 \quad \forall \delta \mathbf{v}_i. \qquad (8) $$
<figure><img src="image_7.png" /><figcaption>Figure 5. Discretization in space and time: (a) fluid-structure interface in the space-time domain and (b) time-discontinuous shape functions.</figcaption></figure>

## 3.2. Weak form of the Lagrangian layer

The weighted residual formulation of the Lagrangian interface layer is obtained from Equation (4) by choosing the Lagrangian interface traction $\delta \mathbf{t}_L$ as test function and integrating over the space time slab $\mathcal{R}^n$. Performing an integration by parts in time on the term $\int_{\mathcal{R}^n} \delta \mathbf{t}_L \cdot \boldsymbol{\alpha}_L \boldsymbol{\beta}_L^{-1} \dot{\mathbf{t}}_L \, d\mathcal{R}$, the second-order time derivative is avoided and linear shape functions can be used in time direction. Since a time-discontinuous Galerkin scheme is used, the continuity of state variables is weakly enforced by introducing additional jump terms to the integral form of the interface equation.

Using relation (5), the resulting weak form of the Lagrangian part of the two-layer interface model then reads as

$$ \int_{\mathcal{R}^n} \delta \mathbf{t}_L \cdot \boldsymbol{\beta}_L^{-1} \dot{\mathbf{t}}_L \, d\mathcal{R} - \int_{\mathcal{R}^n} \delta \mathbf{t}_L \cdot \mathbf{v}_i \, d\mathcal{R} + \int_{\mathcal{R}^n} \delta \mathbf{t}_L \cdot \mathbf{v}_S \, d\mathcal{R} \quad (9a) $$

$$ - \int_{\mathcal{R}^n} \delta \mathbf{t}_L \cdot \gamma_L \boldsymbol{\beta}_L^{-1} \dot{\mathbf{v}}_i \, d\mathcal{R} + \int_{\mathcal{R}^n} \delta \mathbf{t}_L \cdot \gamma_L \boldsymbol{\beta}_L^{-1} \dot{\mathbf{v}}_S \, d\mathcal{R} \quad (9b) $$

$$ - \int_{\mathcal{R}^n} \delta \dot{\mathbf{t}}_L \cdot \boldsymbol{\alpha}_L \boldsymbol{\beta}_L^{-1} \dot{\mathbf{t}}_L \, d\mathcal{R} + \int_{\Sigma_{n+1}} \delta \mathbf{t}_L(t_{n+1}^-) \cdot \boldsymbol{\alpha}_L \boldsymbol{\beta}_L^{-1} \dot{\mathbf{t}}_L(t_{n+1}^-) \, d\Sigma \quad (9c) $$

$$ - \int_{\Sigma_n} \delta \mathbf{t}_L(t_n^+) \cdot \boldsymbol{\alpha}_L \boldsymbol{\beta}_L^{-1} \dot{\mathbf{t}}_L(t_n^+) \, d\Sigma + \int_{\Sigma_n} \delta \mathbf{t}_L(t_n^+) \cdot \boldsymbol{\beta}_L^{-1} [\mathbf{t}_L(t_n^+) - \mathbf{t}_L(t_n^-)] \, d\Sigma \quad (9d) $$

$$ - \int_{\Sigma_n} \delta \mathbf{t}_L(t_n^+) \cdot \gamma_L \boldsymbol{\beta}_L^{-1} [\mathbf{v}_i(t_n^+) - \mathbf{v}_i(t_n^-)] \, d\Sigma + \int_{\Sigma_n} \delta \mathbf{t}_L(t_n^+) \cdot \gamma_L \boldsymbol{\beta}_L^{-1} [\mathbf{v}_S(t_n^+) - \mathbf{v}_S(t_n^-)] \, d\Sigma \quad (9e) $$

$$ - \int_{\Sigma_n} \delta \dot{\mathbf{t}}_L(t_n^+) \cdot \boldsymbol{\alpha}_L \boldsymbol{\beta}_L^{-1} [\mathbf{t}_L(t_n^+) - \mathbf{t}_L(t_n^-)] \, d\Sigma = 0 \quad \forall \delta \mathbf{t}_L. \quad (9f) $$

## 3.3. Weak form of the Eulerian layer

The weighted residual formulation of the Eulerian layer is obtained by multiplying Equation (6) by the Eulerian interface traction $\delta \mathbf{t}_E$ and integrating over the space time slab $\mathcal{R}^n$. Using relation (7), the resulting weak form of the Eulerian part of the two-layer interface model then reads as

$$ \int_{\mathcal{R}^n} \delta \mathbf{t}_E \cdot \gamma_E^{-1} \mathbf{t}_E \, d\mathcal{R} - \int_{\mathcal{R}^n} \delta \mathbf{t}_E \cdot \mathbf{v}_F \, d\mathcal{R} + \int_{\mathcal{R}^n} \delta \mathbf{t}_E \cdot \mathbf{v}_i \, d\mathcal{R} = 0 \quad \forall \delta \mathbf{t}_E. \quad (10) $$

## 3.4. Transformation to the global basis

The velocities for fluid and solid bodies are usually described in the global frame of reference, while the interface conditions and rheological models considered in the previous derivations are
defined in a local interface coordinate system. A symmetric transformation matrix from the local interface basis to the global frame of reference is introduced for both velocity and stress vectors.

# 4. SPACE-TIME DISCRETIZATION OF THE RHEOLOGICAL INTERFACE

## 4.1. Discretization of the weak form in space and time

The approximation space for the velocity at the coupling boundaries is determined by the chosen approximations of the solid and fluid domains, respectively. In other words, the discretized continuous solid boundary velocity $v_S^h$ and fluid boundary velocity $v_F^h$ are given by the space-time finite element discretization in the domain as is shown in Figure 5(a). The additional velocity at the interface $v_i^h$ as well as the tractions $t_L^h$ and $t_E^h$ use the same approximation space as structure and fluid velocity leading to an LBB-stable discrete formulation, see Figure 6. In the following, an approximation of the interface field variables in terms of space-time ansatz functions is chosen:

$$v_i^h(\mathbf{x}, t) = N_i(\mathbf{x}, t)\hat{v}_i, \quad t_L^h(\mathbf{x}, t) = N_L(\mathbf{x}, t)\hat{t}_L, \quad t_E^h(\mathbf{x}, t) = N_E(\mathbf{x}, t)\hat{t}_E, \qquad (11)$$

where $\hat{v}_i$, $\hat{t}_L$, and $\hat{t}_E$ are the vectors of nodal values of the approximation related to the shape function given by $N_i$, $N_L$, and $N_E$ for the interface, the Lagrangian layer and the Eulerian layer, respectively.

## 4.2. Verification of the rheological interface model

Verification of the developed rheological interface is performed taking into account the purely time-dependent setup depicted in Figure 7. This setup consists of the Lagrangian interface layer and an attached mass $m_s$. The Eulerian layer is fixed such that $v_E = 0 = v_i$. The Lagrangian traction $t_L = 2e_1 H(4-t)$ describes an instant discontinuous loading function. In this work, $H(t)$ denotes the Heaviside function. Several combinations of material parameters are studied according to Table II,

<figure><img src="image_8.png" /><figcaption>Figure 6. Spacial discretization at the rheological interface. Note that the interface has no thickness in the model while it has one on the figure for the sake of clarity.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7. Verification test setup with loading function.</figcaption></figure>
<table><caption>Table II. Data values for verification tests.</caption><thead><tr><th>Device</th><th>m<sub>S</sub></th><th>η</th><th></th><th>E<sub>1</sub></th><th></th><th>E<sub>2</sub></th></tr></thead><tbody><tr><td>Hooke</td><td>1</td><td>0</td><td></td><td>1</td><td>10</td><td>100</td><td>—</td></tr><tr><td>Kelvin</td><td>1</td><td>1</td><td></td><td>1</td><td>10</td><td>100</td><td>—</td></tr><tr><td>Newton</td><td>1</td><td>10</td><td>100</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Maxwell</td><td>0</td><td>100</td><td></td><td>2.5</td><td>5</td><td>10</td><td>—</td></tr><tr><td>Zener(K)-1</td><td>0</td><td>5</td><td></td><td>2.5</td><td>5</td><td>10</td><td>10</td></tr><tr><td>Zener(K)-2</td><td>0</td><td>5</td><td></td><td>10</td><td></td><td>2.5</td><td>5</td><td>10</td></tr></tbody></table>

representing five different rheological devices namely Hooke, Kelvin-Voigt, Newton, Maxwell, and Zener(K) devices. The time step size is chosen as $\Delta t = 0.01$. As shown in Figure 8, numerical results obtained by the proposed implemented method are in very good agreement with closed-form solutions to the time-dependent behavior associated to fundamental rheological device types. The time-discontinuous Galerkin method enables accurate representation of jumps in the temporal solution to the Maxwell device, see Figure 8(d) and the Zener-(K,M) device, see Figures 8(e) and (f), for the prescribed discontinuous loading function.

# 5. RHEOLOGICAL INTERFACE MODEL IN A FLUID-STRUCTURE PROBLEM

## 5.1. Model of the fluid flow

For an incompressible fluid of density $\rho$ the strong forms of the balance of momentum and mass conservation

$$ \rho(\mathbf{v}_{F,t} + \mathbf{v}_F \cdot \nabla \mathbf{v}_F - \mathbf{b}) - \nabla \cdot \mathbf{T} = \mathbf{0} \quad \text{in } \mathcal{Q}^F = \Omega^F \times I \quad (12) $$

$$ \nabla \cdot \mathbf{v}_F = 0 \quad \text{in } \mathcal{Q}^F \quad (13) $$

are defined on the current space-time domain $\mathcal{Q}^F$ of the flow and are completed by Dirichlet boundary conditions $\bar{\mathbf{v}}_F$ for the velocity and Neumann boundary conditions $\bar{\mathbf{t}}$ for the tractions on $P^F$ as well as initial conditions $\mathbf{v}_{F0} = \mathbf{v}_F(t_0)$, fulfilling the incompressibility constraint at $t_0$. Volume forces (e.g. gravity) can be introduced by the vector $\mathbf{b}$. For a Newtonian fluid, the stress tensor consists of the static pressure $p$ and the linear strain rate tensor $\mathbf{D}$ with the kinematic viscosity $\mu$,

$$ \mathbf{T} = 2\mu\mathbf{D} - p\mathbf{I}, \quad \mathbf{D} = \frac{1}{2}((\nabla\mathbf{v}_F) + (\nabla\mathbf{v}_F)^T). \quad (14) $$

The weighted residual formulation of the incompressible and viscous fluid using the Galerkin method reads as $\forall \delta\mathbf{v}_F, \delta p, \delta\mathbf{t}$

$$ \int_{\mathcal{Q}^F} \delta \mathbf{v}_F \cdot \rho (\mathbf{v}_{F,t} + \mathbf{v}_F \cdot \nabla \mathbf{v}_F - \mathbf{b}) d\mathcal{Q} + 2\mu \int_{\mathcal{Q}^F} \mathbf{D}(\delta \mathbf{v}_F) : (\mathbf{D}(\mathbf{v}_F) d\mathcal{Q} - \int_{\mathcal{Q}^F} \nabla \cdot \delta \mathbf{v}_F p \, d\mathcal{Q} \quad (15a) $$

$$ + \int_{\mathcal{Q}^F} \delta p \nabla \cdot \mathbf{v}_F d\mathcal{Q} \quad (15b) $$

$$ + \int_{\Omega^F} \delta \mathbf{v}_F (t_n^+) \cdot \rho (\mathbf{v}_F(t_n^+) - \mathbf{v}_F(t_n^-)) d\Omega \quad (15c) $$

$$ - \int_{P_v^F} \delta \mathbf{v}_F \cdot \mathbf{t} \, dP + \int_{P_v^F} \delta \mathbf{t} \cdot (\mathbf{v}_F - \bar{\mathbf{v}}_F) \, dP - \int_{P_t^F} \delta \mathbf{v}_F \cdot \bar{\mathbf{t}} \, dP \quad (15d) $$

$$ + \sum_e \int_{\mathcal{Q}_e^F} (\rho \delta \mathbf{v}_{e,t} + \rho \mathbf{v} \cdot \nabla(\delta \mathbf{v}) - \nabla \cdot (\delta \mathbf{T})) \cdot \tau_m \frac{1}{\rho} (\rho \mathbf{v}_{e,t} + \rho \mathbf{v} \cdot \nabla \mathbf{v} - \nabla \cdot \mathbf{T} - \rho \mathbf{b}) d\mathcal{Q} = 0. \quad (15e) $$
<figure><img src="image_10.png" /></figure>

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

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

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

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

<figure><img src="image_15.png" /><figcaption>Figure 8. Results for the verification of the rheological interface: (a) displacements for Hooke device; (b) displacements for Kelvin–Voigt device; (c) velocities for Newton device; (d) displacement for Maxwell device; (e) displacement for Zener-K device, $E_1$ varia; and (f) Displacement for Zener-K device, $E_2$ varia.</figcaption></figure>

Line (15a) represents the weak form of conservation of momentum for a Newtonian fluid with volume force vector **b**. The incompressibility constraint (15b) is weighted with the pressure. Line (15c) ensures the consistent transfer of momentum from the previous time slab end $t_n^-$ to the current one $t_n^+$. Prescribed velocities $\bar{v}_F$ and tractions $\bar{t}$ at Dirichlet and Neumann boundaries are present in variational form (15d) after partial integration. The weighted residual form is stabilized by a Galerkin/least squares term [32] of the momentum balance in line (15e). The Galerkin/least squares stabilization suppresses numerical oscillations in solutions to hyperbolic differential equations
by the introduction of additional numerical diffusion and allows the application of equal order approximations of velocities and pressure for the incompressible flow field. The stabilization parameter $\tau_m$ is determined for each element e using the definition given by Tezduyar et al. [24].

The approximation of the fluid field variables in terms of a space-time discretization using ansatz functions depending on the spatial coordinates and time is

$$ \mathbf{v}_F^h(\mathbf{x}, t) = \mathbf{N}_F(\mathbf{x}, t)\hat{\mathbf{v}}_F \quad p^h(\mathbf{x}, t) = \mathbf{N}_F(\mathbf{x}, t)\hat{p} \quad (16) $$

and follows the detailed implementation in [26]. The deforming fluid domain is taken into account by deformable space-time finite elements and a pseudo-structure approach for mesh motion.

## 5.2. Model of spring-supported rigid body

For investigation of vibration phenomena of less deforming bodies like sections of cables or bridge decks in aeroelastic applications, it is often sufficient to model an elastically supported rigid body instead of a deformable solid. The momentum balance of a rigid body is given by

$$ \mathbf{M}\dot{\mathbf{v}}_S + \mathbf{f} = \mathbf{p}(t) \quad \text{in } I, \quad (17) $$

where $\mathbf{v}_S$ contains the center of gravity velocity as well as the angular velocity of the solid, $\mathbf{M}$ is the mass matrix combined with the inertia matrix, $\mathbf{p}(t)$ the vector of resulting external forces, and $\mathbf{f}$ is the vector of linear spring forces governed by Hooke's law in a rate formulation

$$ \mathbf{K}^{-1}\dot{\mathbf{f}} - \mathbf{v}_S = \mathbf{0} \quad \text{in } I, \quad (18) $$

where $\mathbf{K}$ is the stiffness matrix of the spring forces. Herein, the acceleration of a material point is expressed by the temporal derivative of velocity $\mathbf{v}_S$ and gives rise to uniform and consistent formulation in terms of velocities at the boundary of the body. In the following application, the rotational degree of freedom is not considered but could be easily incorporated by inertia $\theta$ and angular velocity $\dot{\phi}$. For solution of the initial value problem, conditions $\mathbf{v}_{S0} = \mathbf{v}_S(t_0)$ on initial velocity are necessary.

Following [26], the time-discontinuous Galerkin formulation of the above model equations in the time slab $I_n$ reads as

$$ \int_{I_n} \delta \mathbf{v}_S \cdot (\mathbf{M} \dot{\mathbf{v}}_S + \mathbf{f} - \mathbf{p}) dt + \delta \mathbf{v}_S(t_n^+) \mathbf{M} (\mathbf{v}_S(t_n^+) - \mathbf{v}_S(t_n^-)) \quad (19a) $$

$$ + \int_{I_n} \delta \mathbf{f} \cdot (\mathbf{K}^{-1} \dot{\mathbf{f}} - \mathbf{v}_S) dt + \delta \mathbf{f}(t_n^+) \mathbf{K}^{-1} (\mathbf{f}(t_n^+) - \mathbf{f}(t_n^-)) = 0 \quad \forall \delta \mathbf{v}_S, \delta \mathbf{f}. \quad (19b) $$

The approximation of the rigid body variables uses ansatz functions in time only

$$ \mathbf{v}_S^h(t) = \mathbf{N}_S(t)\hat{\mathbf{v}}_S \quad \mathbf{f}^h = \mathbf{N}_S(t)\hat{\mathbf{f}} \quad (20) $$

and the velocity state at the outer boundary of the rigid body is obtained by first order projection of $\mathbf{v}_S$ at the center of gravity.

# 6. APPLICATION TO ONE-DIMENSIONAL TEST PROBLEM

## 6.1. Setup

The one-dimensional problem presented in Figure 9 validates the formulation for a fluid-structure interaction problem. This problem also enables a detailed description of the implementation. The Lagrangian layer of the interface is a Hooke device ($\beta \neq 0$) while the Eulerian layer is a Newton device ($\gamma \neq 0$). The structure is a spring (stiffness $k$)-supported mass $m$ which is fixed on the other extremity. The fluid (density $\rho$) is contained in a channel (length $l$) and is supposed to be incompressible. The coordinate system $x$ is attached to the channel such that $x=0$ is the contact
<figure><img src="image_16.png" /><figcaption>Figure 9. One-dimensional problem.</figcaption></figure>

point between the fluid and the Eulerian layer of the interface. The traction $\bar{t}$ is imposed on the right boundary at $x=l$. The whole system is initially at rest. The velocities of the interface, structure, and fluid are denoted by $v_i(t)$, $v_S(t)$, and $v_F(x, t)$, respectively. The pressure in the fluid is denoted $p(x, t)$. The traction in the Lagrangian layer, Eulerian layer and in the spring are denoted by $t_L$, $t_E$, and $f$, respectively.

## 6.2. Weak form

The weak form follows the description given in Section 3 for the interface, in Section 5.2 for the solid, and in Section 5.1 for the fluid. The notation is the same but degenerated to the one-dimensional case.

The weak form of this specific one-dimensional fluid-structure system reads for the spring-supported mass as

$$ \int_{I_n} \delta v_S(m \dot{v}_S + f - t_L) dt + \delta v_S(t_n^+) m (v_S(t_n^+) - v_S(t_n^-)) = 0 \quad \forall \delta v_S \quad (21) $$

$$ \int_{I_n} \delta f \left( \frac{1}{k} \dot{f} - v_S \right) dt + \delta f(t_n^+) \frac{1}{k} (f(t_n^+) - f(t_n^-)) = 0 \quad \forall \delta f \quad (22) $$

for the simplified fluid (no convective term is taken into account)

$$ \int_{[0 \ l] \times I_n} (\delta v_F \rho \dot{v}_F \delta v_{F,x} p) dx dt - \int_{I_n} (\delta v_F(0) t_E \delta v_F(l) \bar{t}) dt + \int_{[0 \ l]} \delta v_F(t_n^+) \rho (v_F(t_n^+) - v_F(t_n^-)) dx = 0 \quad \forall \delta v_F, \quad (23) $$

$$ \int_{[0 \ l] \times I_n} \delta p v_{F,x} dx dt = 0 \quad \forall \delta p, \quad (24) $$

and for the two-layer rheologic interface

$$ \int_{I_n} \delta v_i (t_E - t_L) dt = 0 \quad \forall \delta v_i, \quad (25) $$

$$ \int_{I_n} \delta t_L \frac{1}{\beta} \dot{t}_L dt - \int_{I_n} \delta t_L v_i dt + \int_{I_n} \delta t_L v_S dt + \delta t_L(t_n^+) \frac{1}{\beta} (t_L(t_n^+) - t_L(t_n^-)) = 0 \quad \forall \delta t_L, \quad (26) $$

$$ \int_{I_n} \delta t_E \frac{1}{\gamma} t_E dt - \int_{I_n} \delta t_E v_f dt + \int_{I_n} \delta t_E v_i dt = 0 \quad \forall \delta t_E. \quad (27) $$

## 6.3. Space-time discretization

Physical variables are discretized by linear shape functions between $t_n$ and $t_{n+1}$ denoted by $t$ superscript. Space in the fluid domain is discretized by quadratic shape functions ($q$ superscript) for velocity and by linear shape functions ($l$ superscript) for pressure, the LBB condition is thus
satisfied and stabilization terms are not needed. This discretization scheme leads to the following shape functions and element matrices

$$N_1^t(t) = \frac{t_{n+1} - t}{\Delta t}; \quad N_2^t(t) = \frac{t - t_n}{\Delta t} \qquad (28)$$

$$N_1^x(x) = 1 - \frac{x}{l}; \quad N_2^x(x) = \frac{x}{l} \qquad (29)$$

$$N_1^q(x) = \frac{1}{l^2}(2x-l)(x-l); \quad N_2^q(x) = \frac{1}{l^2}4x(l-x); \quad N_3^q(x) = \frac{1}{l^2}x(2x-l) \qquad (30)$$

$$\mathbf{N}_t = [N_1^t \ N_2^t]; \quad \mathbf{N}_l = [N_1^x \ N_2^x]; \quad \mathbf{N}_q = [N_1^q \ N_2^q \ N_3^q] \qquad (31)$$

$$\mathbf{N}_p = [N_1^t \mathbf{N}_l \ N_2^t \mathbf{N}_l]; \quad \mathbf{N}_v = [N_1^t \mathbf{N}_q \ N_2^t \mathbf{N}_q] \qquad (32)$$

$$\mathbf{A} = \int_{t_n}^{t_{n+1}} \mathbf{N}_t^T \frac{\partial}{\partial t} \mathbf{N}_t \, dt; \quad \mathbf{B} = \int_{t_n}^{t_{n+1}} \mathbf{N}_t^T \mathbf{N}_t \, dt; \quad \mathbf{C} = \int_{t_n}^{t_{n+1}} \int_0^l \mathbf{N}_v^T \frac{\partial}{\partial t} \mathbf{N}_v \, dx \, dt \qquad (33)$$

$$\mathbf{D} = \int_{t_n}^{t_{n+1}} \int_0^l \frac{\partial}{\partial x} \mathbf{N}_v^T \mathbf{N}_p \, dx \, dt; \quad \mathbf{E} = \int_0^l \mathbf{N}_q^T \mathbf{N}_q \, dx; \quad \mathbf{F} = \int_{t_n}^{t_{n+1}} \int_0^l \mathbf{N}_v^T \mathbf{N}_v \, dx \, dt \qquad (34)$$

with the vectors of unknown coefficients of the approximation in space and time

$$\hat{\mathbf{v}}_i = [v_i^n \ v_i^{n+1}]^T; \quad \hat{\mathbf{t}}_E = [t_E^n \ t_E^{n+1}]^T; \quad \hat{\mathbf{t}}_L = [t_L^n \ t_L^{n+1}]^T \qquad (35)$$

$$\hat{\mathbf{v}}_S = [v_S^n \ v_S^{n+1}]^T; \quad \hat{\mathbf{f}} = [f^n \ f^{n+1}]^T \qquad (36)$$

$$\hat{\mathbf{v}}_F = [v_{F,0}^n \ v_{F,\frac{1}{2}}^n \ v_{F,l}^n \ v_{F,0}^{n+1} \ v_{F,\frac{1}{2}}^{n+1} \ v_{F,l}^{n+1}]^T \qquad (37)$$

$$\hat{\mathbf{p}} = [p_0^n \ p_l^n \ p_0^{n+1} \ p_l^{n+1}]^T \qquad (38)$$

where the superscripts n and n+1 are relative to time t_n and t_{n+1}, respectively, and the subscripts 0, l/2, and l are used to define the fluid velocity at x=0, x=l/2, and x=l, respectively. The approximation of the different variables of the problem can be summarized as

$$v_i = \mathbf{N}_i \hat{\mathbf{v}}_i; \quad t_E = \mathbf{N}_i \hat{\mathbf{t}}_E; \quad t_L = \mathbf{N}_i \hat{\mathbf{t}}_L \qquad (39)$$

$$v_S = \mathbf{N}_i \hat{\mathbf{v}}_S; \quad f = \mathbf{N}_i \hat{\mathbf{f}} \qquad (40)$$

$$v_F = \mathbf{N}_v \hat{\mathbf{v}}_F; \quad p = \mathbf{N}_p \hat{\mathbf{p}}. \qquad (41)$$

The discretized weak form for the interface is

$$\delta \hat{\mathbf{v}}_i^T \mathbf{B} \hat{\mathbf{t}}_E - \delta \hat{\mathbf{v}}_i^T \mathbf{B} \hat{\mathbf{t}}_L = 0 \quad \forall \delta \hat{\mathbf{v}}_i, \qquad (42)$$

$$\delta \hat{\mathbf{t}}_L^T \frac{1}{\beta} \mathbf{A} \hat{\mathbf{t}}_L - \delta \hat{\mathbf{t}}_L^T \mathbf{B} \hat{\mathbf{v}}_i + \delta \hat{\mathbf{t}}_L^T \mathbf{B} \hat{\mathbf{v}}_S + \delta t_L^n \frac{1}{\beta} t_L^n = \delta t_L^n t_L^{n-} \quad \forall \delta \hat{\mathbf{t}}_L, \qquad (43)$$

$$\delta \hat{\mathbf{t}}_E^T \frac{1}{\gamma} \mathbf{B} \hat{\mathbf{t}}_E - \delta \hat{\mathbf{t}}_E^T \mathbf{B} \hat{\mathbf{v}}_F + \delta \hat{\mathbf{t}}_E^T \mathbf{B} \hat{\mathbf{v}}_i = 0 \quad \forall \delta \hat{\mathbf{t}}_E. \qquad (44)$$

The discretized weak form for the solid is

$$\delta \hat{\mathbf{v}}_S^T m \mathbf{A} \hat{\mathbf{v}}_S + \delta \hat{\mathbf{v}}_S^T \mathbf{B} \hat{\mathbf{f}} - \delta \hat{\mathbf{v}}_S^T \mathbf{B} \hat{\mathbf{t}}_L + \delta v_S^n m v_S^n = \delta v_S^n m v_S^{n-} \quad \forall \delta \hat{\mathbf{v}}_S, \qquad (45)$$

$$\delta \hat{\mathbf{f}}^T \frac{1}{k} \mathbf{A} \hat{\mathbf{f}} - \delta \hat{\mathbf{f}}^T \mathbf{B} \hat{\mathbf{v}}_S + \delta f^n \frac{1}{k} f^n = \delta f^n \frac{1}{k} f^{n-} \quad \forall \delta f. \qquad (46)$$
<table><caption>Table III. Parameters for the one-dimensional case.</caption><thead><tr><th>β</th><th>γ</th><th>k</th><th>m</th><th>ρ</th><th>l</th><th>Ω</th><th>τ̃<sub>0</sub></th></tr></thead><tbody><tr><td>100</td><td>0.01</td><td>500</td><td>50</td><td>1000</td><td>1</td><td>1.2</td><td>20 000</td></tr></tbody></table>

The discretized weak form for the fluid is

$$ \delta \hat{v}_F^T \rho C \hat{v}_F - \delta \hat{v}_F^T D \hat{p} \qquad (47) $$

$$ -[\delta v_{F,0}^n \ \delta v_{F,0}^{n+1}] \mathbf{B} \hat{t}_E + [\delta v_{F,0}^n \delta v_{F,\frac{l}{2}}^n \ \delta v_{F,l}^n] \rho \mathbf{E} [v_{F,0}^n \ v_{F,\frac{l}{2}}^n \ v_{F,l}^n]^T \qquad (48) $$

$$ = [\delta v_{F,l}^n \ \delta v_{F,l}^{n+1}] \hat{\mathbf{B}} \hat{t} + [\delta v_{F,0}^n \ \delta v_{F,\frac{l}{2}}^n \ \delta v_{F,l}^n] \rho \mathbf{E} [v_{F,0}^{n-} \ v_{F,\frac{l}{2}}^{n-} \ v_{F,l}^{n-}]^T \quad \forall \delta \hat{v}_F, \qquad (49) $$

$$ \delta \hat{\mathbf{p}}^T D^T \hat{v}_F = 0 \quad \forall \delta \hat{\mathbf{p}}. \qquad (50) $$

The final discretized system of Equations (42)–(50) consists of a total of 20 unknowns given by the vectors $\hat{v}_i$, $\hat{t}_E$, $\hat{t}_L$, $\hat{v}_S$, $\hat{\mathbf{f}}$, $\hat{v}_F$, and $\hat{\mathbf{p}}$. This linear symmetric system is solved for each time slab.

## 6.4. Verification of the method by Newmark reference solution

The pressure and the velocity in the fluid domain are solved analytically writing the momentum balance and continuity equations

$$ -p_{,x} = \rho \dot{v}_F \qquad (51) $$

$$ v_{F,x} = 0 \qquad (52) $$

with $p(x=l, t)=-\bar{l}(t)$, $p(x, t=0)=0$, and $v_F(x, t=0)=0$. This leads to the solution

$$ p(x,t) = \rho \dot{v}_F(l-x) - \bar{l}(t), \qquad (53) $$

where $v_F$ is only a function of time.

The momentum balance applied successively to the solid, to the interface and the fluid connection to the Eulerian interface layer lead to the following time dependant discrete system:

$$ \begin{bmatrix} m & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \rho l \end{bmatrix} \begin{bmatrix} \dot{v}_S \\ \dot{v}_i \\ \dot{v}_F \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 \\ 0 & \gamma & -\gamma \\ 0 & -\gamma & \gamma \end{bmatrix} \begin{bmatrix} v_S \\ v_i \\ v_F \end{bmatrix} + \begin{bmatrix} k+\beta & -\beta & 0 \\ -\beta & \beta & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} x_S \\ x_i \\ x_F \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \bar{l}(t) \end{bmatrix}, \qquad (54) $$

where $v_F$ denotes the fluid velocity at the interface, $x_S$, $x_i$, and $x_F$ are respectively the solid, interface, and fluid positions. In the last equation the fluid acts as an added mass to the Eulerian interface layer. This system is solved by the Newmark method with the standard value $\beta_{\text{Newmark}} = 0.25$ and $\gamma_{\text{Newmark}} = 0.5$, leading to a second-order implicit scheme in conjunction with a time step $\Delta t = 0.01$ to produce a qualified reference solution.

The parameters for verification of the space-time approach are given in Table III. The space-time computation is performed using the same time step $\Delta t = 0.01$ as the Newmark-based reference computation. The results for a given traction $\bar{l}(t) = \bar{l}_0 H(t)$ are shown in Figures 10(a), (c), and (e). and those for the prescribed harmonic traction $\bar{l}(t) = \bar{l}_0 \sin(\Omega t)$ in Figures 10(b), (d), and (f). The results obtained by the proposed approach are in very good agreement with the reference solution.
<figure><img src="image_17.png" /></figure>

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

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

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

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

<figure><img src="image_22.png" /><figcaption>Figure 10. Results of the verification for the one-dimensional problem: (a) fluid velocity $v_f(0)$ at the interface, $\bar{t}(t)=\bar{t}_0 H(t)$; (b) fluid velocity $v_f(0)$ at the interface, $\bar{t}(t)=\bar{t}_0 \sin(\Omega t)$; (c) interface velocity $v_i$, $\bar{t}(t)=\bar{t}_0 H(t)$; (d) interface velocity $v_i$, $\bar{t}(t)=\bar{t}_0 \sin(\Omega t)$; (e) structure velocity $v_S$, $\bar{t}(t)=\bar{t}_0 H(t)$; and (f) structure velocity $v_S$, $\bar{t}(t)=\bar{t}_0 \sin(\Omega t)$.</figcaption></figure>

It is expected that the time-discontinuous Galerkin method will outperform the Newmark method if discontinuous loading functions are applied to the system.

## 6.5. Investigation on the influence of Eulerian interface damping

The eigenfrequency of the spring mass system is given by $\omega_s = \sqrt{k/m} = 3.16$. In order to study the influence of the damping interface parameter $\gamma$, the load is chosen as $\bar{t}(t)=\bar{t}_0 \sin(\Omega t)$ with
<figure><img src="image_23.png" /><figcaption>Figure 11. Influence of Eulerian interface damping for the one-dimensional problem.</figcaption></figure>

$\Omega=0.95\omega_s$ so that the excitation frequency is very close to the eigenfrequency of the structural system. The Eulerian interface damping parameter takes the values $\gamma=1, 0.1$, and $0.01$ while the other numerical values for this test are the same as in the previous verification.

The structure velocity is plotted for the three values of $\gamma$ in Figure 11. A large value of damping ($\gamma=1$) is close to a perfect interface, meaning that energy is transferred correctly from the fluid to the structure while a small value ($\gamma=0.01$) is close to a slip boundary condition. The results show that smaller damping values lead to reduced oscillation amplitudes, which is in good agreement with the expectations.

## 7. APPLICATION TO VORTEX-INDUCED VIBRATIONS

In this section the two-layer interface model is applied to the two-dimensional problem of flow-induced vibration of a spring-supported cylinder in a cross flow (Figure 12). This test case reveals two interesting and examinable phenomena: von Kármán vortex streets and the lock-in effect. The vortices are alternately detaching from both sides of the body and then convected with the flow. Coupled to the periodic vortex shedding, periodic pressure fluctuations occur on the surface of the body. The alternation leads to a harmonic excitation of the body, hence the spring-supported cylinder may start moving orthogonal to the flow direction. If the vortex shedding frequency $f$ is close to the eigenfrequency $f_0$ of the structural system, resonance and vibrations at large amplitudes may develop. In certain ranges of Reynolds number a reverse interaction from solid to fluid can be observed. The frequency of eddy detachment then adapts to the eigenfrequency of the structure, called the *lock-in effect*.

The diameter of the cylinder is $d=0.16$ cm and the inflow on the left boundary is denoted by $v_x=v_\infty$. At the outlet, the pressure gradient is set to zero, i.e. the pressure is set to the surrounding reference pressure $p=0$. The cylinder is spring-supported in the transverse flow direction. The stiffness of the spring is $k=6793g \cdot s^{-2}$ and the density of the cylinder is given by $\rho_s=148.2g \cdot cm^{-3}$. This parameter setup determines the eigenfrequency of the structural system to $f=7.6$ Hz. The material parameters of water at a temperature of $T=20°C$ are used for the fluid: $\mu=0.01002g \cdot s^{-1} \cdot cm^{-1}$ and $\rho=0.9982g \cdot cm^{-3}$. This setup has also been used by Walhorn [33], so that the simulations with rheological devices can be compared to simulations with strong coupling conditions. Moreover, there are experimental examinations for this case, coupling the detaching frequency to the Reynolds number of the flow [34].
<figure><img src="image_24.png" /><figcaption>Figure 12. Spring-supported rigid cylinder in a flow.</figcaption></figure>

<table><caption>Table IV. Inflow velocities for different Reynolds numbers.</caption><thead><tr><th>Re</th><th>90</th><th>100</th><th>110</th><th>120</th><th>130</th></tr></thead><tbody><tr><th>v<sub>∞</sub>[cm·s<sup>-1</sup>]</th><td>5.6464</td><td>6.2738</td><td>6.9012</td><td>7.5286</td><td>8.1559</td></tr></tbody></table>

To trace the eddy detachment frequency, the horizontal fluid velocity at the monitoring point A is observed. The state at the interface is observed by tracking displacements and velocities at the monitoring points B and C. The flow is disturbed by an asymmetric inflow velocity during the first time steps in order to accelerate the development of vortices.

## 7.1. Fixed cylinder without rheological interface

In this first test the cylinder is fixed ($v_s=0$) and the rheological interface model is not active. In order to verify the fluid model, the behavior at different Reynolds numbers ($Re=v_∞·d·ρ/μ$) is investigated by changing the inflow velocity according to Table IV. The developing velocity and pressure field is depicted in Figures 13(a) and (b). The von Kármán vortex street can be clearly resolved, the transportation of the eddies with the fluid velocity is observed. The typical asymmetric and time-dependent stress distribution is present on the surface of the cylinder. The resulting Strouhal number ($Sr=f·d/v_∞$) in comparison with the results of Walhorn [33], Sarrate and Huerta [35] and the relation stated by Roshko [34] is depicted in Figure 13(c). The Strouhal number obtained by simulations on the chosen mesh is slightly above the results of Sarrate and Roshko. The results show that the considered mesh can be used in order to further examine the qualitative behavior of the rheological interface.

## 7.2. Fixed cylinder with active rheological interface

In the following, the test case described above is examined for interface layers with finite stiffness and damping. The impact of rheological interfaces on the development of the von Kármán vortex street is examined. The solid body is still fixed in space while the rheological interface is now active. Consequently, interface deformations are possible and lead to different shear stress distributions at the interface. The vortex shedding frequency $f_{Re=110}≈7.6$ Hz, associated to $Re=110$ without rheological interface effects, matches the eigenfrequency $f_s=7.6$ Hz of the cylindrical rigid body system. The effect of the two components of the rheological interface model on the vortex shedding frequency at the specific Reynolds number $Re=110$ is therefore investigated in this section.

### 7.2.1. Lagrangian part of interface layer.
The Kelvin–Voigt model is applied exemplarily for the Lagrangian layer of the interface. For decreasing damping $γ$ this model approaches the pure Hooke device. At constant stiffness $β=10^4$, the damping parameter is varied in the range
<figure><img src="image_25.png" /><figcaption>Figure 13. Fixed cylinder in a flow: (a) velocities and pressure field with shedding vortices (fine mesh); (b) velocities and pressure field around the cylinder; and (c) comparison of Strouhal number.</figcaption></figure>

$\gamma \in [2 \times 10^2, 6 \times 10^5]$. For increasing damping the eddy shedding frequency converges to the frequency of the direct coupling case, see Figure 14(a). The Lagrangian part of the rheological interface leads to a variation of only 1% in the vortex shedding frequency. At constant damping $\gamma=10^3$, the stiffness parameter is varied in the range $\beta \in [5 \times 10^3, 1 \times 10^5]$. The resulting vortex shedding frequency band is between 7.3 and 7.6 Hz as shown in Figure 14(b). According to theory, direct coupling conditions are implemented for large values of $\beta$. Normal and tangential motions are possible in the rheological interface with active Lagrangian layer. The interface deformation changes the boundary of the flow domain and therefore influences the vortex shedding frequency.

### 7.2.2. Eulerian part of interface layer.
A Newtonian rheological interface is applied in the tangential direction while in normal direction kinematic consistency is realized. In Figure 15(a), the resulting horizontal velocity in control point C is depicted. For large values of tangential damping parameter $\gamma_t$ the resulting velocity converges to zero, corresponding to a no-slip interface condition (direct kinematic coupling of fluid and structure velocities). For increasing damping coefficients, a linear relation between damping and tangential velocity is observed. In this region the impact of pressure and shear stress from the fluid on the cylinder is nearly constant. Below $\gamma \approx 1$ the behavior changes and the tangential velocity converges to the free-stream velocity of the surrounding flow with slip interface conditions ($\gamma \to 0$). The impact on the frequency of detaching eddies is shown in Figure 15(b). For a damping coefficient larger than $\gamma=10$ the vortex shedding frequency is not changed. For $\gamma \in [1, 10]$ the frequency increases by 7%. In a narrow transition region between $\gamma=1$
<figure><img src="image_26.png" /><figcaption>Figure 14. Influence of the Lagrangian part of the rheological interface on the vortex shedding frequency: (a) results for constant stiffness and (b) results for constant damping.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 15. Eulerian part of rheological interface with Newton device in tangential direction only: (a) horizontal velocity in point C and (b) eddy frequency.</figcaption></figure>

and $\gamma=0.5$ vortex shedding does not take place anymore. This transition region is very small and the eddy detaching behavior is very sensitive.

### 7.3. Spring-supported cylinder with rheologic interface at $Re=110$

The rigid cylinder is now free to move and the vortex shedding frequency associated with $Re=110$ is again close to the eigenfrequency of the spring mass system such that flow-induced vibrations can be expected. To check the influence of a rheological interface on the resonance behavior, the investigation of the tangential Newtonian interface is repeated for the full fluid-structure system.

The displacement of the cylinder is plotted for various damping parameters in Figure 16. For large damping parameters ($\gamma>5$) the tangential Newtonian device meet no slip interface conditions for the fluid-structure interaction system. Therefore energy is constantly transferred from the flow to the rigid cylinder due to active shear stress and resonance can be observed. If the damping parameter $\gamma$ is decreased ($\gamma=1, 2$) a transition to slip interface conditions is realized and the eddy detachment frequency increases up to 15%. This shift between the eigenfrequency of the structure and the frequency of the exciting flow force now avoids the resonance case and leads to limit cycle oscillations of the cylinder at still large amplitudes. For $\gamma<0.5$ the eddies vanish completely as already seen for the rigid cylinder in the previous section. The result is a full decoupling of flow shear stresses and structural forces since the pressure field is symmetric.
<figure><img src="image_28.png" /><figcaption>Figure 16. Influence of rheological damping on resonance behavior.</figcaption></figure>

This investigation shows that already a very simple rheological device is able to influence the
behavior of a fluid-structure interaction system.

8. CONCLUSIONS

A rheological interface model for fluid-structure interaction problems is introduced in order to investigate its influence on flow-induced vibrations of the coupled system. The interface model is composed of two characteristic layers, a Lagrangian for the solid-like part and a Eulerian for the fluid-like part. The formulation introduces interface velocity and interface tractions for both layers in addition to fluid and solid variables. In this work the fluid-structure system is assumed to consist of a spring-supported rigid body and an incompressible viscous fluid. The velocity-based weak form of the equations for the whole coupled system including the rheological interface is discretized using space-time finite elements. This approach leads to a monolithic and uniform finite element description of the coupled system including the proposed rheologic interface model.

The interface model is first validated and then applied to a simplified one-dimensional generic problem. It is shown that by introduction of a proper damping interface, oscillations of the coupled system are reduced when it is harmonically excited close to the eigenfrequency of the structural part.

The influence of various combinations of rheologic interface models to flow-induced vibrations of a cylinder is studied. At the critical Reynolds number, parameter studies of the rheologic model are performed. It is shown that a rheologic interface is able to change the vortex shedding frequency and therefore can modify the sensitivity to resonance behavior of the flow-immersed cylindrical rigid body.

An important aspect of the developed model is the determination of correct model parameters within the two different interface layers. These parameters can be obtained by either experiments or numerical analysis of the real multidimensional interface setup and dimensional reduction. Moreover, it is known that for vibro-acoustic problems interface model parameters are frequency-dependent, hence it is expectable to observe dependency on the Reynolds number also for free flow fluid-structure interaction problems. This would lead to a non-linear and time-dependent rheological interface model, which is a topic for further investigation.

ACKNOWLEDGEMENTS

This publication has been realized within a cooperation between the Institut für Statik at the Technische Universität Braunschweig and the Structural Mechanics and Coupled Systems Laboratory of the Conser-vatoire National des Arts et Métiers (CNAM) in Paris. The collaboration is supported by the PROCOPE programme between France (PAI-EGIDE) and Germany (DAAD).
1. Allard J-F. *Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials*. Elsevier Applied Science: London, New York, 1993.

2. Atalla N, Panneton R, Debergue P. A mixed displacement-pressure formulation for poroelastic materials. *Journal of Acoustical Society of America* 1998; **104**(3):1444-1452.

3. Goransson P. A 3-D, symmetric, finite element formulation of the Biot equations with application to acoustic wave propagation through an elastic porous medium. *International Journal for Numerical Methods in Engineering* 1998; **41**(1):167-192.

4. Davidsson P, Sandberg G. A reduction method for Structure-acoustic and poroelastic-acoustic problems using interface-dependent Lanczos vectors. *Computer Methods in Applied Mechanics and Engineering* 2006; **195** (17-18):1933-1945.

5. Larbi W, Deü J-F, Ohayon R. A new finite element formulation for internal acoustic problems with dissipative walls. *International Journal for Numerical Methods in Engineering* 2006; **68**(3):381-399.

6. Deü J-F, Larbi W, Ohayon R. Vibration and transient response of structural acoustic interior coupled systems with dissipative interface. *Computer Methods in Applied Mechanics and Engineering* 2008; **197**(51-52):4894-4905.

7. Bermúdez A, Rodríguez R. Modelling and numerical solution of elastoacoustic vibrations with interface damping. *International Journal for Numerical Methods in Engineering* 1999; **46**(10):1763-1779.

8. Cebral JR. Loose coupling algorithms for fluid-structure interaction. Ph.D. Thesis, George Mason University, Fairfax, Virginia, 1996.

9. Farhat C, Lesoinne M. Two efficient staggered procedures for the serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems. *Computer Methods in Applied Mechanics and Engineering* 2000; **182**:499-515.

10. Löhner R, Yang C, Cebral J, Baum JD, Luo H, Pelessone D. Fluid-structure interaction using a loose coupling algorithm and adaptive unstructured grids. *Computational Fluid Dynamics Review*. Wiley: New York, 1995; 755-776.

11. Wall WA, Mok DP, Ramm E. Partitioned analysis approach for the transient coupled response of viscous fluids and flexible structures. In *Proceedings of the European Conference on Computational Mechanics*, ECCM '99, Wunderlich W (ed.), Munich, 1999.

12. Le Tallec P, Mouro J. Fluid-structure interaction with large structural displacements. *Computer Methods in Applied Mechanics and Engineering* 2001; **190**:3039-3067.

13. Steindorf J. Partitionierte Verfahren für Probleme der Fluid-Struktur-Wechselwirkung. *Dissertation*, Technische Universität Braunschweig, 2002.

14. Mok D, Wall WA. Partitioned analysis schemes for the transient interaction of incompressible flows and nonlinear flexible structures. In *Trends in Computational Structural Mechanics*, Wall WA, Bletzinger K-U, Schweizerhof K (eds). CIMNE: Barcelona, 2001.

15. Wüchner R, Kupzok A, Bletzinger K-U. A framework for stabilized partitioned analysis of thin membrane-wind interaction. *International Journal for Numerical Methods in Fluids* 2007; **54**(6-8):945-963.

16. Michler C, van Brummelen E, de Borst R. An interface Newton-Krylov solver for fluid-structure interaction. *International Journal for Numerical Methods in Fluids* 2005; **47**(10-11):1189-1195.

17. Fernandez M, Moubachir M. A Newton method using exact Jacobians for solving fluid-structure coupling. *Computers and Structures* 2005; **83**:127-142.

18. Tezduyar TE, Sathe S, Stein K. Solution techniques for the fully discretized equations in computation of fluid-structure interactions with the space-time formulations. *Computer Methods in Applied Mechanics and Engineering* 2006; **195**:5743-5753.

19. Dettmer WG, Peric D. A computational framework for fluid-structure interaction: finite element formulation and applications. *Computer Methods in Applied Mechanics and Engineering* 2006; **195**:5754-5779.

20. Matthies HG, Niekamp R, Steindorf J. Algorithms for strong coupling procedures. *Computer Methods in Applied Mechanics and Engineering* 2006; **195**:2028-2049.

21. Küttler U, Wall W. Fixed-point fluid-structure interaction solvers with dynamic relaxation. *Computational Mechanics* 2008; **43**(1):61-72.

22. Grohmann B, Dornberger R, Dinkler D. Time-discontinuous stabilized space time finite elements for aeroelasticity. *Proceedings of the Fourth International Symposium on Fluid-Structure Interactions, Aeroelasticity, Flow-induced Vibration and Noise*, Dallas, U.S.A., 1997.

23. Hübner B, Dinkler D. Numerische Analyse aeroelastischer Instabilitäten von Brückenquerschnitten. In *Windeinwirkungen auf Bauwerke*, Udo Peil (ed.). Windtechnologische Gesellschaft: Aachen, 2001.

24. Tezduyar TE, Behr M, Liou J. A new strategy for finite element computations involving moving boundaries and interfaces—the deforming-spatial-domain/space-time procedure: I. The concept and the preliminary numerical tests. *Computer Methods in Applied Mechanics and Engineering* 1992; **94**(3):339-351.

25. Hughes TJR, Hulbert GM. Space-time finite element methods for elastodynamics: formulations and error estimates. *Computer Methods in Applied Mechanics and Engineering* 1988; **66**(3):339-363.

26. Hübner B, Walhorn E, Dinkler D. A monolithic approach to fluid-structure interaction using space-time finite elements. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**(23-26):2069-2086.

27. Walhorn E, Kölke A, Hübner B, Dinkler D. Fluid-structure coupling within a monolithic model involving free surface flows. *Computers and Structures* 2005; **83**(25-26):2100-2111.
28. Zilian A, Legay A. The enriched space-time finite element method (EST) for simultaneous solution of fluid-structure interaction. *International Journal for Numerical Methods in Engineering* 2008; **75**(3):305-334.

29. Zilian A, Dinkler D, Vehre A. Projection-based reduction of fluid-structure interaction systems using monolithic space-time modes. *Computer Methods in Applied Mechanics and Engineering* 2009; **198**(47-48):3795-3805.

30. Grau P, Berg G, Giessmann E-J. Rheologische Untersuchungen fester Stoffe mit Dehnungsratenwechselversuchen. *Technische Mechanik* 1983; **4**:3-54.

31. Pierce AD. Acoustics: An Introduction to its Physical Principles and Applications. McGraw-Hill Education: New York, 1989.

32. Masud A, Hughes TJR. A space-time Galerkin/least-squares finite element formulation of the Navier-Stokes equations for moving domain problems. *Computer Methods in Applied Mechanics and Engineering* 1997; **146**(1-2):91-126.

33. Walhorn E. Ein simultanes Berechnungsverfahren für Fluid-Struktur-Wechselwirkungen mit finiten Raum-Zeit-Elementen. Ph.D. Thesis, Institut für Structural Analysis, TU Braunschweig, Braunschweig, 2002.

34. Roshko A. *On the Development of Turbulent Wakes from Vortex Streets*. California Institute of Technology: Pasadena, CA, 1953.

35. Sarrate J, Huerta A. ALE formulation for wake flows accepting secondary objects. *First European Conference on Computational Mechanics: Solids, Structures, and Coupled Problems in Engineering*, München, Germany, 1999.