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

HAL
open science

# Embedding Domain Technique for a Fluid-Structure Interaction Problem

Cornel Marius Murea, Andrei Halanay

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

## To cite this version:

Cornel Marius Murea, Andrei Halanay. Embedding Domain Technique for a Fluid-Structure Interaction Problem. 25th System Modeling and Optimization (CSMO), Sep 2011, Berlin, Germany. pp.358-367, 10.1007/978-3-642-36062-6_36. hal-01347557

HAL Id: hal-01347557

https://inria.hal.science/hal-01347557

Submitted on 21 Jul 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Embedding domain technique for a
fluid-structure interaction problem

Cornel Marius Murea¹ and Andrei Halanay²

¹ Laboratoire de Mathématiques, Informatique et Applications,
Université de Haute Alsace, France
cornel.murea@uha.fr
http://www.edp.lmia.uha.fr/murea/

² Department of Mathematics 1,
University Politehnica of Bucharest, Romania,
halanay@mathem.pub.ro

**Abstract.** We present a weak formulation for a steady fluid-structure interaction problem using an embedding domain technique with penalization. Except of the penalizing term, the coefficients of the fluid problem are constant and independent of the deformation of the structure, which represents an advantage of this approach. A second advantage of this model is the fact that the continuity of the stress at the fluid-structure interface does not appear explicitly. Numerical results are presented.

**Keywords:** fluid-structure interaction; embedding domain;

# 1 A steady fluid-structure interaction problem

The present paper is devoted to the study of the numerical behavior of an elastic structure immersed in a viscous incompressible fluid. We use Stokes equation to model the flow motion. The displacement of the structure under the flow motion will be modeled by linear elasticity equations, under the small deformations assumption. In this paper, we study the steady case.

Let $D \subset \mathbb{R}^2$ be a bounded open domain with boundary $\partial D$. Let $\Omega_0^S$ be the undeformed structure domain, and suppose that its boundary admits the decomposition $\partial \Omega_0^S = \Gamma_D \cup \Gamma_0$, where $\Gamma_0$ is a relatively open subset of the boundary. On $\Gamma_D$ we impose zero displacement for the structure. We assume that $\Omega_0^S \subset D$.

Suppose that the structure is elastic and denote by $\mathbf{u} = (u_1, u_2) : \Omega_0^S \to \mathbb{R}^2$ its displacement. A particle of the structure with initial position at the point $\mathbf{X}$ will occupy the position $\mathbf{x} = \varphi(\mathbf{X}) = \mathbf{X} + \mathbf{u}(\mathbf{X})$ in the deformed domain $\Omega_u^S = \varphi(\Omega_0^S)$.

We assume that $\Omega_u^S \subset D$ and the fluid occupies $\Omega_u^F = D \setminus \bar{\Omega}_u^S$. We set $\Gamma_u = \varphi(\Gamma_0)$ and we suppose that $\Gamma_u$ does not touch the container wall, i.e. $\partial D \cap \Gamma_u = \emptyset$. We recall that $\Gamma_0$ is a relatively open subset. The boundary $\Gamma_u$ represents the moving fluid-structure interface. The boundary of the deformed
structure is $\partial\Omega_u^S = \Gamma_D \cup \Gamma_u$. In the case when $\bar{\Omega}_u^S \subset D$, the fluid-structure geometrical configuration is represented in Figure 1 and the boundary of the fluid domain admits the decomposition $\partial\Omega_u^F = \partial D \cup \Gamma_D \cup \Gamma_u$.

<figure><img src="image_4.png" /><figcaption>Fig. 1. Geometrical configuration.</figcaption></figure>

Also, for the first numerical test, we consider the case when $\Gamma_D \subset \partial D$ as in Figure 2 and the boundary of the fluid domain admits the decomposition $\partial\Omega_u^F = (\partial D \setminus \Gamma_D) \cup \Gamma_u$.

## 2 Weak formulation using an embedding domain technique with penalization

We introduce the tensor $\epsilon(\mathbf{w}) = \frac{1}{2}(\nabla\mathbf{w} + (\nabla\mathbf{w})^T)$ and we assume that the fluid is Newtonian and the Cauchy stress tensor is given by $\sigma^F(\mathbf{v}, p) = -p\mathbf{I} + 2\mu^F\epsilon(\mathbf{v})$, where $\mu^F > 0$ is the viscosity of the fluid and $\mathbf{I}$ is the unit matrix. We assume that the structure verifies the linear elasticity equation, under the assumption of small deformations. The stress tensor of the structure written in the Lagrangian framework is $\sigma^S(\mathbf{u}) = \lambda^S(\nabla \cdot \mathbf{u})\mathbf{I} + 2\mu^S\epsilon(\mathbf{u})$, where $\lambda^S, \mu^S > 0$ are the Lamé coefficients.

We present in an informal and intuitive manner the ideas behind our approximation approach using embedding domain technique with penalization. In the fluid domain, Stokes equations are solved:

$$-\nabla \cdot \sigma^F(\mathbf{v}, p) = \mathbf{f}^F, \quad \text{in } \Omega_u^F \qquad (1)$$

$$\nabla \cdot \mathbf{v} = 0, \quad \text{in } \Omega_u^F \qquad (2)$$

We introduce two more equations concerning the fluid fields, but written on the deformed structure domain:

$$-\nabla \cdot \sigma^F(\mathbf{v}, p) + \frac{1}{\varepsilon}\mathcal{P}(\mathbf{v}) = \mathbf{f}^F, \quad \text{in } \Omega_u^S \qquad (3)$$

$$\nabla \cdot \mathbf{v} = 0, \quad \text{in } \Omega_u^S \qquad (4)$$
where $\varepsilon > 0$ is a penalization parameter,

$$ \mathcal{P}(\mathbf{v}) = (|v_1|^{\alpha-1} \operatorname{sgn}(v_1), |v_2|^{\alpha-1} \operatorname{sgn}(v_2)) \quad (5) $$

where $\mathbf{v} = (v_1, v_2)$ and $1 < \alpha < 2$ is a real number.

**Remark 1.** This choice of the penalization term is justified in [2], in order to obtain existence of the fluid-structure interaction problem. For the steady case, the role of the penalization term is to obtain very small values of the fluid velocity in the structure domain. If we take other values for $\alpha$ or for a $H^1$ penalization, we do not get more regularity of the solution.

Let $\chi_u^S$ be the characteristic function of $\Omega_u^S$. Combining (1) and (3), it follows that

$$ -\nabla \cdot \sigma^F(\mathbf{v}, p) + \frac{1}{\varepsilon} \chi_u^S \mathcal{P}(\mathbf{v}) = \mathbf{f}^F, \quad \text{in } D. \quad (6) $$

Similarly, we have from (2) and (4)

$$ \nabla \cdot \mathbf{v} = 0, \quad \text{in } D. \quad (7) $$

In view of the equation (3), the “fictitious” fluid velocity and pressure defined on the structure domain $\Omega_u^S$ depend on $\varepsilon$. In the following, we denote by $\mathbf{v}_\varepsilon$ and $p_\varepsilon$ the fluid velocity and pressure defined all over the domain $D$.

Let us introduce the bi-linear forms

$$ a_S(\mathbf{u}, \mathbf{w}^S) = \int_{\Omega_0^S} (\lambda^S(\nabla \cdot \mathbf{u})(\nabla \cdot \mathbf{w}^S) + 2\mu^S \varepsilon(\mathbf{u}) : \varepsilon(\mathbf{w}^S)) \, d\mathbf{X} $$

$$ a_F(\mathbf{v}, \mathbf{w}) = \int_D 2\mu^F \varepsilon(\mathbf{v}) : \varepsilon(\mathbf{w}) \, d\mathbf{x} $$

$$ b_F(\mathbf{w}, p) = - \int_D (\nabla \cdot \mathbf{w}) \, p \, d\mathbf{x} $$

and the Hilbert spaces

$$ W^S = \{ \mathbf{w}^S \in (H^1(\Omega_0^S))^2; \mathbf{w}^S = 0 \text{ on } \Gamma_D \}, $$

$$ W = (H_0^1(D))^2, $$

$$ Q = L_0^2(D) = \{q \in L^2(D); \int_D q \, dx = 0\}. $$

We assume for the moment that $\mathbf{f}^F \in (L^2(D))^2$, $\mathbf{f}^S \in (L^2(\Omega_0^S))^2$ and $\mathbf{g} \in (H^{1/2}(\partial D))^2$, such that $\int_{\partial D} \mathbf{g} \cdot \mathbf{n}^F ds = 0$.

For a given $\mathbf{u} \in (W^{1,\infty}(\Omega_0^S))^2$, such that $\|\mathbf{u}\|_{1,\infty,\Omega_0^S} < 1$ and $\mathbf{u} = 0$ on $\Gamma_D$, we define:

- fluid velocity $\mathbf{v}_\varepsilon \in (H^1(D))^2$, $\mathbf{v}_\varepsilon = \mathbf{g}$ on $\partial D$,

- fluid pressure $p_\varepsilon \in Q$,
– structure displacement $\mathbf{u}_{\varepsilon} \in W^S$,

as the solution of the following weakly coupled system of PDE’s:

$$a_F(\mathbf{v}_\varepsilon, \mathbf{w}) + b_F(\mathbf{w}, p_\varepsilon) + \frac{1}{\varepsilon} \int_D \tilde{H}_u \mathcal{P}(\mathbf{v}_\varepsilon) \cdot \mathbf{w} \, d\mathbf{x} = \int_D \mathbf{f}^F \cdot \mathbf{w} \, d\mathbf{x}, \quad \forall \mathbf{w} \in W \quad (8)$$

$$b_F(\mathbf{v}_\varepsilon, q) = 0, \quad \forall q \in Q \quad (9)$$

$$a_S(\mathbf{u}_\varepsilon, \mathbf{w}^S) = \int_{\Omega_0^S} \mathbf{f}^S \cdot \mathbf{w}^S \, d\mathbf{X} + \int_{\Omega_0^S} J(\sigma^F(\mathbf{v}_\varepsilon, p_\varepsilon) \circ \varphi) \mathbf{F}^{-T} : \nabla_{\mathbf{X}} \mathbf{w}^S \, d\mathbf{X} \\ + \frac{1}{\varepsilon} \int_{\Omega_0^S} J \tilde{H}_u \mathcal{P}(\mathbf{v}_\varepsilon \circ \varphi) \cdot \mathbf{w}^S \, d\mathbf{X} \\ - \int_{\Omega_0^S} J(\mathbf{f}^F \circ \varphi) \cdot \mathbf{w}^S \, d\mathbf{X}, \quad \forall \mathbf{w}^S \in W^S \quad (10)$$

where $\varphi(\mathbf{X}) = \mathbf{X} + \mathbf{u}(\mathbf{X})$, $\mathbf{F}(\mathbf{X}) = \mathbf{I} + \nabla_{\mathbf{X}} \mathbf{u}(\mathbf{X})$, $J(\mathbf{X}) = \det \mathbf{F}(\mathbf{X})$.

The equations (8) and (9) are obtained from (6) and (7). The coefficient $\tilde{H}_u$ in (8) is a regularization of the characteristic function of $\Omega_u^S$, which is necessary in order to prove the continuity of the solution with respect to the structure displacement.

**Remark 2.** From the structure equation $-\nabla \cdot \sigma^S(\mathbf{u}_\varepsilon) = \mathbf{f}^S$, in $\Omega_0^S$ using Green's formula, we obtain for all $\mathbf{w}^S = 0$ on $\Gamma_D$ that

$$a_S(\mathbf{u}_\varepsilon, \mathbf{w}^S) = \int_{\Omega_0^S} \mathbf{f}^S \cdot \mathbf{w}^S \, d\mathbf{X} + \int_{\Gamma_0} \sigma^S(\mathbf{u}_\varepsilon) \mathbf{n}^S \cdot \mathbf{w}^S \, dS.$$

We can prove (see [2]) that the sum of the last three terms in (10) is equal to the fluid forces acting on the structure which is also equals to $\int_{\Gamma_0} \sigma^S(\mathbf{u}_\varepsilon) \mathbf{n}^S \cdot \mathbf{w}^S \, dS$. In fact, from (10) and the above weak formulation of the structure, we can get that the boundary condition at the interface concerning the continuity of the stress is verified in a weak sense (see [2]). The second boundary condition at the interface is the continuity of the velocity, i.e. $\mathbf{v} = 0$ on $\Gamma_u$ in the steady case. This is obtained by using the penalization term in the structure domain.

Define the nonlinear operator

$$T_\varepsilon : \{ \mathbf{u} \in (W^{1,\infty}(\Omega_0^S))^2 ; \|\mathbf{u}\|_{1,\infty,\Omega_0^S} < 1, \mathbf{u} = 0 \text{ on } \Gamma_D \} \to (W^{1,\infty}(\Omega_0^S))^2$$

by $T_\varepsilon(\mathbf{u}) = \mathbf{u}_\varepsilon$. A solution of the penalized fluid-structure interaction problem will be, by definition, a fixed point of $T_\varepsilon$. In [2], we discuss the existence of a solution of the penalized fluid-structure interaction problem. The convergence of $\mathbf{u}_\varepsilon, \mathbf{v}_\varepsilon, p_\varepsilon$ when $\varepsilon$ goes to 0 is also analyzed.
## 3 Partitioned procedures based on fixed point iterations

The penalized term $P(v)$ is non-linear in $v$ for $\alpha \neq 2$. But for $\alpha = 2$, we have $P(v) = v$. Now, the fluid problem, at the **Step 2** of the algorithm below, becomes linear and, for a given $u_{\epsilon}^k$, it has a unique solution.

For $1 < \alpha < 2$, we can prove the existence of a fixed point for the nonlinear operator $T_{\epsilon}$ defined at the end of the previous section, but not for $\alpha = 2$. We can also replace $\tilde{H}_u$ in (8) by $\chi_u^S$ the characteristic function of $\Omega_u^S$ in order to simplify the computation. The regularization of the characteristic function was necessary in order to prove the continuity of the solution with respect to the structure displacement.

Under the assumption of small displacements for the structure, we can approach the Jacobian determinant $J$ by 1 and the gradient of the deformation $F$ by the identity matrix $I$. The structure problem at the **Step 3** is linear and, for given $v_{\epsilon}^k$ and $p_{\epsilon}^k$, it has a unique solution.

**Algorithm**

**Step 1.** Given the initial displacement of the structure $u^0 \in W^S$, compute the characteristic function $\chi_{u^0}^S$, put $k := 0$.

**Step 2.** Find the velocity $v_{\epsilon} \in (H^1(D))^2$, $v_{\epsilon} = g$ on $\partial D$ and the pressure $p_{\epsilon}^k \in Q$ by solving the fluid problem

$$a_F(v_{\epsilon}^k, w) + b_F(w, p_{\epsilon}^k) + \frac{1}{\epsilon} \int_D \chi_{u_{\epsilon}^k}^S v_{\epsilon}^k \cdot w \, dx = \int_D f^F \cdot w \, dx, \quad \forall w \in W \\ b_F(v_{\epsilon}^k, q) = 0, \quad \forall q \in Q.$$

**Step 3.** Find the new displacement of the structure $u_{\epsilon}^{k+1} \in W^S$ by solving

$$a_S(u_{\epsilon}^{k+1}, w^S) = \int_{\Omega_0^S} (f^S - f^F) \cdot w^S \, dx + \int_{\Omega_0^S} 2\mu^F \epsilon (v_{\epsilon}^k) : \epsilon(w^S) \, dx \\ - \int_{\Omega_0^S} (\nabla \cdot w^S) p_{\epsilon}^k \, dx + \frac{1}{\epsilon} \int_{\Omega_0^S} (v_{\epsilon}^k \circ \varphi_{\epsilon}^k) \cdot w^S \, dx \quad \forall w^S \in W^S$$

where $\varphi_{\epsilon}^k(X) = X + u_{\epsilon}^k(X)$.

**Step 4.** Stopping test: if $\|u_{\epsilon}^k - u_{\epsilon}^{k+1}\|_{0, \Omega_0^S} \le tol$, then **Stop**.

**Step 5.** Compute the characteristic function $\chi_{u_{\epsilon}^{k+1}}^S$, put $k := k + 1$ and **Go to Step 2**.

It is possible to consider Navier-Stokes equations for the fluid domain, then at the **Step 2** we have to solve a non-linear system. In this case we can use the non-linear penalization term for $\alpha \neq 2$. If we use a non-linear model for the structure, too, we have to use more accurate approximations for the Jacobian determinant $J$ and the gradient of the deformation $F$.

## 4 Numerical tests

The numerical tests have been produced using the software FreeFem++ [3].
## 4.1 Test 1. Shell in steady-state cross flow

First, we have performed numerical simulation using a 2D model adapted from [1] (see Figure 2) where we have changed the physical parameters of the fluid and of the structure.

<figure><img src="image_5.png" /><figcaption>Fig. 2. Geometrical configuration for the Test 1.</figcaption></figure>

The dimensions of a rectangular elastic structure are: height $l = 3$ m, thickness $h = 0.125$ m. The computational domain of the fluid $D$ is a rectangle of height $H = 5$ m and length $L = 12$ m. The distance between the left side of the fluid and the left side of the structure is $2$ m. The lower left corner of Figure 2 is $(x_1 = 0, x_2 = 0)$.

We denote by $\Sigma_1, \Sigma_3$ the left and the right vertical boundaries of $D$ and by $\Sigma_2, \Sigma_4$ the bottom and the top boundaries of $D$, respectively.

The mechanical proprieties of the structure (polybutadiene) are: Young modulus $E^S = 1.6 \times 10^6$ N/m$^2$, Poisson's ratio $\nu^S = 0.49$, the applied volume forces on the structure $f^S: \Omega_0^S \to \mathbb{R}^2, f^S = (0,0) N/m^3$.

The dynamic viscosity of the fluid (glycerin) is $\mu^F = 1.14 N \cdot s/m^2$.

The inflow velocity profile on $\Sigma_1$ is

$$v_1(x_1, x_2) = V \times 1.5 \frac{(2H x_2 - x_2^2)}{H^2} \text{ m/s, } V = 1, \quad v_2(x_1, x_2) = 0.$$

The other boundary conditions are: $v = 0$ (no-slip) on $\Sigma_2$, $v \cdot n^F = 0$ (slip) on $\Sigma_4$ and $v \times n^F = 0$, $n^F \cdot (\sigma^F(v,p)n^F) = 0$ (the tangential velocity and the normal traction are zero) on $\Sigma_3$.

We use a fixed mesh for the fluid domain of 13096 triangles and 6719 vertices. The mesh of the structure domain has 188 triangles and 145 vertices. The fluid and structure meshes are not compatible, for example, a vertex on the structure boundary is not necessary a vertex on the fluid mesh. For the approximation of the fluid velocity and pressure we have employed the triangular finite elements $\mathbb{P}_1+\text{bubble}$ and $\mathbb{P}_1$ respectively, also called “mini” finite elements. The finite element $\mathbb{P}_1$ was used in order to solve the structure problem. The characteristic function was approached by $\mathbb{P}_0$ finite elements.
<figure><img src="image_6.png" /><figcaption>Fig. 3. The fluid pressure [Pa].</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. The fluid velocity around the final position of the structure. In each point of the grid, there is an arrow giving the direction of the velocity and the length of the arrow is proportional to the euclidean norm of the velocity. The maximal value for the horizontal component $v_1$ is 2.82 m/s and for the vertical component $v_2$ is 1.18 m/s</figcaption></figure>
We have performed the simulation using the **Algorithm** described in the previous section. For the stopping criterion at the **Step 4**, we have used the tolerance $tol = 0.2 \times 10^{-4}$. The penalization parameter is $\epsilon = 10^{-3}$. The stopping criterion holds after 5 iterations of the fixed point algorithm. The maximal structural displacement is $0.14 m$.

The fluid velocity into the fictitious domain is very small

$$ \|\mathbf{v}_{\epsilon}\|_{0,\Omega_{u_{\epsilon}}^{S}} = \sqrt{\int_{D} \chi_{u_{\epsilon}}^{S} \mathbf{v}_{\epsilon} \cdot \mathbf{v}_{\epsilon} d\mathbf{x}} = 0.0703384. $$

## 4.2 Test 2. Flexible appendix in a flow

We have adapted the benchmark from [4]. Originally, the structure was placed horizontally, parallel to the flow, but the displacements in this case are very small. We have placed the structure vertically, transversely to the flow, see Figure 5.

<figure><img src="image_8.png" /><figcaption>Fig. 5. Geometrical configuration.</figcaption></figure>

The structure is composed by a rectangular flexible appendix attached to a fixed circle. The circle center is positioned at $(0.2, 0.2) m$ measured from the left top corner of the channel. The circle has the radius $r = 0.5 m$ and the rectangular appendix is of length $\ell = 0.35 m$, thickness $h = 0.02 m$. The Young modulus is $E^S = 1.6 \times 10^6 N/m^2$ and Poisson's ratio is $\nu^S = 0.49$ (polybutadiene).

The channel has the length $L = 2.5 m$ and the width $H = 0.75 m$. The fluid dynamic viscosity is $\mu^F = 1.420 N \cdot s/m^2$ (glycerin).

We have used the following boundary conditions: at the inflow the velocity is

$$ v_1(x_1, x_2) = V \times 1.5 \frac{(H x_2 - x_2^2)}{(H/2)^2} \text{ m/s, } V = 1, \quad v_2(x_1, x_2) = 0; $$

at the bottom and the top we have imposed the no-slip boundary condition $\mathbf{v} = \mathbf{0}$ and at the outflow the traction free $\sigma^F(\mathbf{v}, p) \mathbf{n}^F = 0$.

We use a fixed mesh for the fluid domain of 30330 triangles and 15461 vertices and a structure mesh of 128 triangles and 97 vertices. We have employed the
same finite elements as for the Test 1. We have treated by the embedding domain technique only the flexible part of the structure.

The penalization parameter is $\epsilon = 10^{-4}$ and for the stopping criterion $tol = 10^{-8}$. The fixed point algorithm stops after 8 iterations.

The maximal horizontal displacement of the structure is 0.10886 m. The pressure and the velocity of the fluid are presented in Figures 6 and 7. The fluid velocity in the fictitious domain is very small

$$ \|\mathbf{v}_{\epsilon}\|_{0, \Omega_{u_{\epsilon}}^{S}} = \sqrt{\int_{D} \chi_{u_{\epsilon}}^{S} \mathbf{v}_{\epsilon} \cdot \mathbf{v}_{\epsilon} \, d\mathbf{x}} = 0.080920. $$

<figure><img src="image_9.png" /><figcaption>Fig. 6. The fluid pressure [Pa].</figcaption></figure>

## Conclusions

We have presented a fixed point algorithm for solving steady fluid-structure interaction problem. Using the embedding domain technique with penalization, the fluid equations as well as the structure equations are solved in fixed meshes. The fluid and structure meshes could be generated independently. The algorithm can be used for the three dimensions problems.
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** The fluid velocity [m/s] around the final position of the structure. In each point of the grid, there is an arrow giving the direction of the velocity. The length of the arrow is proportional to the euclidean norm of the velocity which is represented in the color bar.</figcaption></figure>

## Acknowledgment

The second author gratefully acknowledges support by Grant CNCS Romania 145/2011.

## References

1. Bathe, K.-J., Ladezma, G.: Benchmark problems for incompressible fluid flows with structural interactions, Comput. & Structures, 85, 628–644 (2007)

2. Halanay, A., Murea, C.M., Tiba, D.: Existence and approximation for a steady fluid-structure interaction problem using fictitious domain approach with penalization, accepted for publication in Mathematics and its Applications

3. Hecht, F.: FreeFem++, http://www.freefem.org

4. Turek, S., Hron, J.: Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow. In: Fluid-Structure Interaction - Modelling, Simulation, Optimization, Bungartz, H.-J., Schfer, M. (eds), 371–385, Lect. Notes Comput. Sci. Eng., 53, Springer, Berlin, (2006)