r4.02.02 Éléments pour le couplage interaction fluide-structure linéaire avec fluide inerte

Résumé :

Les éléments décrits ici permettent de réaliser des calculs des fréquences et de modes propres d’une structure couplée à un fluide. Ils permettent également le calcul de réponse acoustique.

Après la formulation du problème de couplage fluide-structure, ce document décrit la démarche suivie pour implémenter les nouveaux éléments finis.

Le couplage vibro-acoustique

Présentation

Soit une structure élastique définie dans un domaine \({\Omega}_{s}\) qui vibre en présence d’un fluide parfait, non pesant, compressible, en évolution isentropique défini dans un domaine \({\Omega}_{f}\) . \({\Sigma}_{f}\) étant le bord du domaine fluide \({\Omega}_{f}\) et \({\Sigma}_{s}\) étant le bord du domaine solide \({\Omega}_{s}\) , on désigne par \(\Sigma ={\Sigma}_{f}\cap {\Sigma}_{S}\) leur surface commune. On note \({n}_{s}\) la normale extérieure au domaine fluide \({\Omega}_{f}\) .

À un instant donné, l’état du fluide est défini par son champ de pression (scalaire) \(p\) et celui de la structure par son champ de déplacement (vectoriel) \({u}_{s}\) . On considère que le système couplé est soumis à des petites perturbations autour de son état d’équilibre où le fluide et la structure sont au repos. Le problème d’interaction fluide-structure consiste alors à résoudre simultanément deux problèmes:

• l’un dans la structure soumise, sur \(\Sigma\) , à un champ de pression \(p\) imposé par le fluide;

• l’autre dans le fluide soumis à un champ de déplacement \(u\) de la paroi \(\Sigma\) .

Remarque : si le sens de la normale fluide vers structure ou structure vers fluide n’a pas d’importance (c’est une convention), il est néanmoins impératif que cettenormale soit toujours orientée dans le mêmesens. Il est fortement conseillé de garder la convention d’orientation de la structure vers le fluide pour toutes les modélisations d’interface fluide-structureafin de garder une cohérence avec les autres chargements sur la structure (voir en particulier le cas des chargements de type ONDE_PLANE).

Formulation du problème vibro-acoustique

Description de la structure

Hypothèse : la structure est homogène et obéit aux lois de l’élasticité linéaire. Compte tenu de cette hypothèse, on peut écrire les différentes équations suivantes gouvernant l’état de la structure [bib2].

Équation de conservation de la quantité de mouvement

L’équation de conservation de la quantité de mouvement s’écrit, en l’absence de forces volumiques autres que les forces d’inertie, comme suit:

(918)\[{\rho}_{s}{\ddot{u}}_{s}-div({\sigma}_{s}({u}_{s}))=0\text{dans}{\Omega}_{s}\]

Avec \({\rho}_{s}\) la masse volumique, \({u}_{s}\) les déplacements et \({\sigma}_{s}\) le tenseur des contraintes.

Relation de compatibilité

On établit la relation de compatibilité classique sur le tenseur des déformations:

(919)\[{\epsilon}_{s}=\frac{1}{2}\left(\nabla {u}_{s}+\nabla {u}_{s}^{T}\right)\]

Où \({\epsilon}_{s}\) est le tenseur des (petites) déformations.

Loi de comportement en élasticité linéaire isotrope

On suppose que le solide est élastique linéaire, donc:

(920)\[{\sigma}_{s}({u}_{s})=C:{\epsilon}_{s}({u}_{s})\]

\(C\) étant le tenseur d’élasticité, symétrique. Pour le cas isotrope, on peut aussi écrire:

(921)\[{\sigma}_{s}({u}_{s})=\lambda \text{tr}({\epsilon}_{s}({u}_{s}))I+2\mu {\epsilon}_{s}({u}_{s})\]

Avec \(\lambda ` et :math:\)mu ` les coefficients de Lamé.

Conditions limites

Il y a des conditions limites sur la surface de la structure telles que:

(922)\[\begin{split}\lbrace \begin{array}{c}{\sigma}_{s}({u}_{s})\cdot {n}_{s}={f}_{s}\text{sur}\partial {\Omega}_{N}\\ {u}_{s}=\overline{u}\text{sur}\partial {\Omega}_{U}\end{array}\end{split}\]

Où \(\partial {\Omega}_{N}\) est la partie de la surface de la structure sur laquelle les conditions en efforts (surfaciques) \({f}_{s}\) s’appliquent et \(\partial {\Omega}_{U}\) est la partie de la surface de la structure sur laquelle les conditions en déplacement \(\overline{u}\) s’appliquent.

Description du fluide

Hypothèse: le fluide obéit aux lois de l’acoustique linéaire.

Équation de conservation de la quantité de mouvement

L’équation de conservation de la quantité de mouvement s’écrit, en l’absence de sources:

(923)\[{\rho}_{0}{\ddot{u}}_{f}-div({\sigma}_{f})=0\text{dans}{\Omega}_{f}\]

Avec \({\sigma}_{f}\) le tenseur des contraintes du fluide, \({\rho}_{0}\) la masse volumique du fluide à l’état de repos et \({u}_{f}\) le champ de déplacement du fluide.

Équation de conservation de la masse

Au premier ordre et en l’absence de sources acoustiques, l’équation de conservation de la masse s’exprime par la relation:

(924)\[{\dot{\rho}}_{f}+{\rho}_{0}div\left({\dot{u}}_{f}\right)=0\]

Loi de comportement

Hypothèse: le fluide obéit aux lois de l’acoustique linéaire, donc:

(925)\[{\sigma}_{f}=-pI\]

Le fluide est supposé en évolution barotrope (la pression \(p\) est, pour le fluide donné, une fonction connue de la seule masse volumique):

(926)\[p={\rho}_{f}{c}_{0}^{2}\]

où \({c}_{0}\) est la célérité du son dans le fluide au repos.

Formulation en pression

On doit maintenant combiner ces équations pour obtenir l’équation d’état du fluide, il y a plusieurs manières de le faire. Si on prend la pression \(p\) du fluide comme variable d’état, en combinant () et (), on obtient:

(927)\[{\rho}_{0}div({\ddot{u}}_{f})+\Delta p=0\text{dans}{\Omega}_{f}\]

En dérivant par rapport au temps l’équation () et en la combinant à ():

(928)\[\frac{1}{{c}_{0}^{2}}\ddot{p}+{\rho}_{0}div\left({\ddot{u}}_{f}\right)=0\]

Si on injecte () dans (), on obtient l’équation de d’Alembert:

(929)\[\frac{1}{{c}_{0}^{2}}\ddot{p}-\Delta p=0\text{dans}{\Omega}_{f}\]

Formulation à deux variables

On associe à la variable de pression \(p\) une variable supplémentaire appelée potentiel de déplacement du fluide et notée :math:`varphi ` , on fait l’hypothèse d’absence de terme rotationnel dans les mouvements du fluide tel que:

(930)\[{u}_{f}=\nabla \varphi\]

Voir [bib1], [bib4], [bib5].

En prenant l’équation de conservation de la quantité de mouvement () et en y injectant () après dérivation, qui est une ré-écriture de l’équation de Bernoulli avec \({\dot{u}}_{f}^{2}=0\) (hypothèse de linéarité):

(931)\[\nabla \left(p+{\rho}_{0}\ddot{\varphi}\right)=0\]

Ce qui nous donne une relation directe entre la pression et le potentiel du déplacement du fluide:

(932)\[\frac{1}{{\rho}_{0}{c}_{0}^{2}}p+\frac{1}{{c}_{0}^{2}}\ddot{\varphi}=0\]

D’où la nouvelle forme d’équation de propagation des ondes:

(933)\[\frac{1}{{c}_{0}^{2}}\ddot{p}+{\rho}_{0}\Delta \ddot{\varphi}=0\]

Formulation en potentiel de vitesse

On peut également introduire le potentiel de vitesse du fluide, noté \(\psi\) et tel que:

(934)\[{\dot{u}}_{f}=\nabla \psi\]

Ce qui nous donne comme relation avec la pression \(p\) du fluide:

(935)\[p=-{\rho}_{0}\dot{\psi}\]

Cette dernière formulation est plus naturelle lorsque l’on doit imposer une condition limite sur le bord en vitesse, comme une onde incidente.

Description de l’interaction fluide-structure

A l’interface fluide-structure \(\Sigma\) , le fluide étant non visqueux, il n’adhère pas à la paroi. On écrit donc la continuité des contraintes normales:

(936)\[{\sigma}_{s}\cdot {n}_{s}=-p{n}_{s}\text{sur}\Sigma\]

Et la continuité des déplacements normaux:

(937)\[{u}_{s}\cdot {n}_{s}={u}_{f}\cdot {n}_{s}\text{sur}\Sigma\]

Comme le fluide est considéré initialement au repos, la continuité des déplacements normaux est équivalente à la continuité des vitesses normales:

(938)\[{\dot{u}}_{s}\cdot {n}_{s}={\dot{u}}_{f}\cdot {n}_{s}\text{sur}\Sigma\]

Et la continuité des accélérations normales:

(939)\[{\ddot{u}}_{s}\cdot {n}_{s}={\ddot{u}}_{f}\cdot {n}_{s}\text{sur}\Sigma\]

En appliquant () dans ():

(940)\[{\ddot{u}}_{f}=-\frac{1}{{\rho}_{0}}\nabla p\]

L’équation de continuité des accélérations normales () s’écrit:

(941)\[{\ddot{u}}_{s}\cdot {n}_{s}=-\frac{1}{{\rho}_{0}}\nabla p\cdot {n}_{s}\text{sur}\Sigma\]

Pour le cas où l’on exprime l’équation d’état du fluide par la pression et le potentiel de déplacement du fluide, les conditions IFS seront:

(942)\[ \begin{align}\begin{aligned}{\sigma}_{s}\cdot {n}_{s}={\rho}_{0}\ddot{\varphi}{n}_{s}\text{sur}\Sigma\\{\ddot{u}}_{s}\cdot {n}_{s}=\nabla \ddot{\varphi}\cdot {n}_{s}\text{sur}\Sigma\end{aligned}\end{align} \]

Et si on utilise le potentiel de vitesse du fluide:

(943)\[ \begin{align}\begin{aligned}{\sigma}_{s}\cdot {n}_{s}={\rho}_{0}\dot{\psi}{n}_{s}\text{sur}\Sigma\\{\ddot{u}}_{s}\cdot {n}_{s}=\nabla \varphi \cdot {n}_{s}\text{sur}\Sigma\end{aligned}\end{align} \]

Ces écritures permettent de constater que le sens de la normale \({n}_{s}\) est bien uniquement une affaire de convention d’écriture. Quelle que soit cette orientation, les équations de continuité auront strictement la même forme.

Condition de rayonnement à l’interface extérieure

Dans le cas où la structure est entourée par le fluide, la propagation des ondes dans le fluide doit respecter les conditions de rayonnement à la frontière extérieure de ce milieu fluide \(\partial {\Sigma}_{R}\) , cette condition appelée condition de Sommerfeld stipule que les ondes ne sont pas réfléchies par la surface extérieure du fluide. La condition de rayonnement s’écrit:

(944)\[\underset{r\to \infty }{\lim}{r}^{d-1/2}\left(\nabla p\cdot {n}_{f}+\frac{\dot{p}}{{c}_{0}}\right)=0\]

Cette condition étant asymptotique, son application à une modélisation par éléments finis n’est pas directe. On utilisera une approximation de la condition, notée BGT (comme le nom des auteurs, voir [Bib6]). L’approximationd’ordre zéros’écrit:

(945)\[ \begin{align}\begin{aligned}\nabla p\cdot {n}_{f}+\frac{{\rho}_{0}}{{Z}_{C}}\dot{p}=0\text{sur}\partial {\Sigma}_{R}\\\textrm{avec}\\{Z}_{C}={\rho}_{0}{c}_{0}\end{aligned}\end{align} \]

Cette condition est exacte pour une onde plane. Pour approcher asymptotiquement le comportement à l’infini d’une onde sphérique, pour un rayon \(R\) donné, on utilise une approximation d’ordre 1:

(946)\[ \begin{align}\begin{aligned}\nabla p\cdot {n}_{f}+{\rho}_{0}\left(\frac{1}{{Z}_{C}}\dot{p}+\frac{1}{{Z}_{R}}p\right)=0\text{sur}\partial {\Sigma}_{R}\\\textrm{avec}\\{Z}_{C}={\rho}_{0}{c}_{0}` et :math:`{Z}_{R}={\rho}_{0}R\end{aligned}\end{align} \]

Décomposition de la vitesse du fluide

La vitesse du fluide \({\dot{u}}_{f}\) qui peut s’écrire :math:`{dot{u}}_{f}=nabla psi ` peut être décomposée en trois composantes:

(947)\[\psi ={\psi}^{\mathit{inc}}+{\psi}^{\mathit{rad}}+{\psi}^{\mathit{ref}}\]

Le premier terme est le terme d’incidence, le deuxième celui de radiation et le troisième celui de réflexion. Cette décomposition correspond à la même sur la pression du fluide:

(948)\[p={p}^{\mathit{inc}}+{p}^{\mathit{rad}}+{p}^{\mathit{ref}}\]

Formulation du problème couplé

En définitive, la formulation du problème de vibro-acoustique en termes de déplacements pour la structure et de pression dans le fluide conduit aux équations du problème harmonique (\(P\) ) dans le milieu solide:

(949)\[{C}_{ijkl}\cdot {u}_{k,\mathit{lj}}+{\omega}^{2}{\rho}_{S}{u}_{i}=0\text{dans}{\Omega}_{S}\]

L’équation de propagation des ondes dans le milieu fluide:

(950)\[\Delta p+{k}^{2}p=0\text{dans}{\Omega}_{f}\]

Avec les deux équations de couplage fluide-structure. En contraintes:

(951)\[{C}_{ijkl}\cdot {u}_{k,l}\cdot {n}_{i}=-p{\delta}_{ij}\cdot {n}_{i}\text{sur}\Sigma\]

Et en vitesses:

(952)\[{u}_{i}\cdot {n}_{i}=\frac{1}{{\rho}_{0}{\omega}^{2}}\frac{\partial p}{\partial n}\text{sur}\Sigma\]

En utilisant la décomposition (), on a la condition () à l’interface qui s’écrit:

(953)\[\begin{split}\lbrace \begin{array}{c}\nabla {\psi}^{\mathit{ref}}\cdot {n}_{s}=-{v}^{\mathit{inc}}\cdot {n}_{s}\\ \nabla {\psi}^{\mathit{rad}}\cdot {n}_{s}={\dot{u}}_{s}\cdot {n}_{s}\end{array}\text{sur}\Sigma\end{split}\]

Modèles éléments finis du problème couplé

On résout le problème couplé en utilisant la méthode des éléments finis à partir de la formulation faible du problème.

Les inconnues nodales sont notées \({U}_{s}\) , \(P\) , \(\Phi ` et :math:\)Psi ` .

Soit \(({N}_{i}^{f})\) les fonctions de forme éléments finis pour le fluide et \(({N}_{i}^{s})\) les fonctions de forme éléments finis pour la structure. On utilisera les mêmes fonctions de forme pour \(P\) , \(\Phi ` et :math:\)Psi ` .

Formulation en u-p

Cette formation est la plus «naturelle» car elle utilise les variables physiques du problème. L’ensemble à résoudre est donc:

(954)\[\begin{split}\lbrace \begin{array}{c}{\rho}_{s}{\ddot{u}}_{s}-div({\sigma}_{s}({u}_{s}))=0\text{dans}{\Omega}_{s}\\ {\sigma}_{s}({u}_{s})\cdot {n}_{s}={f}_{s}\text{sur}\partial {\Omega}_{N}\\ {u}_{s}=\overline{u}\text{sur}\partial {\Omega}_{U}\\ {\sigma}_{s}\cdot {n}_{s}=-p{n}_{s}\text{sur}\Sigma \\ {\ddot{u}}_{s}\cdot {n}_{s}=-\frac{1}{{\rho}_{0}}\nabla p\cdot {n}_{s}\text{sur}\Sigma \\ \frac{1}{{c}_{0}^{2}}\ddot{p}-\Delta p=0\text{dans}{\Omega}_{f}\end{array}\end{split}\]

L’approximation par éléments finis du problème complet conduit au système suivant:

(955)\[\begin{split}\left[\begin{array}{cc}{M}_{s}& 0\\ -{\rho}_{0}{C}^{T}& {M}_{f}\end{array}\right]\left\lbrace \begin{array}{c}{\ddot{U}}_{s}\\ \ddot{P}\end{array}\right\rbrace +\left[\begin{array}{cc}{K}_{s}& C\\ 0& {K}_{f}\end{array}\right]\left\lbrace \begin{array}{c}{U}_{s}\\ P\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{F}_{s}\\ 0\end{array}\right\rbrace\end{split}\]

Les quantités \({M}_{s}\) , \({K}_{s}\) et \({F}_{s}\) sont relatives à la modélisation de la structure et nous renvoyons le lecteur vers les documentations adéquates. Pour le fluide, on a:

(956)\[ \begin{align}\begin{aligned}{({M}_{f})}_{ij}={\int}_{{\Omega}_{f}}\frac{1}{{c}_{0}^{2}}.{N}_{i}^{f}{N}_{j}.d\Omega\\\textrm{et}\\{({K}_{f})}_{ij}={\int}_{{\Omega}_{f}}\nabla {N}_{i}^{f}\cdot \nabla {N}_{j}.d\Omega\end{aligned}\end{align} \]

Et pour le couplagefluide-structure:

(957)\[{(C)}_{ij}={\int}_{\Sigma}{N}_{j}^{f}\left({N}_{i}^{s}\cdot {n}_{s}\right).d\Sigma\]

Les matrices \({M}_{s}\) et \({M}_{f}\) sont définies positives, les matrices \({K}_{s}\) et \({K}_{f}\) sont semi-définies positives (il faut ajouter les conditions limites en déplacement et pression). Le système résultant () n’est pas symétrique.

Pour obtenir une matrice de masse et de rigidité fluide rigoureusement nulle, par convention, il suffit de mettre \({\rho}_{0}={c}_{0}=0\) (soit les coefficients RHO et CELE_R dans DEFI_MATERIAU).

Pour la condition BGT, le système est modifié et vaut:

(958)\[\begin{split}\left[\begin{array}{cc}{M}_{s}& 0\\ -{\rho}_{0}{C}^{T}& {M}_{f}\end{array}\right]\left\lbrace \begin{array}{c}{\ddot{U}}_{s}\\ \ddot{P}\end{array}\right\rbrace +\left[\begin{array}{cc}0& 0\\ 0& \frac{{\rho}_{0}}{{Z}_{C}}Q\end{array}\right]\left\lbrace \begin{array}{c}{\dot{U}}_{s}\\ \dot{P}\end{array}\right\rbrace +\left[\begin{array}{cc}{K}_{s}& C\\ 0& {K}_{f}+\frac{{\rho}_{0}}{{Z}_{R}}Q\end{array}\right]\left\lbrace \begin{array}{c}{U}_{s}\\ P\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{F}_{s}\\ 0\end{array}\right\rbrace\end{split}\]

On ajoute donc un terme en amortissement avec la matrice \(Q\) telle que:

(959)\[{(Q)}_{ij}={\int}_{{\Sigma}_{R}}{N}_{i}^{f}{N}_{j}^{f}.d\Gamma\]

Formulation en u-p-phi

Le système précédent étant non-symétrique, la recherche de modes propres nécessite un solveur modal QEP, en général moins robuste que les solveurs modaux symétriques. C’est pour ça qu’on peut également employer une formulation de type U-P-PHI.

L’ensemble à résoudre est donc:

(960)\[\begin{split}\lbrace \begin{array}{c}{\rho}_{s}{\ddot{u}}_{s}-div({\sigma}_{s}({u}_{s}))=0\text{dans}{\Omega}_{s}\\ {\sigma}_{s}({u}_{s})\cdot {n}_{s}={f}_{s}\text{sur}\partial {\Omega}_{N}\\ {u}_{s}=\overline{u}\text{sur}\partial {\Omega}_{U}\\ {\sigma}_{s}\cdot {n}_{s}={\rho}_{0}\ddot{\varphi}{n}_{s}\text{sur}\Sigma \\ {\ddot{u}}_{s}\cdot {n}_{s}=\nabla \ddot{\varphi}\cdot {n}_{s}\text{sur}\Sigma \\ \frac{1}{{\rho}_{0}{c}_{0}^{2}}p+\frac{1}{{c}_{0}^{2}}\ddot{\varphi}=0\text{dans}{\Omega}_{f}\\ \frac{1}{{c}_{0}^{2}}\ddot{p}+{\rho}_{0}\Delta \ddot{\varphi}=0\text{dans}{\Omega}_{f}\end{array}\end{split}\]

L’approximation par éléments finis du problème complet conduit au système suivant:

(961)\[\begin{split}\left[\begin{array}{ccc}{M}_{s}& 0& -{\rho}_{0}C\\ 0& 0& {M}_{f}\\ -{\rho}_{0}{C}^{T}& {M}_{f}& -{\rho}_{0}{K}_{f}\end{array}\right]\left\lbrace \begin{array}{c}{\ddot{U}}_{s}\\ \ddot{P}\\ \ddot{\Phi}\end{array}\right\rbrace +\left[\begin{array}{ccc}{K}_{s}& 0& 0\\ 0& \frac{1}{{\rho}_{0}}{M}_{f}& 0\\ 0& 0& 0\end{array}\right]\left\lbrace \begin{array}{c}{U}_{s}\\ P\\ \Phi \end{array}\right\rbrace =\left\lbrace \begin{array}{c}{F}_{s}\\ 0\\ 0\end{array}\right\rbrace\end{split}\]

Les matrices sont les mêmes que dans le paragraphe précédent. On voit que le système est symétrique mais qu’il contient des zéros sur la diagonale. En particulier le terme de masse relatif à la variable \(p\) est problématique car il va rendre les calculs modaux plus délicats. En effet, si l’on considère que les valeurs propres du problème sont de la forme \(\sqrt{(\frac{K}{M})}\) , on voit que le terme nul va nécessiter:

• soit de calculer les valeurs propres dans une bande de fréquences excluant le zéro, ce qui nécessitera en particulier d’avoir une bonne idée a priori de cette bande de fréquences (en général, un pré-calcul sans couplage avec le fluide est utilisé pour ça);

• soit de modifier le système global en ajoutant une valeur sur la diagonale qui soit suffisamment grande pour rendre le système inversible, et pas trop pour ne pas perturber les résultats

De plus, la présence de la nouvelle variable :math:`Phi ` nécessite de définir des conditions limites sur cette variable.

Pour la condition BGT, le système est modifié et vaut:

(962)\[\begin{split}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& -\frac{{\rho}_{0}^{2}}{{Z}_{C}}Q\end{array}\right]\left\lbrace \begin{array}{c}{\stackrel{⃛}{U}}_{s}\\ \stackrel{⃛}{P}\\ \stackrel{⃛}{\Phi}\end{array}\right\rbrace +\left[\begin{array}{ccc}{M}_{s}& 0& -{\rho}_{0}C\\ 0& 0& {M}_{f}\\ -{\rho}_{0}{C}^{T}& {M}_{f}& -{\rho}_{0}{K}_{f}-\frac{{\rho}_{0}^{2}}{{Z}_{R}}Q\end{array}\right]\left\lbrace \begin{array}{c}{\ddot{U}}_{s}\\ \ddot{P}\\ \ddot{\Phi}\end{array}\right\rbrace +\left[\begin{array}{ccc}{K}_{s}& 0& 0\\ 0& \frac{1}{{\rho}_{0}}{M}_{f}& 0\\ 0& 0& 0\end{array}\right]\left\lbrace \begin{array}{c}{U}_{s}\\ P\\ \Phi \end{array}\right\rbrace =\left\lbrace \begin{array}{c}{F}_{s}\\ 0\\ 0\end{array}\right\rbrace\end{split}\]

On voit que l’ajout de cette condition rend le système d’ordre trois. Or, il n’existe pas de schéma robuste pour cet ordre, dans code_aster, on modifie donc le système en le rendant d’ordre deux par un terme d’amortissement et utilisant la relation définie dans le système d’équations , c’est-à-dire \(\frac{1}{{\rho}_{0}{c}_{0}^{2}}p+\frac{1}{{c}_{0}^{2}}\ddot{\varphi}=0\text{dans}{\Omega}_{f}\) . Le système devient alors:

(963)\[\begin{split}\left[\begin{array}{ccc}{M}_{s}& 0& -{\rho}_{0}C\\ 0& 0& {M}_{f}\\ -{\rho}_{0}{C}^{T}& {M}_{f}& -{\rho}_{0}{K}_{f}-\frac{{\rho}_{0}^{2}}{{Z}_{R}}Q\end{array}\right]\left\lbrace \begin{array}{c}{\ddot{U}}_{s}\\ \ddot{P}\\ \ddot{\Phi}\end{array}\right\rbrace +\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& \frac{{\rho}_{0}}{{Z}_{C}}Q& 0\end{array}\right]\left\lbrace \begin{array}{c}{\dot{U}}_{s}\\ \dot{P}\\ \dot{\Phi}\end{array}\right\rbrace +\left[\begin{array}{ccc}{K}_{s}& 0& 0\\ 0& \frac{1}{{\rho}_{0}}{M}_{f}& 0\\ 0& 0& 0\end{array}\right]\left\lbrace \begin{array}{c}{U}_{s}\\ P\\ \Phi \end{array}\right\rbrace =\left\lbrace \begin{array}{c}{F}_{s}\\ 0\\ 0\end{array}\right\rbrace\end{split}\]

On remarque que la nouvelle matrice d’ordre 1 devient non-symétrique.

Formulation en u-psi

La troisième formulation proposée dans code_aster est celle utilisant la variable \(\psi\) (potentiel de vitesse) pour le fluide. L’ensemble à résoudre est:

(964)\[\begin{split}\lbrace \begin{array}{c}{\rho}_{s}{\ddot{u}}_{s}-div({\sigma}_{s}({u}_{s}))=0\text{dans}{\Omega}_{s}\\ {\sigma}_{s}({u}_{s})\cdot {n}_{s}={f}_{s}\text{sur}\partial {\Omega}_{N}\\ {u}_{s}=\overline{u}\text{sur}\partial {\Omega}_{U}\\ {\sigma}_{s}\cdot {n}_{s}={\rho}_{0}\dot{\psi}{n}_{s}\text{sur}\Sigma \\ {\dot{u}}_{s}\cdot {n}_{s}=\nabla \psi \cdot {n}_{s}\text{sur}\Sigma \\ \frac{1}{{c}_{0}^{2}}\ddot{\psi}-\Delta \psi =0\text{dans}{\Omega}_{f}\end{array}\end{split}\]

L’approximation par éléments finis du problème complet conduit au système suivant:

(965)\[\begin{split}\left[\begin{array}{cc}{M}_{s}& 0\\ 0& -{\rho}_{0}{M}_{f}\end{array}\right]\left\lbrace \begin{array}{c}{\ddot{U}}_{s}\\ \ddot{\Psi}\end{array}\right\rbrace +\left[\begin{array}{cc}0& -{\rho}_{0}C\\ -{\rho}_{0}{C}^{T}& 0\end{array}\right]\left\lbrace \begin{array}{c}{\dot{U}}_{s}\\ \dot{\Psi}\end{array}\right\rbrace +\left[\begin{array}{cc}{K}_{s}& 0\\ 0& -{\rho}_{0}{K}_{f}\end{array}\right]\left\lbrace \begin{array}{c}{U}_{s}\\ \Psi \end{array}\right\rbrace =\left\lbrace \begin{array}{c}{F}_{s}\\ 0\end{array}\right\rbrace\end{split}\]

Ce système est symétrique mais il contient un terme en vitesse, ce qui nécessite en particulier un solveur modal quadratique.

Pour la condition BGT, le système est modifié et vaut:

(966)\[\begin{split}\left[\begin{array}{cc}{M}_{s}& 0\\ 0& -{\rho}_{0}{M}_{f}\end{array}\right]\left\lbrace \begin{array}{c}{\ddot{U}}_{s}\\ \ddot{\Psi}\end{array}\right\rbrace +\left[\begin{array}{cc}0& -{\rho}_{0}C\\ -{\rho}_{0}{C}^{T}& \frac{-{\rho}_{0}^{2}}{{Z}_{C}}Q\end{array}\right]\left\lbrace \begin{array}{c}{\dot{U}}_{s}\\ \dot{\Psi}\end{array}\right\rbrace +\left[\begin{array}{cc}{K}_{s}& 0\\ 0& -{\rho}_{0}({K}_{f}+\frac{{\rho}_{0}}{{Z}_{C}}Q)\end{array}\right]\left\lbrace \begin{array}{c}{U}_{s}\\ \Psi \end{array}\right\rbrace =\left\lbrace \begin{array}{c}{F}_{s}\\ 0\end{array}\right\rbrace\end{split}\]

Discrétisation par éléments finis

Si l’on impose des conditions limites de type vitesse imposée ou impédance de paroi imposée au fluide, on est conduit à résoudre le système matriciel suivant:

(967)\[\begin{split}\left[\begin{array}{ccc}K& 0& 0\\ 0& \frac{{M}_{f}}{{\rho}_{0}{c}^{2}}& 0\\ 0& 0& 0\end{array}\right]-{\omega}^{2}\left[\begin{array}{ccc}M& 0& {\rho}_{0}{M}_{\Sigma}\\ 0& 0& \frac{{M}_{\mathrm{fl}}}{{c}^{2}}\\ {\rho}_{0}{M}_{\Sigma}^{T}& \frac{{M}_{\mathrm{fl}}^{T}}{{c}^{2}}& {\rho}_{0}H\end{array}\right]+i{\omega}^{3}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& \frac{{\rho}_{0}^{2}}{Z}Q\end{array}\right]=\left[\begin{array}{}0\\ i\omega \\ V\end{array}\right]\end{split}\]

\(\mathrm{Q}\) étant la matrice obtenue à partir de la forme bilinéaire \({\int}_{{\Sigma}_{z}}{\Phi}^{2}\mathit{dS}\) et \(\mathrm{V}\) , le vecteur obtenu à partir de \({\int}_{{\Sigma}_{V}}{\rho}_{0}{v}_{0}\Phi \mathit{dS}\) .

\(\mathrm{Q}\) étant la matrice obtenue à partir de la forme bilinéaire \({\int}_{{\Sigma}_{z}}{\Phi}^{2}\mathit{dS}\) et \(\mathrm{V}\) , le vecteur obtenu à partir de \({\int}_{{\Sigma}_{V}}{\rho}_{0}{v}_{0}\Phi \mathit{dS}\) .

Utilisation

Les éléments décrits précédemment appartiennent, pour la partie fluide, à la modélisation “3D_FLUIDE”, “2D_FLUIDE” ou “AXIS_FLUIDE” du phénomène MECANIQUE et, pour l’interface fluide-structure, à la modélisation “FLUI_STRU” ou “2D_FLUI_STRU”. Pour définir les éléments absorbants fluides, qui prennent en compte la condition BGT, il faut affecter les modélisations “3D_FLUI_ABSO”, “2D_FLUI_ABSO” ou “AXIS_FLUI_ABSO”.

Pour choisir la formulation (\(u-p\) , \(u-p-\phi\) ou \(u-\psi\) ), on utilisera le mot-clef FORMULATION dans AFFE_MODELE.

Ils conduisent à des éléments volumiques ou surfaciques, pour la partie fluide, et à des éléments surfaciques ou linéiques pour l’interface fluide-structure ou pour les éléments absorbants.

Flambement d’Euler, raideur géométrique

On considère le problème de flambement d’Euler d’un solide élastique en interaction avec un domaine fluide, modélisé avec l’hypothèse vibro-acoustique. On modélise des champs fluctuants. Le flambement d’Euler est donc traité par la recherche de l’annulation des fréquences propres de vibrations du système couplé, en incluant la présence de la raideur géométrique associée à l’état statique précontraint dans le solide, contrôlée par un paramètre de chargement critique à déterminer.

On admet dans un premier temps que la raideur géométrique associée à l’état statique initial dans le fluide est négligeable. On propose donc de définir une matrice de raideur géométrique nulle dans le domaine fluide et sur les interfaces fluide/structure. Cependant cette méthode d’analyse néglige le fait que la pression fluide est un chargement suiveur pour le solide, ce qui introduit donc un terme non linéaire, qui a souvent des effets notables sur les charges critiques de flambement. Il est donc préférable de traiter le problème par une analyse dynamique non linéaire associée à une analyse de stabilité, en pilotant le chargement mécanique : consulter le document U2.06.11.

Bibliographie

1. 20. BALANNEC, S. COURTIER-ARNOUX, E. LUZZATO, P. THOMAS : Confrontation d’outils numériques du département AMV sur un test de couplage fluide‑structure. Rapport interne E.D.F. D.E.R. HP-61/91.012.

2. 16. GERMAIN : Introduction à la mécanique des milieux continus. Masson.

3. R.J. GIBERT : Vibrations des structures. Interactions avec les fluides. Sources d’excitation aléatoires. Collection de la Direction des Etudes et Recherche d’E.D.F. n°69.

4. 18. OHAYON, R. VALID : True symmetric formulations of free vibrations of fluid‑structure interaction. Applications and extensions. International conference on numerical method for coupled systems.

5. 16. THOMAS : Formulation de couplage fluide-structure pour l’analyse modale dans SIVA. Application aux éléments tridimensionnels. Rapport interne E.D.F. D.E.R. HP-35.82/259.

6. 1. Bayliss, M. Gunzburger, E. Turkel, Boundary Conditions for the Numerical Solution of Elliptic Equations in Exterior Regions, SIAM J. Appl. Math., 42(2), 430–451

7. 20. Rayleigh, The theory of sound, 1894