# Préconditionnement d'algorithmes d'interaction fluide-structure

Miguel Angel Fernández, Jimmy Mullaert, Yvon Maday

► To cite this version:

Miguel Angel Fernández, Jimmy Mullaert, Yvon Maday. Préconditionnement d'algorithmes d'interaction fluide-structure. 10e colloque national en calcul des structures, May 2011, Giens, France. hal-00592797

HAL Id: hal-00592797

https://hal.science/hal-00592797v1

Submitted on 3 May 2011

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Préconditionnement d'un algorithme partitionné d'interaction fluide-structure

M.A. Fernández<fnref n="1" />, Y. Maday<fnref n="2" />, J. Mullaert<fnref n="3" />

**Résumé** — La simulation numérique de l'interaction entre un fluide incompressible circulant dans un solide pose des difficultés. L'effet de masse ajoutée perturbe la convergence des algorithmes Dirichlet-Neumann dans la résolution des schémas de couplage implicites (ou semi-implicites). La relaxation de la condition de couplage cinématique ou de l'incompressibilité à l'interface permet de rendre ces itérations moins sensibles à cet effet. Nous proposons de regrouper ces méthodes dans un cadre de préconditionnement non-linéaire commun et d'y ajouter des nouveaux algorithmes.

**Mots clés** — interaction fluide-structure, approche partitionnée, préconditionnement, masse ajoutée.

## 1 Introduction

Dans cet article, on étudie la simulation numérique d'un problème d'interaction fluide-structure mettant en oeuvre un fluide newtonien incompressible et une structure élastique. En pratique, il peut par exemple s'agir d'une portion de grosse artère dans laquelle s'écoule du sang (voir [4]). Pour ce type de problème, il est connu que les algorithmes de couplage Dirichlet-Neumann explicites peuvent être inconditionnellement instables suivant la quantité de masse ajoutée, dictée par la géométrie du problème et le rapport des densités fluide et solide (voir [2]).

Pour contourner cette difficultés, nous partons d'un schéma dans lequel le couplage est implicite et qui nécessite donc des itérations à l'intérieur de chaque pas de temps. Les itérations Dirichlet-Neumann doivent être relaxées pour assurer leur convergence et celle-ci est très fortement affectée par l'effet de masse ajoutée (voir [2, 8]). Nous proposons donc de préconditionner ce schéma dans le cadre d'itérations de défauts corrigés pour en améliorer les propriétés de convergence. Ceci permet de rassembler certaines méthodes itératives existantes dans un même cadre et d'en proposer de nouvelles.

Nous commencerons donc par présenter le problème continu, sa discrétisation en temps pour ensuite introduire quatre différent préconditionneurs dont les propriétés sont illustrées numériquement dans la dernière section. Le choix du paramètre de relaxation est aussi analysé pour certains de ces algorithmes.

## 2 Le problème continu

On considère un domaine de référence $\Omega = \Omega^f \cup \Omega^s$ et on note $\Sigma \stackrel{\text{def}}{=} \partial\Omega^f \cap \partial\Omega^s$ l'interface entre le fluide et le solide et $\Gamma$ les bords d'entrée et de sortie du fluide. On désigne par $\Gamma^d$ et $\Gamma^n$ des parties de la frontières du solide en dehors de l'interface fluide-structure. On introduit le déplacement $d^f : \Omega^f \times \mathbb{R}^+ \to \mathbb{R}^d(d = 2,3)$ du domaine fluide et la carte ALE $\mathcal{A} \stackrel{\text{def}}{=} I_{\Omega^f} + d^f$, de telle sorte que le domaine fluide mobile s'écrive $\Omega^f(t) = \mathcal{A}(\Omega^f, t)$ En pratique, on obtient $d^f$ par un relèvement du déplacement solide $d : \Omega^s \times \mathbb{R}^+ \to \mathbb{R}^d$ de telle sorte que $d^f = \text{Ext}(d|_{\Sigma})$. On note $\Sigma(t) \stackrel{\text{def}}{=} \partial\Omega^s(t) \cap \partial\Omega^f(t)$, la position de l'interface fluide-structure à l'instant $t$.

Le problème d'interaction fluide-structure que l'on considère dans cet article consiste à trouver la vitesse fluide $u : \Omega^f \times \mathbb{R}^+ \to \mathbb{R}^d$, la pression $p : \Omega^f \times \mathbb{R}^+ \to \mathbb{R}$ et le déplacement solide $d : \Omega^s \times \mathbb{R}^+ \to \mathbb{R}^d$ qui vérifient

<footnote n="1">¹ INRIA, CRI Paris-Rocquencourt, France, miguel.fernandez@inria.fr</footnote>

<footnote n="2">² LJLL, Université Pierre et Marie Curie, France, maday@ann.jussieu.fr</footnote>

<footnote n="3">³ INRIA, CRI Paris-Rocquencourt, France, jimmy.mullaert@inria.fr</footnote>
$$
\left\{
\begin{array}{ll}
\rho^f \partial_t u|_{\mathcal{A}} + \rho^f(u-w) \cdot \nabla u - \nabla \cdot \sigma^f(u,p) = 0 & \text{dans } \Omega^f(t), \\
\nabla \cdot u = 0 & \text{dans } \Omega^f(t), \\
\sigma^f(u,p)n^f = -\bar{p}n^f & \text{sur } \Gamma, \\
\rho^s \partial_{tt} d - \nabla \cdot \Pi(d) = 0 & \text{dans } \Omega^s, \\
d = 0 & \text{sur } \Gamma^d, \\
\Pi(d)n^s = 0 & \text{sur } \Gamma^n, \\
\Omega^f(t) = (I_{\Omega^f} + d^f)(\Omega^f), \quad d^f = \text{Ext}(d|_{\Sigma}), \quad w = \partial_t d^f & \text{dans } \Omega^f, \\
u = \partial_t d & \text{sur } \Sigma(t), \\
\Pi(d)n^s = -J^f \sigma^f(u,p)(F^f)^{-T} n^f & \text{sur } \Sigma.
\end{array}
\right.
\tag{1}
$$

Ici, une pression $\bar{p}$ est imposée sur le bord d’entrée-sortie $\Gamma$, $\rho^f$ et $\rho^s$ représentent respectivement les densités fluide et solide. La dérivée ALE est notée $\partial_t|_{\mathcal{A}}$ et le tenseur des contraintes de Cauchy $\sigma^f(u,p) \stackrel{\text{def}}{=} -p\boldsymbol{I} + 2\mu\epsilon(u)$ avec $\mu$ la viscosité dynamique, $\epsilon(u) \stackrel{\text{def}}{=} \frac{1}{2}(\nabla u + \nabla u^T)$ le tenseur du taux de déformation, $\Pi(d)$ le premier tenseur de Piola-Kirchhoff de la structure, $F^f \stackrel{\text{def}}{=} \nabla \mathcal{A}$ le gradient de déformation du domaine fluide et $J^f \stackrel{\text{def}}{=} \det F^f$ son Jacobien. $n^f$ et $n^s$ sont les normales unitaires des domaines fluide et solide.

3 La discrétisation en temps

Nous obtenons un problème discret en temps à partir du problème continu présenté ci-dessus en
utilisant pour le fluide un schéma d'Euler implicite avec un traitement semi-implicite du terme convectif
non linéaire et pour le solide un schéma du point milieu qui présente l'avantage de conserver l'énergie
élastique, tout comme le problème continu. Dans la suite, on notera $\tau$ le pas de temps discret, $D_{\tau}X^{n+1} \stackrel{\text{def}}{=} (X^{n+1} - X^n)/\tau$ la dérivée discrète et $X^{n+1/2} \stackrel{\text{def}}{=} (X^{n+1} + X^n)/2$ l'approximation du point milieu.

Pour simplifier les notations par la suite, on introduit l'opérateur Neumann-Dirichlet, qui fait corres-
pondre à un couple (u, p) le déplacement de l'interface $\lambda = d|_{\Sigma}$, solution du problème suivant :

$$
\left\{
\begin{array}{ll}
\frac{\rho^s}{\tau}(\dot{d} - \dot{d}^n) - \nabla \cdot \Pi((d+d^n)/2) = 0, & \text{dans } \Omega^s, \\
\frac{1}{\tau}(d-d^n) = \frac{1}{2}(\dot{d} + \dot{d}^n), & \text{dans } \Omega^s, \\
d = 0 & \text{sur } \Gamma^d, \\
\Pi((d+d^n)/2)n^s = 0, & \text{sur } \Gamma^n, \\
\Pi((d+d^n)/2)n^s = -J^{f,n+1}\sigma^f(u,p)(F^{f,n+1})^{-T}n^s, & \text{sur } \Sigma.
\end{array}
\right.
$$

Ainsi, on peut écrire simplement par $\lambda = S(u, p)$, le déplacement de l'interface induit par une contrainte fluide exercée sur le solide. Le changement d'un modèle solide se traduit uniquement par un changement de l'opérateur $S$, ce qui permet de s'affranchir du modèle solide utilisé dans le reste de l'article.

On utilise un schéma de couplage dans lequel la continuité des efforts et de la vitesse sont traités de manière implicite et la géométrie de manière explicite. On suppose que $\Omega^n$, $u^n$, $p^n$, $d^n$ and $\dot{d}^n$ sont donnés et on propose de calculer $\Omega^{n+1}$, $u^{n+1}$, $p^{n+1}$, $d^{n+1}$ et $\dot{d}^{n+1}$ de la manière suivante :

On définit l'espace $V \stackrel{\text{def}}{=} (H_{\Sigma}^1(\Omega^f))^3 \times L^2(\Omega^f) \times (L^2(\Sigma))^3$ et l'opérateur

$$
A : (H^1(\Omega^f))^3 \times L^2(\Omega^f) \times (L^2(\Sigma))^3 \to V',
$$
défini par

$$
\begin{equation}
\begin{split}
\langle A(\mathbf{u}, p, \boldsymbol{\lambda}), (\mathbf{v}, q, \boldsymbol{\xi}) \rangle_{V', V} &\stackrel{\text{def}}{=} \frac{\rho^f}{\tau} \left( \int_{\Omega^{f,n+1}} \mathbf{u} \cdot \mathbf{v} - \int_{\Omega^{f,n}} \mathbf{u}^n \cdot \mathbf{v} \right) - \rho^f \int_{\Omega^{f,n+1}} (\nabla \cdot \mathbf{w}^{n+1}) \mathbf{u} \cdot \mathbf{v} \right. \\
&\qquad \left. + \rho^f \int_{\Omega^{f,n+1}} (\mathbf{u}^n - \mathbf{w}^{n+1}) \cdot \nabla \mathbf{u} \cdot \mathbf{v} + \int_{\Omega^{f,n+1}} \boldsymbol{\sigma}(\mathbf{u}, p) : \boldsymbol{\epsilon}(\mathbf{v}) \right) \\
&\qquad + \int_{\Omega^{f,n+1}} q \nabla \cdot \mathbf{u} + \int_{\Sigma^{n+1}} \left( \mathbf{u} - \frac{\boldsymbol{\lambda} - \mathbf{d}^n}{\tau} \right) \cdot \boldsymbol{\xi} + \int_{\Gamma} \bar{p} \mathbf{v} \cdot \mathbf{n}^f,
\end{split}
\tag{2}
\end{equation}
$$

pour tout $(v, q, \xi) \in V$. Finalement, pour simplifier les notation, notre problème semi-discret en temps
revient à trouver $(u^{n+1}, p^{n+1}) \in (H^1(\Omega^f))^3 \times L^2(\Omega^f)$ qui annule le résidu suivant :

$$
R(\boldsymbol{u}^{n+1}, \boldsymbol{p}^{n+1}) \stackrel{\text{def}}{=} A(\boldsymbol{u}^{n+1}, \boldsymbol{p}^{n+1}, S(\boldsymbol{u}^{n+1}, \boldsymbol{p}^{n+1})) = 0 \quad \text{dans } V'. \quad (3)
$$

4 Stratégies de préconditionnement

Pour annuler un résidu $R(u, p)$, on va utiliser des itérations de type "défauts corrigés" (voir [10]), qui
utilisent comme préconditionneur un résidu simplifié $\tilde{R}$ proche de celui que l'on souhaite annuler, mais
bien plus facile à résoudre. Le fil conducteur de cet article sera l'itération suivante :

$$
\tilde{R}(\mathbf{u}_{k+1}, p_{k+1}) = \tilde{R}(\mathbf{u}_k, p_k) - R(\mathbf{u}_k, p_k) \quad \forall k \ge 0, \quad (4)
$$

On verra que cela revient à relaxer la condition de continuité des vitesse ou l'incompressibilité du fluide, facteurs reconnus de l'effet de masse ajoutée.

4.1 Algorithme 1 : relaxation de la condition de continuité des vitesses

Dans cette partie, nous allons utiliser un modèle de solide plus simple que celui du problème couplé initial. On définit l'opérateur linéaire $\tilde{S}$, destiné à remplacer l'opérateur non linéaire Dirichlet-Neumann, et qui fait intervenir un réel positif $\alpha$ qui représente la compliance de notre solide simplifié.

$$
\tilde{S}(u, p) \stackrel{\text{def}}{=} -\alpha \sigma(u, p) n^f, \quad (5)
$$

Ce modèle relie linéairement la contrainte que le fluide exerce sur l'interface Σ et le déplacement
de cette dernière. Le solide se comporte donc comme un ressort sans masse de raideur 1/α. Dans la
pratique, on peut ajuster ce paramètre pour que le solide simplifié ressemble le plus possible au modèle
solide complet dans son comportement, ce qui pose un problème d'estimation. Nous verrons qu'il existe
une autre interprétation de ce paramètre et que l'on peut lui donner une valeur naturelle en fonction des
données physiques et géométriques du problème. Finalement, on introduit le résidu simplifié suivant,
dans lequel on a remplacé le modèle solide initial par notre solide simplifié.

$$
\tilde{R}_I(u, p) = A_I(u, p, \tilde{S}(u, p)) \tag{6}
$$

Lorsque l'on écrit l'itération de "défauts corrigés" correspondante, on trouve un problème qui s'écrit
de la manière suivante :

$$
\left\{
\begin{array}{ll}
\rho^f \frac{u_{k+1} - u^n}{\tau} \Big|_{\mathcal{A}} + \rho^f (u^n - w^{n+1}) \cdot \nabla u_{k+1} - \nabla \cdot \sigma^f(u_{k+1}, p_{k+1}) = 0 & \text{dans } \Omega^{f,n+1}, \\[1em]
\nabla \cdot u_{k+1} = 0 & \text{dans } \Omega^{f,n+1}, \\[1em]
\sigma^f(u_{k+1}, p_{k+1}) n^f + \frac{\tau}{\alpha} u_{k+1} = \frac{\tau}{\alpha} \frac{S(u_k, p_k) - d^n}{\tau} + \sigma^f(u_k, p_k) n^f & \text{sur } \Sigma^{n+1}, \\[1em]
\sigma(u_{k+1}, p_{k+1}) n^f = -\bar{p} n^f & \text{sur } \Gamma.
\end{array}
\right.
\tag{7}
$$

On voit là apparaître une relaxation de l'équation de continuité des vitesses en une condition de
Robin purement fluide. Cette itération est équivalente à l'itération Robin-Neumann proposée dans [1],
mais il est intéressant de noter que la formulation (7) ne fait pas intervenir les efforts de la structure
à l'interface (les échanges se font comme dans une itération Dirichlet-Neumann). Comme on le verra
dans les exemples numériques (voir aussi [1]), cette relaxation permet de réduire le nombre d'itérations
nécessaires à la convergence ainsi que la sensibilité de la vitesse de convergence à l'effet de masse
ajoutée.

4.1.1 Le choix du paramètre α

Dans l'algorithmé présenté ci-dessus, on a fait intervenir un réel positif α. En pratique, si on donne
à ce paramètre une valeur trop petite, l'itération ne converge pas. Inversement, si on donne de grandes
valeurs, l'itération converge, mais très lentement. On aimerait donc connaître la valeurs optimale à donner
à ce paramètre en fonction des données physiques et géométriques du problème. On peut également en
savoir plus sur la valeur critique en deçà de laquelle l'itération ne converge plus. L'étude d'un modèle
simplifié réalisée dans [1] permet d'extraire ces valeurs :

$$
\alpha_c = \frac{1}{2} \frac{\tau^2 (\rho^f L - \tanh(\pi R/L) (\beta \tau^2 + \rho^s H))}{L \rho^f (\beta \tau^2 + \rho^s H)},
$$

$$
\alpha_{\text{opt}} = \frac{\tau^2}{\beta \tau^2 + \rho^s H}.
$$

C'est la valeur $\alpha_{\text{opt}}$ qui sera utilisée dans les expérimentations numériques présentées à la fin de l'article.

4.2 Algorithmes 2 : relaxation de la condition d'incompressibilité du fluide

Dans cette section, on propose de changer de préconditionneur pour un modèle inspiré d'un schéma
de projection ([5]). On propose donc de définir l'opérateur fluide simplifié suivant :

$$
\left\{
\begin{array}{l}
\rho^f \left. \frac{u^{n+1} - u^n}{\tau} \right|_{\mathcal{A}} + \rho^f (u^n - w^{n+1}) \cdot \nabla u^{n+1} - 2\mu \nabla \cdot \epsilon(u^{n+1}) = -\nabla p^n \quad \text{dans } \Omega^{f,n+1} \\[1em]
u^{n+1} = \frac{d^n - d^{n-1}}{\tau} \quad \text{sur } \Sigma^{n+1} \qquad (8) \\[1em]
2\mu\epsilon(u^{n+1})n = 0 \quad \text{sur } \Gamma
\end{array}
\right.
$$

$$
\left\{
\begin{array}{ll}
\displaystyle \frac{\tau}{\rho^f} \Delta p^{n+1} = \nabla \cdot u^{n+1} & \text{dans } \Omega^{f,n+1} \\[2ex]
\displaystyle \frac{\tau}{\rho^f} \frac{\partial p^{n+1}}{\partial n} = -\left(\frac{\lambda - 2d^n + d^{n-1}}{\tau}\right) \cdot n^f & \text{sur } \Sigma^{n+1} \\[2ex]
\displaystyle \frac{\tau}{\rho^f} \frac{\partial p^{n+1}}{\partial n} = \gamma(p^{n+1} - \bar{p}) & \text{sur } \Gamma
\end{array}
\right.
\qquad (9)
$$

Ce schéma traite la continuité des vitesses de manière explicite. Remarquons que nous aurions pu
écrire $u^{n+1} = 0$ sur $\Sigma$ sans que cela ne change l'itération préconditionnée. Ce phénomène est une consé-
quence de la structure de l'algorithme des défauts corrigés. Après préconditionnement, nous aboutissons
à l'algorithme suivant :

$$
\left\{
\begin{array}{l}
\rho^f D_{\tau|\mathcal{A}} u_{k+1} + \rho^f (u^n - w^{n+1}) \cdot \nabla u_{k+1} - 2\mu \nabla \cdot \epsilon(u_{k+1}) = -\nabla p_k \quad \text{dans } \Omega^{f,n+1}, \\[1em]
u_{k+1} = \frac{S(u_k, p_k) - d^n}{\tau} \quad \text{sur } \Sigma^{n+1}, \qquad (10) \\[1em]
2\mu\epsilon(u_{k+1})n^f = (p_k - \bar{p})n^f \quad \text{sur } \Gamma,
\end{array}
\right.
$$

$$
\left\{
\begin{array}{ll}
\displaystyle \frac{\tau}{\rho^f} \Delta p_{k+1} = \frac{\tau}{\rho^f} \Delta p_k + \nabla \cdot u_{k+1} & \text{dans } \Omega^{f,n+1}, \\[2ex]
\displaystyle \frac{\tau}{\rho^f} \frac{\partial p_{k+1}}{\partial n} + \frac{\alpha}{\tau} p_{k+1} = \frac{\tau}{\rho^f} \frac{\partial p_k}{\partial n} + \frac{\alpha}{\tau} p_k & \text{sur } \Sigma^{n+1}, \\[2ex]
\displaystyle \frac{\tau}{\rho^f} \frac{\partial p_{k+1}}{\partial n} + \gamma p_{k+1} = \frac{\tau}{\rho^f} \frac{\partial p_k}{\partial n} + \gamma p_k & \text{sur } \Gamma.
\end{array}
\right.
\qquad (11)
$$
Remarquons que nous avons à la fois simplifié le solveur solide, mais aussi modifié la partie fluide, dans le but de permettre une résolution séparée de la vitesse et de la pression du fluide. Dans un premier temps, on résout la première partie pour trouver $u_{k+1}$, ce qui nous permet de résoudre la deuxième partie pour en déduire $p_{k+1}$. Nous verrons dans les exemples numériques que cela permet un gain de temps passé pour une itération, mais aussi un gain en nombre d'itération par rapport à l'algorithme précédent. De plus, il est intéressant de remarquer que le couplage entre la pression et le solide ne se fait qu'au travers du terme $\nabla \cdot u_{k+1}$ puisque la condition sur l'interface $\Sigma$ disparait à la convergence.

**Remarque 1** Dans le préconditionneur, nous avons remplacé la condition de Dirichlet $p = \bar{p}$ sur l'interface d'entrée/sortie par une condition de Robin pour des raisons techniques. Le paramètre $\gamma$ introduit influe peu sur le comportement de l'itération et peut sans problèmes être fixé à 100 et c'est ce qui est fait dans la suite de l'article.

## 4.3 Algorithme 3. Un autre algorithme dans le cas discret

Dans cette partie, on propose de simplifier encore plus le problème initial en négligeant la convection et la viscosité dans le préconditionneur. On peut justifier cette idée en remarquant que l'effet de masse ajoutée provient du caractère incompressible du fluide. On considère donc le préconditionneur suivant, dont la partie qui concerne la vitesse n'a de sens que dans le cadre d'une discrétisation en espace de type éléments finis :

$$ 
\begin{aligned}
\langle \tilde{A}_I(u, p, \lambda), (v, q, \xi) \rangle_{\tilde{V}'_h, \tilde{V}_h} = & \frac{\rho^f}{\tau} \left( \int_{\Omega^{f,n+1}} u \cdot v - \int_{\Omega^{f,n}} u^n \cdot v \right) - \int_{\Omega^{f,n+1}} p \nabla \cdot v + \int_{\Gamma} \bar{p} v \cdot n \\ 
& + \int_{\Sigma^{n+1}} \left( u - \frac{\lambda - d^n}{\tau} \right) \cdot \xi + \int_{\Omega^{f,n+1}} \frac{\tau}{\rho^f} \nabla p \cdot \nabla q \\ 
& + \int_{\Sigma^{n+1}} \left( \frac{\lambda - d^n}{\tau} - u^n \right) \cdot n q + \int_{\Gamma} \frac{\tau}{\gamma \rho^f} (p - \bar{p}) q.
\end{aligned}
$$

Comme dans [9], nous avons une équation d'incompressibilité qui fait intervenir le déplacement à l'interface. C'est précisément ce qui permet d'apporter au solveur fluide une information sur le comportement de l'interface plus précise qu'une simple condition de Dirichlet. Finalement, l'itération préconditionnée consiste à résoudre d'abord le problème suivant sur la pression :

$$ 
\left\{ 
\begin{array}{ll} 
\frac{\tau}{\rho^f} \Delta p_{k+1} = \frac{\tau}{\rho^f} \Delta p_k + \nabla \cdot u_k & \text{dans } \Omega^f, \\ 
\frac{\tau}{\rho^f} \frac{\partial p_{k+1}}{\partial n} + \frac{\alpha}{\tau} p_{k+1} = \frac{\tau}{\rho^f} \frac{\partial p_k}{\partial n} + \frac{\alpha}{\tau} p_k - \frac{\rho^f}{\tau} \left( u_k - \frac{S(u_k, p_k) - d^n}{\tau} \right) \cdot n^f & \text{sur } \Sigma, \\ 
\frac{\tau}{\rho^f} \frac{\partial p_{k+1}}{\partial n} + \frac{\tau \gamma}{\rho^f} p_{k+1} = \frac{\tau}{\rho^f} \frac{\partial p_k}{\partial n} + \frac{\tau \gamma}{\rho^f} p_k & \text{sur } \Gamma. 
\end{array} 
\right. 
\quad (12) 
$$

Et ensuite retrouver la vitesse fluide $u_{k+1}$ en résolvant le problème discret suivant

$$ 
\begin{aligned}
& \frac{\rho^f}{\tau} \left( \int_{\Omega^{f,n+1}} u_{k+1} \cdot v - \int_{\Omega^{f,n}} u^n \cdot v \right) - \rho^f \int_{\Omega^{f,n+1}} (\nabla \cdot w^{n+1}) u_k \cdot v + \rho^f \int_{\Omega^{f,n+1}} (u^n - w^{n+1}) \cdot \nabla u_k \cdot v \\ 
& + \int_{\Omega^{f,n+1}} \sigma(u_k, p_{k+1}) : \epsilon(v) + \int_{\Gamma} \bar{p} v \cdot n^f + \int_{\Sigma^{n+1}} \left[ u_{k+1} - \frac{\alpha(p_{k+1} - p_k)n^f + S(u_k, p_k) - d^n}{\tau} \right] \cdot \xi = 0 
\end{aligned} 
\quad (13) 
$$

Ce nouvel algorithme n'a de sens que dans le cadre d'une discrétisation en espace de type éléments finis, mais il cumule les avantages des deux précédents ;

1. Il permet une résolution séparée des inconnues vitese et pression.

2. Il fait intervenir une relaxation de la condition d'incompressibilité et de continuité des vitesses
## 4.4 Algorithme 4. La compressibilité artificielle

Dans cette section, on propose de préconditionner notre schéma en modifiant l'équation d'incompressibilité comme présenté dans [3] (voir aussi [6]). Encore une fois, le but est de relaxer cette condition, mais d'une autre façon que dans l'algorithme précédent. On introduit alors un réel positif c qui relie linéairement la divergence à la dérivée temporelle de la pression.

$$ \nabla \cdot \mathbf{u}^{n+1} = -c \frac{p^{n+1} - p^n}{\tau} \qquad (14) $$

Nous aboutissons au schéma préconditionné suivant :

$$
\left\{
\begin{array}{ll}
\rho^f \left. \frac{u_{k+1} - u^n}{\tau} \right|_{\mathcal{A}} + \rho^f (u^n - w^{n+1}) \cdot \nabla u_{k+1} - \nabla \cdot \sigma^f (u_{k+1}, p_{k+1}) = 0 & \text{dans } \Omega^{f,n+1}, \\[1em]
\nabla \cdot u^{k+1} = -c \frac{p^{k+1} - p^n}{\tau} & \text{dans } \Omega^{f,n+1} \\[1em]
u^{k+1} = \frac{S(u_k, p_k) - d_n}{\tau} n & \text{sur } \Sigma^{n+1}
\end{array}
\right.
\quad (15)
$$

Ce schéma, présenté dans la littérature (voir, p.e., [3, 6]), présente des similitudes avec les autres schémas. Il s'obtient par préconditionnement du problème initial et fait intervenir un paramètre dont l'ajustement a priori pose problème. Nous présenterons une analyse qui permet de décrire comment la vitesse de convergence varie avec c et quelques valeurs remarquables.

### 4.4.1 Le choix du paramètre c

Comme pour le paramètre $\alpha$ qui intervient dans la condition de Robin des trois premiers algorithmes, on peut définir la valeur critique $c_{crit}$ en deçà de laquelle l'itération ne converge plus et la valeur optimale pour laquelle la convergence est la plus rapide. Comme dans [2, 1], on prend un modèle simplifié qui permet de mener l'analyse. On considère le problème suivant :

$$
\left\{
\begin{array}{ll}
\rho^f \frac{u^{k+1} - u_n}{\tau} + \nabla p^{k+1} = 0 & \text{dans } \Omega^f \\[1em]
\nabla \cdot u^{k+1} = -c \frac{p^{k+1} - p^n}{\tau} \delta_{|\Sigma} & \text{dans } \Omega^f \\[1em]
u^{k+1} = \frac{d^k - d_n}{\tau} n & \text{sur } \Sigma \\[1em]
\rho^s H \frac{d^{k+1} - 2d_n + d_{n-1}}{\tau^2} + \beta d^{k+1} = p^{k+1} & \text{sur } \Sigma
\end{array}
\right.
\qquad (16)
$$

Ce problème reprend les caractéristiques fondamentales du problème complet d'interaction fluide-structure présenté au début de l'article, mais en négligeant quelques termes pour en simplifier l'analyse : la convection, la géométrie mobile, la viscosité. De plus, le solide est ici modélisé par une membrane, ce qui permet d'incorporer son équation d'évolution dans le problème fluide comme condition aux bords. Nous allons résoudre ce problème sur notre géométrie de test 2D composé d'un domaine fluide rectangle de longueur $L$ et de largeur $R$ surmonté d'un solide occupant la frontière $y = R$. Une condition de symétrie est appliquée en $y = 0$, axe de la conduite.

En écrivant l'itération sur la pression, on fait apparaître un opérateur dont le rayon spectral $\sigma$ dépend de $c$. Plus précisément, on a

$$ \sigma = \left| \frac{c - \frac{\tau^2}{\rho^s H + \beta \tau^2}}{c + \frac{\beta \tau^2}{\rho^f} \tanh(\beta R)} \right| $$

Nous obtenons donc, sur le modèle de la partie précédente, les valeurs suivantes :
$$ egin{cases} c_{\text{crit}} = \frac{1}{2} \left( \frac{\beta \tau^2}{\rho^f} \tanh \left( \frac{\pi R}{L} \right) + \frac{1}{\frac{\rho^s H}{\tau^2} + \beta} \right) \\ c_{\text{opt}} = \frac{\tau^2}{\rho^s H + \beta \tau^2} \end{cases} \quad (17) $$

Il est remarquable que l'on retrouve le paramètre optimal de l'Algorithme 1 obtenu dans [1] (voir la section 4.1.1). Cela suggère que l'essentiel du comportement d'un solide fin composé d'une masse et d'une raideur d'ordre zéro est compris dans cette valeur ou bien que les deux algorithmes ne sont que deux manières équivalentes de préconditionner le même problème.

# 5 Résultats numériques et conclusion

Nous avons testé tous les algorithmes à l'aide du logiciel FreeFEM++ dans la géométrie 2D simple présentée dans la section 4.4.1 avec $L = 10\,\text{cm}$ et $R = 0.5\,\text{cm}$, avec cette fois un solide 2D d'épaisseur $H = 0.1\,\text{cm}$. Les conditions aux limites se compose d'une pression nulle en sortie ($x = L$) et une pression $\bar{p} = 10^4\,\text{dyne/cm}^2$ pendant une durée de $5 \times 10^{-3}\,\text{s}$ en entrée ($x = 0$). Les densités fluide et solide sont fixées à $\rho^f = 1\,\text{g/cm}^3$ et $\rho^s = 1.1\,\text{g/cm}^3$. Le fluide a une viscosité $\mu = 0.035\,\text{Po}$ et la structure a pour module d'Young $E = 3 \times 10^8\,\text{Pa}$ et un coefficient de Poisson de $\nu = 0.3$. Nous utilisons un pas de discrétisation temporelle $\tau = 10^{-4}\,\text{s}$ et une discrétisation spatiale par éléments finis de pas $h = 0.1\,\text{cm}$.

Nous avons ajouté à notre test les itérations Dirichlet/Neumann relaxées avec la formule d'Aitken (voir [7] pour une description des différentes méthodes de relaxation dynamique). Cet algorithme est connu pour être sensible à l'effet de masse ajoutée, dicté par la longueur du domaine et le rapport des densités fluide et solide ([2]).

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Nombre moyen d'itérations par pas de temps.</figcaption></figure>

On constate que comparativement aux itérations Dirichlet-Neuman, les algorithmes préconditionnés sont plus efficaces pour des valeurs physiologiques de $L$ et $\rho^s$. De plus, leur efficacité en terme de nombre d'itération par pas de temps est mois sensible à la variation des paramètre physiques, particulièrement pour les algorithmes 3 et 4. Notons toutefois que pour des valeurs très faibles (et non physiologiques) de $\rho^s$, les algorithmes préconditionnés deviennent moins performants.

Si on s'intéresse maintenant au temps CPU écoulé pour calculer 10 pas de temps, on s'attend à ce que les algorithmes qui permettent une résolution séparée des inconnues fluide et solide donnent de meilleurs résultats lorsque le pas du maillage $h$ diminue. C'est effectivement le cas comme le montre le tableau 1.

Ainsi, l'algorithme 4 qui présentait les meilleurs performances en terme de nombre d'itération par pas de temps devient ici moins efficace et les algorithmes 2 et 3 deviennent les meilleurs.

En conclusion, nous avons présenté un modèle d'interaction fluide-structure et différents algorithmes itératifs permettant de s'affranchir de l'effet de masse ajoutée lorsque les paramètres physiques et géométriques varient. Certains ont pour effet de relaxer la condition de continuité des vitesses à l'interface,
<table><caption>TABLE 1 – Temps de calcul pour 20 pas de temps (sans dimension)</caption><thead><tr><th>h</th><th>0.04</th><th>0.05</th><th>0.08</th><th>0.1</th><th>0.2</th></tr></thead><tbody><tr><th>Algorithme 1</th><td>34</td><td>22,6</td><td>9,4</td><td>6,5</td><td>1,8</td></tr><tr><th>Algorithme 2</th><td>21,2</td><td>14,5</td><td>5,5</td><td>3,8</td><td>1</td></tr><tr><th>Algorithme 3</th><td>22,2</td><td>15,4</td><td>6,8</td><td>4,6</td><td>1,4</td></tr><tr><th>Algorithme 4</th><td>49</td><td>32</td><td>12</td><td>7,6</td><td>1,5</td></tr></tbody></table>

d'autres relaxent la condition d'incompressibilité. Tous introduisent un paramètre supplémentaire de relaxation (α et c), dont la valeur à priori n'a rien d'évident. Une valeur de ce paramètre mal ajustée peut considérablement augmenter le nombre d'itération nécessaires à la convergence, voir même compromettre celle-ci. Nous avons étudié la manière dont l'efficacité de l'algorithmes varie avec ce paramètre en étudiant analytiquement un modèle simplifié. Nous avons ensuite confronté les résultats obtenus à la réalité expérimentale pour les valider à posteriori. Finalement, les valeurs qui ressortent de cette étude ont permis de mener une analyse comparative de tous les algorithmes proposés en comparaison avec les itérations Dirichlet-Neumann relaxés par la méthode d'Aitken. Il ressort de cette étude que l'introduction de la relaxation permet d'une part une efficacité supérieure en terme de nombre d'itérations par pas de temps, et d'autre part une relative stabilité de ces performances lorsque les paramètres physiques et géométriques varient. De plus, certain algorithmes permettent une résolution séparée des inconnues vitesse et pression, ce qui diminue le temps de calcul. Nous proposons donc d'une façon unifiée toute une gamme d'algorithmes robustes, donc certains ont déjà été publiés dans la littérature, en apportant des améliorations significatives (découplage, moindre sensibilité à la masse ajoutée).

Références

[1] S. Badia, F. Nobile, C. Vergara. *Fluid-structure partitioned procedures based on Robin transmission conditions*, J. Comp. Phys. 7027-7051, vol 227, 2008.

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

[3] J. Degroote, A. Swillens, P. Bruggeman, R. Haelterman, P. Segers, J. Vierendeels. *Simulation of fluid-structure interaction with the interface artificial compressibility method*, Int. J. Numer. Meth. Biomed. Engng., 276-289, vol 26, 2010.

[4] L. Formaggia, A. Quarteroni, A. Veneziani (Eds.), *Cardiovascular Mathematics. Modeling and simulation of the circulatory system*, Modeling, Simulation and Applications, vol 1, Springer, 2009.

[5] J.-L. Guermond, P. Minev, J. Shen. *An Overview of Projection methods for incompressible flows*, Comput. Methods Appl. Mech. Engrg. 6011-7045, vol 195, 2006.

[6] E. Järvinen, P. Räback, M. Lyly, J.-P. Salenius. *A method for partitioned fluid-structure interaction computation of flow in arteries*, Medical Engineering & Physics, 917-923, vol 30, 2008.

[7] U. Küttler, W.A. Wall. *Fixed-point fluid-structure interaction solvers with dynamic relaxation*, Comput. Mech. 61-72, vol 43, 2008.

[8] P. Le Tallec and J. Mouro. *Fluid structure interaction with large structural displacements*. Comput. Meth. Appl. Mech. Engrg., 3039-3067, vol 190, 2001.

[9] N. Poussineau. *Réduction variationnelle d'un couplage fluide-structure : application à l'hémodynamique*, Thèse de l'Université Paris 6, 2007.

[10] H.J. Stetter. *The defect correction principle and discretization methods*, Numer. Math. 425-443, vol 29, 1978.