# Steady Fluid-Structure Interaction Using Fictitious Domain

Andrei Halanay, Cornel Marius Murea

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

## To cite this version:

Andrei Halanay, Cornel Marius Murea. Steady Fluid-Structure Interaction Using Fictitious Domain. 26th Conference on System Modeling and Optimization (CSMO), Sep 2013, Klagenfurt, Austria. pp.128-137, 10.1007/978-3-662-45504-3_12 . hal-01286405

HAL Id: hal-01286405

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

Submitted on 10 Mar 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Steady fluid-structure interaction using fictitious domain

Andrei Halanay<fnref n="1" />* and Cornel Marius Murea<fnref n="2" />

**Abstract.** We present a formulation for a steady fluid-structure interaction problem using fictitious domain technique with penalization. Numerical results are presented.

**Keywords:** fluid-structure interaction, fictitious domain

## 1 Setting for a steady fluid-structure interaction problem

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 $\bar{\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 $\bar{\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)$, then the boundary of the deformed structure is $\partial\Omega_u^S = \Gamma_D \cup \Gamma_u$ and the boundary of the fluid domain admits the decomposition $\partial\Omega_u^F = \partial D \cup \Gamma_D \cup \Gamma_u$. The fluid-structure geometrical configuration is represented in Figure 1.

Generally, the fluid equations are described using Eulerian coordinates, while for the structure equations, the Lagrangian coordinates are employed. The gradients with respect to the Eulerian coordinates $\mathbf{x} \in \Omega_u^S$ of a scalar field $q: D \to \mathbb{R}$ or a vector field $\mathbf{w} = (w_1, w_2) : D \to \mathbb{R}^2$ are denoted by $\nabla q$ and $\nabla \mathbf{w}$. The divergence operators with respect to the Eulerian coordinates of a vector field $\mathbf{w} = (w_1, w_2) : D \to \mathbb{R}^2$ and of a tensor $\sigma = (\sigma_{ij})_{1\le i,j\le 2}$ are denoted by $\nabla \cdot \mathbf{w}$ and $\nabla \cdot \sigma$.

<footnote n="1">¹ Department of Mathematics 1, University Politehnica of Bucharest
313 Splaiul Independenței, RO-060042, Bucharest, Romania
halanay@mathem.pub.ro</footnote>

<footnote n="2">² Laboratoire de Mathématiques, Informatique et Applications
Université de Haute Alsace
4-6, rue des Frères Lumière, 68093 Mulhouse Cedex, France
cornel.murea@uha.fr
http://www.edp.lmia.uha.fr/murea/</footnote>

<footnote>* Supported by Grant ID-PCE 2011-3-0211, Romania.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Geometrical configuration.</figcaption></figure>

When the derivatives are with respect to the Lagrangian coordinates $\mathbf{X} = \varphi^{-1}(\mathbf{x}) \in \Omega_0^S$, we use the notations: $\nabla_{\mathbf{X}}\mathbf{u}$, $\nabla_{\mathbf{X}} \cdot \mathbf{u}$, $\nabla_{\mathbf{X}} \cdot \sigma$.

If $\mathbf{A}$ is a square matrix, we denote by $\det \mathbf{A}$, $\mathbf{A}^{-1}$, $\mathbf{A}^T$ its determinant, the inverse and the transposed matrix, respectively. We write $\text{cof} \mathbf{A} = (\det \mathbf{A}) (\mathbf{A}^{-1})^T$ the co-factor matrix of $\mathbf{A}$. We write $\mathbf{A}^{-T} = (\mathbf{A}^{-1})^T$.

We denote by $\mathbf{F}(\mathbf{X}) = \mathbf{I} + \nabla_{\mathbf{X}}\mathbf{u}(\mathbf{X})$ the gradient of the deformation and by $J(\mathbf{X}) = \det \mathbf{F}(\mathbf{X})$ the Jacobian determinant, where $\mathbf{I}$ is the unit matrix.

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(v,p) = -p\mathbf{I} + 2\mu^F\epsilon(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.

The problem is to find the structure displacement $\mathbf{u}: \bar{\Omega}_0^S \to \mathbb{R}^2$, the fluid velocity $\mathbf{v}: \bar{\Omega}_u^F \to \mathbb{R}^2$ and the fluid pressure $p: \bar{\Omega}_u^F \to \mathbb{R}$ such that:

$$-\nabla_{\mathbf{X}} \cdot \sigma^S(\mathbf{u}) = \mathbf{f}^S, \quad \text{in } \Omega_0^S \qquad (1)$$

$$\mathbf{u} = 0, \quad \text{on } \Gamma_D \qquad (2)$$

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

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

$$\mathbf{v} = 0, \quad \text{on } \partial D \qquad (5)$$

$$\mathbf{v} = 0, \quad \text{on } \Gamma_D \qquad (6)$$

$$\mathbf{v} = 0, \quad \text{on } \Gamma_u \qquad (7)$$

$$\omega(\sigma^F(v,p)\mathbf{n}^F) \circ \varphi = -\sigma^S(\mathbf{u})\mathbf{n}^S, \quad \text{on } \Gamma_0 \qquad (8)$$

where $\mathbf{f}^S: \Omega_0^S \to \mathbb{R}^2$ are the applied volume forces on the structure and $\mathbf{n}^S$ is the structure unit outward vector normal to $\partial\Omega_0^S$. Similarly, we define $\mathbf{f}^F: \Omega_u^F \to \mathbb{R}^2$ and $\mathbf{n}^F$ the fluid unit outward vector normal to $\partial\Omega_u^S$.
We point out that the stress tensor of the structure is defined on the un-deformed structure domain $\Omega_0^S$, while the Cauchy stress tensors of the fluid is defined in the deformed domain $\Omega_u^F$.

We have used the notation $\omega(\mathbf{X}) = \|J\mathbf{F}^{-T}\mathbf{n}^S\|_{\mathbb{R}^2} = \|\mathrm{cot}(\mathbf{F})\mathbf{n}^S\|_{\mathbb{R}^2}$ for $\mathbf{X}$ on $\partial\Omega_0^S$, which is a kind of Jacobian determinant for the change of variable formula for integral over surface.

The equations (1), (2) concern the structure, while (3)-(6) concern the fluid. The equations (7), (8) represent the boundary conditions on the moving fluid-structure interface. The fluid and the structure domains $\Omega_u^F, \Omega_u^S$ depend on the structure displacement $u$ which is unknown.

## 2 Parametrization and regularization of the characteristic function

The regularization of the characteristic function of the deformed structure domain is necessary in order to prove the continuity of the solution with respect to the structure displacement, [3].

Denote by $\|\cdot\|_{1,\infty,\Omega}$ the usual norm of the Sobolev space $W^{1,\infty}(\Omega)$ and by $\|\cdot\|_{m,\Omega}$ the usual norm of $H^m(\Omega)$, $m \ge 0$ with the convention $H^0(\Omega) = L^2(\Omega)$.
For every $0 < \delta < 1$, there exists $0 < \eta_\delta < 1$ such that

$$1 - \delta \le \det (\mathbf{I} + \nabla \mathbf{u}) \le 1 + \delta, \quad \text{a.e. } \mathbf{x} \in \Omega_0^S \quad (9)$$

for all $\mathbf{u} \in (W^{1,\infty}(\Omega_0^S))^2$ that satisfy $\|\mathbf{u}\|_{1,\infty,\Omega_0^S} \le \eta_\delta$.

We define

$$B_{\delta} = \{u \in W^{1,\infty}(\Omega_{0}^{S})^{2}; \|u\|_{1,\infty,\Omega_{0}^{S}} \le \eta_{\delta}, u = 0 \text{ on } \Gamma_{D}\}. \quad (10)$$

Let $j \in W^{1,\infty}(D)$ be a parametrization of $\Omega_0^S \subset D$, i.e. :

$$j(\mathbf{x}) > 0, \mathbf{x} \in \Omega_0^S, \quad j(\mathbf{x}) < 0, \mathbf{x} \in D \setminus \bar{\Omega}_0^S, \quad j(\mathbf{x}) = 0, \mathbf{x} \in \partial\Omega_0^S.$$

The parametrization is not necessarily unique.

Let $\mathbf{u} \in B_\delta$ be a given structure displacement. Denote, as before, $\Omega_u^S = \varphi(\Omega_0^S)$, where $\varphi(\mathbf{X}) = \mathbf{X} + \mathbf{u}(\mathbf{X})$. Then $\varphi: \bar{\Omega}_0^S \to \bar{\Omega}_u^S$ is bijective and bilipschitzian and

$$j_u(\mathbf{y}) = \begin{cases} j(\mathbf{x}), & \mathbf{y} = \varphi(\mathbf{x}) \in \Omega_u^S \\ 0, & \mathbf{y} \in \partial\Omega_u^S \\ -\mathrm{dist}(\mathbf{y}, \bar{\Omega}_u^S), & \mathbf{y} \notin \bar{\Omega}_u^S \end{cases}$$

is a parametrization of $\Omega_u^S$, $j_u \in W^{1,\infty}(D)$.

If $H$ is the Heaviside function $H: \mathbb{R} \to \{0, 1\}$,

$$H(r) = \begin{cases} 1, & r \ge 0 \\ 0, & r < 0 \end{cases}$$
then $H(j_u(\cdot))$ is the characteristic function of $\Omega_u^S$.

For $\epsilon > 0$, let $\Omega_0^\epsilon \subset \subset \Omega_0^S$. Since $j : D \to \mathbb{R}$ is Lipschitz continuous and $j > 0$ in $\Omega_0^S$, there is $\mu_\epsilon > 0$ such that $j(\mathbf{x}) \ge \mu_\epsilon > 0$, for all $\mathbf{x} \in \Omega_0^\epsilon$. Consequently,

$$ \mu_{\epsilon} \leq \min_{\mathbf{y} \in \Omega_{u}^{\epsilon}} j_{u}(y), \quad \forall \mathbf{u} \in B_{\delta}. $$

Then we take $H^{\mu_\epsilon}$, the Yosida regularization of $H$

$$ H^{\mu_{\epsilon}}(r) = \begin{cases} 1, & r \ge \mu_{\epsilon} \\ \frac{r}{\mu_{\epsilon}} & 0 \le r < \mu_{\epsilon} \\ 0, & r < 0 \end{cases} $$

and we set $\tilde{H}_u(\mathbf{x}) = H^{\mu_\epsilon}(j_u(\mathbf{x}))$ for all $\mathbf{x} \in D$, which is a Lipschitz regularization of the characteristic function of $\Omega_u^S$. We have constructed $\tilde{H}_u : \bar{D} \to \mathbb{R}$, Lipschitz on $\bar{D}$, $0 \le \tilde{H}_u(\mathbf{x}) \le 1$, for all $\mathbf{x}$ in $\bar{D}$ such that

$$ \tilde{H}_u(\mathbf{x}) = \begin{cases} 0, & \mathbf{x} \in \bar{D} \setminus \Omega_u^S \\ 1, & \mathbf{x} \in \Omega_u^\epsilon \end{cases} \qquad (11) $$

where $\Omega_u^\epsilon \subset \subset \Omega_u^S$.

### 3 Weak formulation using fictitious domain technique with penalization

We assume that $D$ and $\Omega_0^S$ are Lipschitz. 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 \epsilon(\mathbf{u}) : \epsilon(\mathbf{w}^S) \, d\mathbf{X} $$

$$ a_F(\mathbf{v}, \mathbf{w}) = \int_D 2\mu^F \epsilon(\mathbf{v}) : \epsilon(\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 that $\mathbf{f}^F \in (L^2(D))^2$, $\mathbf{f}^S \in (L^2(\Omega_0^S))^2$.
## Weak fluid formulation using fictitious domain

For a given $\mathbf{u} \in B_\delta$, we define: fluid velocity $\mathbf{v}_\epsilon \in W$ and fluid pressure $p_\epsilon \in Q$, as the solution of the following problem:

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

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

## Weak structure formulation

For given $\mathbf{u} \in B_\delta$, $\mathbf{v}_\epsilon \in W$ and $p_\epsilon \in Q$, we define the structure displacement $\mathbf{u}_\epsilon \in W^S$ as the solution of

$$ a_S(\mathbf{u}_\epsilon, \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}_\epsilon, p_\epsilon) \circ \varphi) \mathbf{F}^{-T} : \nabla_X \mathbf{w}^S d\mathbf{X} \\ + \frac{1}{\epsilon} \int_{\Omega_0^S} J \tilde{H}_u(\varphi) ((\mathbf{v}_\epsilon \circ \varphi) \cdot \mathbf{w}^S + (\nabla \mathbf{v}_\epsilon \circ \varphi) \mathbf{F}^{-T} : \nabla_X \mathbf{w}^S) d\mathbf{X} \\ - \int_{\Omega_0^S} J(\mathbf{f}^F \circ \varphi) \cdot \mathbf{w}^S d\mathbf{X}, \forall \mathbf{w}^S \in W^S \quad (14) $$

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

*Remark 1.* From the structure equation $-\nabla \cdot \sigma^S(\mathbf{u}_\epsilon) = \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}_\epsilon, \mathbf{w}^S) = \int_{\Omega_0^S} \mathbf{f}^S \cdot \mathbf{w}^S d\mathbf{X} + \int_{\Gamma_0} \sigma^S(\mathbf{u}_\epsilon) \mathbf{n}^S \cdot \mathbf{w}^S dS. $$

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

For each $i \in \mathbb{N}^*$, there exists an unique eigenvalue $\lambda_i > 0$ and an unique eigenfunction $\phi^i \in W^S$, solution of

$$ a_S(\phi^i, \mathbf{w}^S) = \lambda_i \int_{\Omega_0^S} \phi^i \cdot \mathbf{w}^S d\mathbf{X}, \quad \forall \mathbf{w}^S \in W^S \quad (15) $$

such that

$$ \int_{\Omega_0^S} \phi^i \cdot \phi^j d\mathbf{X} = \delta_{ij}, \quad (16) $$
see [8], chap. 6. We assume that $\lambda_1 \le \lambda_2 \le \dots$. The set $\{\phi^i, i \in \mathbb{N}^*\}$ forms an orthonormal basis of $L^2(\Omega_0^S)$. Let $m \in \mathbb{N}^*$ be given. Let $\mathbf{u}_\varepsilon^m$ be the orthogonal projection of $\mathbf{u}_\varepsilon$ on $\operatorname{span}(\phi^i, i = 1, \dots, m)$ in $L^2(\Omega_0^S)$, so $\mathbf{u}_\varepsilon^m = \sum_{i=1}^m \alpha_i \phi^i$, $\alpha_i \in \mathbb{R}$.

We define

$$B_{\delta}^{m} = \left\{ \mathbf{u} \in (W^{1,\infty}(\Omega_{0}^{S}))^{2}; \mathbf{u} = 0 \text{ on } \Gamma_{D}, \|\mathbf{u}\|_{1,\infty,\Omega_{0}^{S}} < \eta_{\delta}, \right. \\ \left. \mathbf{u} = \sum_{i=1}^{m} \alpha_{i} \phi^{i}, \alpha_{i} \in \mathbb{R} \right\}. \quad (17)$$

We have $B_{\delta}^{m} \subset B_{\delta}$. For each $\mathbf{u} \in B_{\delta}$, we define the nonlinear operator

$$T_{\varepsilon}^{m}(\mathbf{u}) = \mathbf{u}_{\varepsilon}^{m}.$$

A solution of the penalized fluid-structure interaction problem will be, by definition, a fixed point of $T_{\varepsilon}^{m}$ in $B_{\delta}^{m}$.

*Remark 2.* As in [3], we can prove the existence of a solution of the penalized fluid-structure interaction problem $T_{\varepsilon}^{m}(\mathbf{u}_{\varepsilon}^{m}) = \mathbf{u}_{\varepsilon}^{m}$ using the Schauder fixed point theorem. In order to obtain supplementary regularity of the Stokes equations, a non linear penalization term $\frac{1}{\varepsilon} \int_{D} \tilde{H}_{u} sgn(\mathbf{v}_{\varepsilon}) |\mathbf{v}_{\varepsilon}|^{\alpha-1} \cdot \mathbf{w} \, d\mathbf{x}$ was employed in [3], where $\alpha > 2$. In the present paper, we use a linear penalization term in (12), but we are working only with a finite number of eigenfunctions of the structure equations. The behavior of $\mathbf{u}_{\varepsilon}^{m}$ when $\varepsilon$ goes to zero will be studied in [4].

# 4 Fixed point iterations

We replace $\tilde{H}_u$ in (12) and (14) 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 to obtain 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 unit matrix $I$.

## Algorithm 1 by fixed point iterations

**Step 1.** Given the initial displacement of the structure $\mathbf{u}^{m,0} = \sum_{i=1}^{m} \alpha_i^0 \phi^i$, compute the characteristic function $\chi_{u^0}^S$, put $k := 0$.

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

$$a_F(\mathbf{v}_\varepsilon^k, \mathbf{w}) + b_F(\mathbf{w}, p_\varepsilon^k) \\ + \frac{1}{\varepsilon} \int_D \chi_{u_\varepsilon^k}^S (\mathbf{v}_\varepsilon^k \cdot \mathbf{w} + \nabla \mathbf{v}_\varepsilon^k : \nabla \mathbf{w}) \, d\mathbf{x} = \int_D \mathbf{f}^F \cdot \mathbf{w} \, d\mathbf{x}, \quad \forall \mathbf{w} \in W \\ b_F(\mathbf{v}_\varepsilon^k, q) = 0, \quad \forall q \in Q.$$
**Step 3.** Find the new displacement of the structure $\mathbf{u}_{\varepsilon}^{m,k+1} = \sum_{i=1}^{m} \alpha_{i}^{k+1} \phi^{i}$ by solving

$$ \alpha_i^{k+1} \lambda_i = \int_{\Omega_0^S} (\mathbf{f}^S - \mathbf{f}^F) \cdot \phi^i \, d\mathbf{x} \\ + \int_{\Omega_0^S} 2\mu^F \epsilon (\mathbf{v}_{\varepsilon}^k) : \epsilon (\phi^i) \, d\mathbf{x} - \int_{\Omega_0^S} (\nabla \cdot \phi^i) p_{\varepsilon}^k \, d\mathbf{x} \\ + \frac{1}{\varepsilon} \int_{\Omega_0^S} ((\mathbf{v}_{\varepsilon}^k \circ \varphi_{\varepsilon}^k) \cdot \phi^i + (\nabla \mathbf{v}_{\varepsilon}^k \circ \varphi_{\varepsilon}^k) : \nabla \phi^i) \, d\mathbf{x}, \quad i = 1, \dots, m $$

where $\varphi_{\varepsilon}^{k}(\mathbf{X}) = \mathbf{X} + \mathbf{u}_{\varepsilon}^{m,k}(\mathbf{X})$.

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

**Step 5.** Compute the characteristic function $\chi_{\mathbf{u}_{\varepsilon}^{m,k+1}}^{S}$, put $k := k + 1$ and **Go to Step 2**.

## 5 Least squares approach

The previous algorithm converges if the operator $T_{\varepsilon}^{m}$ is a contraction. But, the **Algorithm 1** fails for some physical parameters. For this reason, we introduce a second algorithm which is more robust.

For $\alpha = (\alpha_1, \dots, \alpha_m) \in \mathbb{R}^m$, we can define $\beta = (\beta_1, \dots, \beta_m) \in \mathbb{R}^m$ by

$$ T_{\varepsilon}^{m} \left( \sum_{i=1}^{m} \alpha_{i} \phi^{i} \right) = \sum_{i=1}^{m} \beta_{i} \phi^{i}. $$

We set the cost function $J(\alpha) = \frac{1}{2} \sum_{i=1}^{m} (\alpha_i - \beta_i)^2$ and now the problem to be solved is $\inf_{\alpha \in \mathbb{R}^m} J(\alpha)$. In order to solve the optimization problem, we employ the quasi-Newton iterative method called Broyden, Fletcher, Goldforb, Shano (BFGS) scheme (see for example [2], chap. 9).

**Algorithm 2 by the BFGS method**

**Step 0.** Choose a starting point $\alpha^0 \in \mathbb{R}^m$, an $m \times m$ symmetric positive matrix $H_0$. Set $k = 0$.

**Step 1.** Compute $\nabla J(\alpha^k)$.

**Step 2.** If $\|\nabla J(\alpha^k)\| < tol$ **Stop**.

**Step 3.** Set $\mathbf{d}^k = -H_k \nabla J(\alpha^k)$.

**Step 4.** Determine $\alpha^{k+1} = \alpha^k + \theta_k \mathbf{d}^k$, $\theta_k > 0$ by means of an approximate minimization

$$ J(\alpha^{k+1}) \approx \min_{\theta \ge 0} J(\alpha^k + \theta \mathbf{d}^k). $$

**Step 5.** Compute $\delta_k = \alpha^{k+1} - \alpha^k$.

**Step 6.** Compute $\nabla J(\alpha^{k+1})$ and $\gamma_k = \nabla J(\alpha^{k+1}) - \nabla J(\alpha^k)$.
**Step 7. Compute**

$$H_{k+1} = H_k + \left( 1 + \frac{\gamma_k^T H_k \gamma_k}{\delta_k^T \gamma_k} \right) \frac{\delta_k \delta_k^T}{\delta_k^T \gamma_k} - \frac{\delta_k \gamma_k^T H_k + H_k \gamma_k \delta_k^T}{\delta_k^T \gamma_k}$$

**Step 8.** Update $k = k + 1$ and go to the **Step 2**.

The matrices $H_k$ approach the inverse of the Hessian of $J$. For the inaccurate line search at the **Step 4**, the methods of Goldstein and Armijo were used. If we denote by $g: [0, \infty) \to \mathbb{R}$ the function $g(\theta) = J(\alpha^k + \theta \mathbf{d}^k)$, we determine $\theta_k > 0$ such that: $g(0) + (1 - \lambda) \theta_k g'(0) \le g(\theta_k) \le g(0) + \lambda \theta_k g'(0)$, where $\lambda \in (0, 1/2)$.

In this paper, we compute $\nabla J(\alpha)$ by the Finite Differences Method $\frac{\partial J}{\partial \alpha_k}(\alpha) \approx (J(\alpha + \Delta\alpha_k e_k) - J(\alpha)) / \Delta\alpha_k$, where $e_k$ is the $k$-th vector of the canonical base of $\mathbb{R}^m$ and $\Delta\alpha_k > 0$ is the grid spacing.

Concerning the convergence rate, the fixed point algorithm is slower than the BFGS Method. If the starting point is not sufficiently close to the solution, the fixed point algorithm diverges. On the contrary, the BFGS Method is less sensitive to the choice of the starting point and, in general, it is convergent to a local minimizer from almost any starting point. This is the main advantage, (see [6]).

# 6 Numerical results. Deformation of a tall building under the action of wind

We have performed numerical simulations using a 2D model adapted from [1] (see Figure 2).

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Geometrical configuration for the numerical results</figcaption></figure>

The dimensions of a rectangular tall building are: height $H = 180$ m, length $L = 30$ m. The computational domain of the fluid $D$ is a rectangle of height
$H_1 = 3H$ and length $L_1 = L + 4H$, its left bottom corner is at (0,0). We shall allow nonhomogeneous Dirichlet data in the numerical experiments.

The distance between the left side of the fluid and the left side of the structure is $H$. We denote by $\Sigma_1, \Sigma_3$ the left and the right vertical boundaries and by $\Sigma_2, \Sigma_4$ the bottom and the top boundaries, respectively.

The mechanical proprieties of the building assumed to be an elastic structure are: Young modulus $E^S = 2.3 \times 10^5 \text{ N/m}^2$, Poisson's ratio $\nu^S = 0.25$, the applied volume forces on the structure $\mathbf{f}^S: \Omega_0^S \to \mathbb{R}^2$, $\mathbf{f}^S = (0,0) \text{ N/m}^3$. If the Young modulus is $E^S = 2.3 \times 10^8 \text{ N/m}^2$ as in [1], the displacements of the structure are very small.

The fluid is the air with: dynamic viscosity $\mu^F = 7.03 \times 10^{-2} \text{ N} \cdot \text{s/m}^2$, the applied volume forces on the fluid $\mathbf{f}^S: D \to \mathbb{R}^2$, $\mathbf{f}^F = (0,0) \text{ N/m}^3$. The inflow velocity profile is $g(x_1, x_2) = 100 \left(\frac{x_2}{H}\right)^{0.19} \text{ m/s}$. The considered boundary conditions for the fluid are more natural from the point of view of applications and differ slightly compared with the previous sections. We impose: $\mathbf{v}_\epsilon = \mathbf{g}$ on $\Sigma_1 \cup \Sigma_4$, $\mathbf{v}_\epsilon = 0$ on $\Sigma_2$ and $\sigma^F(\mathbf{v}, p) \mathbf{n}^F = 0$ on $\Sigma_3$.

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** The fixed mesh of the fluid domain (left). The finer zone $[H-L, H+2L] \times [0, H+L]$ of the fluid mesh covers the structure mesh which occupies initially the rectangle $[H, H+L] \times [0, H]$ (zoom, right). 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 (right).</figcaption></figure>

The numerical tests have been produced using the software FreeFem++ [5]. For the approximation of the fluid velocity and pressure we have employed the triangular finite elements $P_1+bubble$ and $P_1$ respectively on a mesh of 34871 triangles and 17550 vertices. The finite element $P_1$ was used in order to solve the structure problem on a mesh of 192 triangles and 125 vertices. The characteristic function was approached by $P_0$ finite element.

We have performed the simulation using the **Algorithm 2** described in the previous section. We have used the initial displacement $\alpha^0 = 0$ at the **Step 0** and the tolerance $tol = 0.0001$ for the stopping criterion at the **Step 2**. The penalization parameter is $\varepsilon = 0.001$ and the number of the eigenfunctions is
$m = 5$. The stopping criterion holds after 6 iterations of the BFGS algorithm, the initial value of the cost function is 34.30 and the final value is $5.9 \times 10^{-13}$. The maximal structural displacement is 0.148 $m$.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Velocity (left) and pressure (right) of the fluid around the deformed structure.</figcaption></figure>

The fluid velocity is almost zero in the deformed structure domain, more precisely $\|\mathbf{v}_\epsilon\|_{1,\Omega_{u_\epsilon}^S} = \sqrt{\int_D \chi_{u_\epsilon}^S (\mathbf{v}_\epsilon \cdot \mathbf{v}_\epsilon + \nabla\mathbf{v}_\epsilon : \nabla\mathbf{v}_\epsilon) d\mathbf{x}} = 0.00555$.

References

1. Braun, A.L., Awruch, A.M.: Aerodynamic and aeroelastic analyses on the CAARC standard tall building model using numerical simulation, Computers and Structures 87, 564-581, (2009).

2. Dennis Jr., J.E., Schnabel, R.B.: Numerical methods for unconstrained optimization and nonlinear equations. Classics in Applied Mathematics, 16. Society for Industrial and Applied Mathematics, Philadelphia, PA, (1996).

3. Halanay, A., Murea, C.M., Tiba, D.: Existence and approximation for a steady fluid-structure interaction problem using fictitious domain approach with penalization, Mathematics and its Applications 5, no. 1-2, 120-147, (2013).

4. Halanay, A., Murea, C.M.: Existence of a steady flow of Stokes fluid past a linear elastic structure using fictitious domain, in preparation

5. Hecht, F.: http://www.freefem.org

6. Murea, C.M.: Numerical simulation of a pulsatile flow through a flexible channel, ESAIM: Math. Model. Numer. Anal., 40, no. 6, 1101-1125, (2006).

7. Murea, C.M., Halanay, A.: Embedding domain technique for a fluid-structure interaction problem, Modeling and Optimization, D. Hömberg, F. Tröltzsch, (Eds.), Springer, pp. 358-367, (2013).

8. Raviart, P.-A., Thomas, J.-M.: Introduction à l'analyse numérique des équations aux dérivées partielles, Dunod, Paris (1998).