# Technique de réduction de modèles appliquée à un problème simple d'interaction fluide-structure

Pascal Schumacher, Alain Combescure, Anthony Gravouil

► To cite this version:

Pascal Schumacher, Alain Combescure, Anthony Gravouil. Technique de réduction de modèles appliquée à un problème simple d'interaction fluide-structure. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717083

HAL Id: hal-01717083

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

Submitted on 25 Feb 2018

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 CC0 1.0 - Universal - International License
# Technique de réduction de modèles appliquée à un problème simple d'interaction fluide-structure

Pascal SCHUMACHER<sup>1,2,*</sup>, Alain COMBESCURE<sup>1</sup>, Anthony GRAVOUIL<sup>1</sup>

## Résumé —

Les temps de résolution de problèmes couplés fluide/structure restent à l'heure actuelle importants, C'est pourquoi une technique de réduction de modèle, du type projection sur une base réduite, est experimentée. Afin d'étudier la faisabilité d'une telle méthode, un modèle couplé unidimensionnel a été réalisé avec un domaine fluide discrétisé par des volumes finis et un domaine solide par des éléments finis. On impose la continuité des vitesses entre les deux milieux de telle manière à préserver une énergie d'interface nulle. Cette méthode permet de rendre plus stable le couplage entre les deux domaines et de garder l'ordre de convergence de chacun des intégrateurs temporels. Le calcul de la base réduite est effectué à posteriori grâce à une décomposition orthogonale aux valeurs propres (POD).

**Mots clés** — Interaction Fluide/Structure, Réduction de modèle, POD.

## 1 Introduction

Afin d'accélérer les temps de calcul, de nouvelles méthodes basées sur la réduction de modèle ont été récemment mises au point (e.g. [1]). L'idée centrale est de rechercher la solution dans un sous espace réduit et d'adapter ce sous espace au fur et à mesure des itérations. Afin d'étudier une possible application de ces méthodes de réduction *a priori* à un problème couplé fluide/structure, une étude de réduction de modèle *a posteriori* a été effectuée sur un problème unidimensionnel.

Le problème considéré est une poutre encastrée en contact avec un fluide compressible. la structure est discretisée à l'aide de la méthode des éléments finis et d'un schéma temporel de Newmark. Les équations du milieu fluide, sous leur forme conservative, sont traitées par la méthode des volumes finis.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – modèle 1D</figcaption></figure>

Afin de relier les deux domaines, une approche duale a été adoptée. La continuité des vitesses à l'interface est imposée sur l'intervalle de temps $[t_n, t_{n+1}]$ de manière à ce que l'énergie d'interface soit exactement nulle (c.f. [2],[3]).

Pour étudier l'efficacité des techniques de réduction de modèle, une base réduite a été calculée *a posteriori* grâce à une décomposition orthogonale aux valeurs propres (POD) car elle fournit par définition la meilleure base de dimension réduite. Ce calcul a été réalisé par une décomposition en valeurs singulières (SVD) de la matrice solution espace temps.

<footnote n="1"><sup>1</sup> Université de Lyon, INSA-Lyon, LaMCoS, CNRS UMR 5259, F-69621, Villeurbanne (France),
{pascal.schumacher, alain.combescure, anthony.gravouil} @insa-lyon.fr</footnote>

<footnote n="2"><sup>2</sup> ANSYS France, Villeurbanne (France),
pascal.schumacher@ansys.com</footnote>

<footnote>* Auteur correspondant</footnote>
# 2 Description du problème solide

## 2.1 Mise en équation

La partie structure est modélisée par une poutre en traction-compression de longueur L et de section S. Elle est encastrée à sa base ($x = 0$). les caractéristiques du materiau sont le module d'Young E et sa masse volumique $\rho_s$. On impose à l'extrémité droite une force $\lambda$. Soit $u$ le champ de déplacement inconnu de la poutre, le problème peut se mettre sous la formulation forte suivante :

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Schéma du modèle structure 1D étudié</figcaption></figure>

$$ egin{cases} 
\rho_s S \frac{\partial^2 u}{\partial t^2} - E S \frac{\partial^2 u}{\partial x^2} = 0 & \text{sur } \Omega_s = ]0, L[ \\ 
ES \frac{\partial u}{\partial x} = -\lambda & \text{sur } \partial\Omega_s^N = \{L\} \\ 
u = 0 & \text{sur } \partial\Omega_s^D = \{0\} 
\end{cases} 
\quad (1) $$

On exprime cette équation sous la formulation faible : il suffit de multiplier par une fonction test $v$ cinématiquement admissible à 0 et de réaliser une intégration par partie pour faire ressortir les conditions aux limites.

$$ \forall v \text{CA0, } \int_{\Omega_s} \rho_s S v \cdot \frac{\partial^2 u}{\partial t^2} d\Omega + \int_{\Omega_s} E S \frac{\partial v}{\partial x} \cdot \frac{\partial u}{\partial x} d\Omega + v(L)\lambda = 0 \quad (2) $$

## 2.2 Discretisation

la structure est discretisée par la méthode des éléments finis. le champ de déplacement solution est approximé par une combinaison linaire finie de fonction que l'on résume par l'équation suivante $u(x,t) = N(x)u(t)$. La même approximation est faite pour la fonction test.

$$ \forall v \quad v^T \cdot \left( \left[ \int_{\Omega_s} \rho_s S N^T N d\Omega \right] \ddot{u} + \left[ \int_{\Omega_s} E S \frac{\partial N^T}{\partial x} \frac{\partial N}{\partial x} d\Omega \right] u + L_s^T \lambda \right) = 0 $$

On obtient finalement l'équation semi-discretisée bien connue que l'on exprime plus succintement à l'aide de la matrice de masse $M$ et de la matrice de raideur $K$.

$$ M\ddot{u} + Ku + L_s^T \lambda = 0 \quad (3) $$

La discretisation temporelle a été effectuée avec un schéma classique de Newmark :

$$ \dot{u}^{n+1} = \dot{u}^n + (1-\gamma)\Delta t\ddot{u}^n + \gamma\Delta t\dot{u}^{n+1} \quad (4) $$

$$ u^{n+1} = u^n + \Delta t \dot{u}^n + \left( \frac{1}{2} - \beta \right) \Delta t^2 \ddot{u}^n + \beta \Delta t^2 \dot{u}^{n+1} \quad (5) $$

# 3 Description du problème fluide

## 3.1 Mise en équation

La partie fluide est un gaz dans une cavité de longueur L et de section S. Le problème est donc modelisé par un fluide accoustique de célérité $c$ et de masse volumique moyenne $\rho_0$. une pression est appliquée à l'extrémité gauche tandis que l'on impose une vitesse sur l'extrémité droite.
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Schéma du modèle fluide 1D étudié</figcaption></figure>

$$
\left\{
\begin{array}{ll}
\displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial (\rho_0 v)}{\partial x} = 0 & \text{sur } \Omega_f = ]L, 2L[ \\
\displaystyle \frac{\partial (\rho_0 v)}{\partial t} + c^2 \frac{\partial \rho}{\partial x} = 0 & \text{sur } \Omega_f \\
\rho = \frac{\lambda}{Sc^2} & \text{sur } \partial \Omega_f^N = \{L\} \\
v = v_{imp} & \text{sur } \partial \Omega_f^D = \{2L\}
\end{array}
\right.
\qquad (6)
$$

On peut reformuler ce système d’équations au dérivées partielles en posant un vecteur d’état $U = [\rho, \rho_0 v]^T$ contenant les variables conservative du problème.

$$
\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} = 0 \qquad (7)
$$

$$
\text{avec } U = \begin{pmatrix} \rho \\ \rho_0 v \end{pmatrix}, \quad F = AU, \quad A = \begin{pmatrix} 0 & 1 \\ c^2 & 0 \end{pmatrix},
$$

3.2 Discretisation

le domaine fluide est discretisé par la méthode des volumes finis. Elle consiste à effectuer un bilan
de flux sur chaque volume de contrôle Kᵢ. Pour cela, il suffit d'intégrer l'équation (7) sur la cellule de
contrôle Kᵢ.

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Exemple de volume de contrôle K<sub>i</sub></figcaption></figure>

$$
\frac{d}{dt} \int_{K_i} U d\Omega + \oint_{\partial K_i} F dS = 0
$$

On approxime l’intégrale de volume à l’aide sa valeur $U_i$ au centre de la cellule. Il ne reste plus qu’à connaître les flux sur $\partial K_{i+1/2}$ et $\partial K_{i-1/2}$.

$$
\Delta x_i \frac{dU_i}{dt} + F_{i+\frac{1}{2}} - F_{i-\frac{1}{2}} = 0 \quad (8)
$$

Nous utiliserons par la suite deux types de schéma courant : différence centrée (Central Differencing) ou différence décentrée amont (Upwind Differencing). La première méthode consiste en une moyenne des valeurs au centre des cellules adjacentes :

$$
F_{i+\frac{1}{2}} = \frac{1}{2} A (U_i + U_{i+1}) \qquad (9)
$$

La deuxième méthode décompose le flux en deux parties : un flux entrant et un flux sortant. l’interpolation de la variable $U$, du centre des cellules au bord, est différente selon le flux sortant ou entrant.

$$
F_{i+\frac{1}{2}} = A^{+}U_{i} + A^{-}U_{i+1} \quad (10)
$$

$$
\text{avec } A^{+} = \frac{1}{2} \begin{pmatrix} c & 1 \\ c^2 & c \end{pmatrix} \quad \text{et} \quad A^{-} = \frac{1}{2} \begin{pmatrix} -c & 1 \\ c^2 & -c \end{pmatrix}
$$
il ne reste plus qu’à introduire un schéma de discretisation temporelle. On intègre l’équation (8) sur l’intervalle $[t_n, t_{n+1}]$ :

$$ \Delta x_i (U_i^{n+1} - U_i^n) + \int_{t_n}^{t_{n+1}} (F_{i+\frac{1}{2}} - F_{i-\frac{1}{2}}) (t) dt = 0 \quad (11) $$

$$ \text{En posant} \quad \int_{t_n}^{t_{n+1}} F(t) dt = \Delta t(1-\theta)F^n + \Delta t\theta F^{n+1} \quad (12) $$

Les divers valeurs de $\theta$ définissent les schémas d’intégration : $\theta = 0$ Euler explicite, $\theta = 1/2$ Crank-Nicholson, $\theta = 1$ Euler implicite.

# 4 algorithme de couplage conservatif en énergie

## 4.1 Energie d'interface nulle

Afin de relier le domaine fluide et le domaine structure, la continuité des vitesses à l'interface a été recherchée. La démarche suivie pour imposer une énergie d'interface nulle est dans la continuité de [2] et [3]. Elle repose sur deux idées : prendre l'inconnue d'interface moyen $\bar{\lambda}_n$ sur l'intervalle de temps $[t_n, t_{n+1}]$ au lieu de $\lambda(t_{n+1})$ ; assurer la continuité des vitesse moyenne sur l'intervalle de temps $[t_n, t_{n+1}]$.

<figure><img src="image_6.png" /><figcaption>Fig. 5 – Discretisation du modèle complet</figcaption></figure>

A partir du bilan d'énergie discretisé pour le schéma de Newmark, on remarque que travail des efforts à l'interface fait intervenir le paramètre d'intégration $\gamma$. Ce n'est plus le cas si l'on pose successivement :

$$ \bar{\lambda}_n = \gamma \lambda_{n+1}^s + (1-\gamma) \lambda_n^s \quad (\text{structure}) \quad (13) $$

$$ \bar{\lambda}_n = \theta \lambda_{n+1}^f + (1 - \theta) \lambda_n^f \quad (\text{fluide}) \quad (14) $$

D'après [3], en considérant que l'effort d'interface est constant sur l'intervalle de temps $[t_n, t_{n+1}]$, il est possible de conclure une énergie d'interface nulle si l'on a assuré la continuité des vitesses d'interface intégrées sur le pas de temps ($v_s^I = \dot{u}_N$ et $v_f^I = v_{1/2}$).

$$ \begin{gathered} \Delta W_I^n = \int_{t_n}^{t_{n+1}} \lambda (v_s^I - v_f^I) dt = \bar{\lambda}_n \int_{t_n}^{t_{n+1}} (v_s^I - v_f^I) dt \\ \Delta W_I^n = 0 \quad \text{si} \quad \int_{t_n}^{t_{n+1}} (v_s^I - v_f^I) dt = 0 \end{gathered} \quad (15) $$

## 4.2 Algorithme de couplage

Si l'on réécrit le problème discretisé complet, un système de la forme suivante est obtenu :

$$ \begin{pmatrix} \mathbf{A}_s & \mathbf{0} & \mathbf{L}_s^T \\ \mathbf{0} & \mathbf{A}_f & \mathbf{L}_f^T \\ \mathbf{L}_s & \mathbf{L}_f & \mathbf{0} \end{pmatrix} \begin{pmatrix} \mathbf{x}_s \\ \mathbf{x}_f \\ \boldsymbol{\lambda} \end{pmatrix} = \begin{pmatrix} \mathbf{b}_s \\ \mathbf{b}_f \\ \mathbf{b}_i \end{pmatrix} \quad (16) $$

Afin de garder une mise en oeuvre partitionnée, le problème est résolu en trois étapes :

- Résolution des problèmes sans liaison :

$$ \mathbf{A}_s \mathbf{x}_s^{sl} = \mathbf{b}_s $$

$$ \mathbf{A}_f \mathbf{x}_f^{sl} = \mathbf{b}_f $$
– Résolution du problème d'interface :

$$ H\lambda = L_s x_s^{sl} + L_f x_f^{sl} - b_i \quad (17) $$

$$ \text{avec } H = L_s A_s^{-1} L_s^T + L_f A_f^{-1} L_f^T \quad (18) $$

– Résolution des problèmes avec liaison :

$$ A_s x_s^{al} + L_s \lambda = \mathbf{0} $$

$$ A_f x_f^{al} + L_f \lambda = \mathbf{0} $$

La solution finale est :

$$ x_s = x_s^{sl} + x_s^{al} $$

$$ x_f = x_f^{sl} + x_f^{al} $$

La seule difficulté de cette méthode est d'obtenir l'opérateur d'interface **H**.

# 5 Réduction de modèle

## 5.1 Projection sur une base réduite

Pour réduire le nombre d'inconnues d'un système, il peut être interessant de rechercher la solution d'un système linéaire (19) dans un sous espace **Q**.

$$ A\mathbf{x} = \mathbf{b} \quad (19) $$

Pour obtenir un système d'équation nous permettant de trouver une solution dans le sous espace **Q**, il faut effectuer une projection de l'équation (19). La plus courante est celle de Galerkin. Si la matrice **A** est symétrique définie positive, elle correspond à la minimisation d'une énergie dans le sous-espace **Q** :

$$ \min_{\mathbf{x} \in \mathbf{Q}} \left( \frac{1}{2} \mathbf{x}^T \mathbf{A} \mathbf{x} - \mathbf{x}^T \mathbf{b} \right) \implies [\mathbf{Q}^T \mathbf{A} \mathbf{Q}] \cdot \tilde{\mathbf{x}} = \mathbf{Q}^T \mathbf{b} \quad (G) \quad (20) $$

Il est aussi possible d'utiliser d'autres types de projection. on peut essayer de minimiser le résidu de l'équation. On obtient la solution dans **Q** au sens des moindres carré.

$$ \min_{\mathbf{x} \in \mathbf{Q}} \| \mathbf{A}\mathbf{x} - \mathbf{b} \|_2^2 \implies [(\mathbf{A}\mathbf{Q})^T (\mathbf{A}\mathbf{Q})] \cdot \tilde{\mathbf{x}} = (\mathbf{A}\mathbf{Q})^T \mathbf{b} \quad (MR) \quad (21) $$

Dans un problème d'interaction fluide/structure, si nous disposons de bases réduites ($Q_s, Q_f$) dans chaque domaine, il est possible d'estimer l'opérateur d'interface à l'aide des modèles réduits. Voici un exemple avec une projection de Galerkin dans chaque domaine :

$$ H = \left[ L_s Q_s (Q_s^T A_s Q_s)^{-1} L_s^T Q_s^T + L_f Q_f (Q_f^T A_f Q_f)^{-1} L_f^T Q_f^T \right] \quad (G) \quad (22) $$

## 5.2 Choix de la base réduite

Toute la difficulté réside dans le choix d'une base réduite **Q**. Il existe plusieurs possibilités notamment une base formée à partir des modes propres du système. Afin d'étudier la réductibilité du problème ainsi que la qualité des projections sur le sous-espace réduit, le choix s'est porté vers une décomposition orthogonale aux valeurs propres (POD). Connaissant la solution à différents instants :

$$ \mathbf{X} = \begin{bmatrix} x_1(t_1) & \cdots & x_1(t_m) \\ \vdots & \ddots & \vdots \\ x_N(t_1) & \cdots & x_N(t_m) \end{bmatrix}, \quad (23) $$
La POD consiste à trouver une suite de vecteur maximisant la colinérité avec l'ensemble de vecteur représenté par **X**. Cet outil statistique permet de trouver les meilleurs sous-espaces réduits pour représenter approximativement **X**. Les vecteurs POD peuvent être calculer en effectuant une décomposition en valeurs singulière (SVD) de la matrice **X**.

$$ \mathbf{X} = \mathbf{USV} = \sum_{k=1}^{\min(N,m)} s_k \mathbf{u}_k \mathbf{v}_k^T, \qquad (24) $$

On peut alors récupérer les premiers vecteurs issus de la POD pour former une base réduite : **Q** = [u₁, ..., uₚ]. On parle alors de réduction de modèle *a posteriori* car il est nécessaire de connaître déjà la solution pour après réduire la dimension du modèle.

Les valeurs singulières $s_k$ ont une importance pour nous décrire la réductibilité de notre problème. En effet s'il l'on sélectionne que les r premiers termes de la somme de (24), la matrice ainsi formée $\tilde{\mathbf{X}}$ représente la meilleure solution dans le sous-espace **Q** pour la norme matricielle de Frobenius (théorème d'Eckart-Young).

$$ \tilde{\mathbf{X}} = \sum_{k=1}^{r} s_k \mathbf{u}_k \mathbf{v}_k^T, \quad r \ll N \\ \| \tilde{\mathbf{X}} - \mathbf{X} \|_F^2 = \sum_{k=r+1}^{\min(N,m)} s_k^2 \qquad (25) $$

Les valeurs singulières nous permettent d'estimer une erreur minimale qui sera faite avec le modèle réduit dans chaque sous domaine. Grâce à la formule (25), plus la suite de valeurs singulières décroît rapidement, plus le problème est réductible.

# 6 Résultats numériques

Pour le problème un 1D, les paramètres materiaux ont été pris pour avoir une célérité identique dans chaque domaine (1 m.s⁻¹). chaque domaine est d'une longueur égal à 1 m. la poutre est discretisée avec 1000 éléments et la cavité fluide avec 500 volumes de contrôle (1000 inconnues dans chaque domaine). A l'instant initial tous les champs sont nuls. On impose une vitesse sur le bord de la cavité. Par la suite, pour le schéma volume fini, le flux a été estimé par un schéma centré et l'intégration temporelle par un schéma de crank-nicholson (θ = 1/2).

<figure><img src="image_7.png" /><figcaption>Fig. 6 – Modèle</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 – Déplacement imposé</figcaption></figure>
Indépendamment du choix des paramètres d'intégration du schéma de Newmark, l'énergie d'interface a été observée nulle comme le montre la figure 8.

<figure><img src="image_9.png" /><figcaption>Fig. 8 – Bilan d'énergie</figcaption></figure>

La comparaison entre l'erreur minimale défini avec (25) et l'erreur réellement effectuée (modèle réduit POD / modèle non réduit) permet de mettre en relief la qualité des projections sur chaque domaine. Les figures 9 et 10 affichent l'erreur lorsque l'on réduit un seul domaine. La figure 11 montre l'erreur sur chaque domaine lorsque l'on réduit les deux domaines.

<figure><img src="image_10.png" /><figcaption>Fig. 9 – Comparaison entre l'erreur minimale et l'erreur réelle pour différente taille de base réduite POD (domaine solide à gauche, domaine fluide à droite)</figcaption></figure>

On peut remarquer que la projection est moins efficace du coté fluide à cause de la présence de matrice non-symétrique. Dans ce cas de figure, une projection au sens des moindres carré (minimisation du résidu) se révèle plus interessante. En observant plus précisément l'erreur commise par le modèle réduit (fluide et structure) sur la figure 11, il est possible d'avoir un modèle précis à $10^{-2}$ en ne retenant que quelques dizaines de vecteurs sur mille. Le problème d'interface semble ne pas être excessivement influencé par la réduction de modèle.
<figure><img src="image_11.png" /><figcaption>Fig. 10 – Comparaison entre l’erreur minimale et l’erreur réelle pour différente taille de base réduite POD (domaine solide à gauche, domaine fluide à droite)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11 – Comparaison entre l’erreur minimale et l’erreur réelle pour différente taille de base réduite POD (domaine solide à gauche, domaine fluide à droite)</figcaption></figure>

Références

[1] D. Ryckelynck. *A priori hyperréduction method : an adaptative approach*, JCP, 346-366, 2005.

[2] N. Mahjoubi, A. Gravouil, A. Combescure, N. Greffet. *A monolithic energy conserving method to couple heterogeneous time integrators with incompatible time steps in structural dynamic*, CMAME, 1069-1086, 2011.

[3] Z. Li, A. Combescure, F. Leboeuf. *Coupling a weak solid formulation with a strong fluid one for transient fluid-structure interaction*, IJNMF, in press 2013.

[4] C. Hirsch *Numerical Simulation of Internal and External. Volume 1*, John Wiley and Sons, 1988.

[5] C. Hirsch *Numerical Simulation of Internal and External. Volume 2*, John Wiley and Sons, 1990.