r4.02.04 Couplage Fluide - Structure avec surface Libre#
Résumé :
On présente ici le couplage fluide / structure dans le cas où le fluide présente une surface libre. Des éléments de surface libre ont été implantés dans Code_Aster pour calculer les modes de ballotement d’un fluide couplé à une structure élastique pour un problème tridimensionnel.
Table des matières
Introduction#
Afin d’étudier le comportement de structures remplies de fluide, on peut être conduit à prendre en compte les phénomènes de ballottement c’est-à-dire ajouter au système couplé fluide‑structure, l’effet de la pesanteur au niveau de la surface libre du liquide. Les structures concernées sont, par exemple, les cuves de centrales nucléaires de la filière rapide, les piscines de stockage du combustible [bib4].
On a donc complété les développements déjà réalisés en couplage fluide-structure [bib3] par l’introduction de nouveaux éléments surfaciques qui prennent en compte, dans leur formulation, l’effet de la pesanteur.
Formulation théorique du problème#
Le problème d’interaction structure-fluide pesant revient à résoudre simultanément trois problèmes :
la structure est soumise à un champ de pression \(P\) imposé par le fluide sur la paroi \(\Sigma\) ;
le fluide est soumis à un champ de déplacement \({X}_{s}\) imposé par la structure sur \(\Sigma\) ;
la pesanteur agit sur la surface libre où \(p=\rho gz\) .
On considérera dans un premier temps que le fluide est non pesant avant d’introduire la gravité au paragraphe [§3.2].
Rappels sur le couplage fluide-structure#
Afin de bien rendre compte de l’interaction fluide-structure, nous allons analyser séparément les équations régissant le comportement du fluide et celles qui régissent celui de la structure, sans considérer dans ce chapitre les conditions aux limites concernant la surface libre.
Description du fluide#
On considère que le système étudié est soumis à de petites perturbations autour de son état d’équilibre où le fluide et la structure sont au repos : ainsi, \(P={p}_{0}+p\) et \({\mathrm{X}}_{\mathrm{s}}={\mathrm{x}}_{\mathrm{s}}({x}_{0}=0)\) . Ce qui permet d’écrire [bib2] :
\(\rho =-{\rho}_{f}div({\mathrm{x}}_{\mathrm{f}})\) d’où \(p=-{\rho}_{f}{c}^{2}div({\mathrm{x}}_{\mathrm{f}})\) .
Avec :
\(p\) la pression fluctuante du fluide,
\(\rho ` la perturbation de masse volumique du fluide, :math:`{\rho}_{f}\) la masse volumique du fluide au repos,
\({x}_{f}(r,t)\) le champ de déplacement d’une particule de fluide.
Le fluide est :
parfait (c’est-à-dire non visqueux)
barotrope :
\(p=\rho {c}^{2}\) éq 3.1.1-1
et irrotationnel : il existe un potentiel des déplacements \(\varphi\) , tel que \(p={\rho}_{f}\frac{{\partial}^{2}\varphi }{\partial {t}^{2}}\)
Le comportement du volume de fluide eulérien est donc décrit par les équations suivantes :
loi de comportement :
\({T}_{ij}=-p{\delta}_{ij}\)
équation de conservation de la quantité de mouvement dans le fluide en l’absence de source:
\(div(\overline{\overline{T}})={\rho}_{f}\frac{{\partial}^{2}{x}_{f}}{\partial {t}^{2}}\) éq 3.1.1-2
équation de la conservation de la masse :
\(\frac{\partial \rho }{\partial t}+{\rho}_{f}div(\frac{\partial {x}_{s}}{\partial t})=0\) éq 3.1.1-3
En combinant les équations de conservation de la quantité de mouvement [éq 3.1.1-2] et de la masse [éq 3.1.1-3] écrites en régime harmonique à la pulsation \(w\) , on obtient, grâce à [éq 3.1.1-1], l’équation de Helmholtz :
\(\Delta p+\frac{{\omega}^{2}}{{c}^{2}}p=0\)
Description de la structure#
On considère que la structure est élastique linéaire et qu’on reste dans le domaine des petites perturbations. Compte tenu de ces hypothèses, on écrit :
la loi de comportement en élasticité linéaire :
\({\sigma}_{ij}={C}_{ijkl}{\varepsilon}_{kl}\)
l’équation de conservation de la quantité de mouvement dans la structure en l’absence de forces volumiques autres que les forces d’inertie :
\(div(\overline{\overline{\sigma}})={\rho}_{s}\frac{{\partial}^{2}{x}_{s}}{\partial {t}^{2}}\)
l’équation de compatibilité sur le tenseur de déformation :
\({\varepsilon}_{kl}=\frac{1}{2}({u}_{k,l}+{u}_{l,k})\)
Description de l’interface fluide-structure#
A l’interface (\(\Sigma\) ) entre le fluide et la structure, comme le fluide n’est pas visqueux, il y a continuité des contraintes normales et des vitesses normales à la paroi, et nullité de la contrainte tangentielle (absence de frottement visqueux). Ces conditions aux limites s’écrivent :
\(\lbrace \begin{array}{c}{\sigma}_{ij}{n}_{i}={T}_{ij}{n}_{i}=-p{\delta}_{ij}{n}_{i}\\ \frac{\partial {x}_{f}}{\partial t}.n=\frac{\partial {x}_{s}}{\partial t}.n\end{array}\)
Formulation du problème couplé#
Finalement, l’équation du problème couplé fluide-structure s’écrit, en prenant \(p\) comme variable décrivant le champ de pression dans le fluide et \({x}_{s}\) le champ des déplacements dans la structure :
\(\lbrace \begin{array}{cc}{C}_{ijkl}{x}_{{s}_{k,\mathrm{lj}}}+{\omega}^{2}{\rho}_{s}{x}_{{s}_{i}}=0& \text{dans}{V}_{s}\\ \Delta p+\frac{{\omega}^{2}}{{c}^{2}}p=0& \text{dans}{V}_{f}\\ {\sigma}_{ij}{n}_{i}={C}_{ijkl}{x}_{{s}_{k,l}}{n}_{i}=-p{\delta}_{ij}{n}_{i}& \text{sur}\Sigma \\ \frac{\partial p}{\partial n}={\rho}_{f}{\omega}^{2}{x}_{{f}_{i}}{n}_{i}& \text{sur}\Sigma \end{array}\)
Les champs de déplacements \({x}_{s}\) pour la structure et de pression \(p\) pour le fluide cherchés minimisent la fonctionnelle :
\(L({x}_{s},p,z)=\frac{1}{2}\underset{{V}_{s}}{\int}\left[{\sigma}_{ij}({x}_{s}){\varepsilon}_{ij}({x}_{s})-{\rho}_{s}{\omega}^{2}{x}_{s}^{2}\right]-\underset{\Sigma}{\int}p{x}_{s}nd\Sigma +\frac{1}{2{\rho}_{f}}\underset{{V}_{f}}{\int}\left[\frac{1}{{\omega}^{2}}(\text{grad}p{)}^{2}-\frac{{p}^{2}}{{c}^{2}}\right]\text{dV}\)
Action de la pesanteur sur la surface libre#
Formulation du problème#
On rappelle ici les équations dynamiques linéarisées décrivant les petits mouvements d’un fluide parfait. On choisit une description eulérienne du fluide :
\(\text{grad}P={\rho}_{f}(\mathrm{g}-\frac{{\partial}^{2}{\mathrm{x}}_{\mathrm{s}}}{\partial {t}^{2}})\) dans \({V}_{f}\)
A l’équilibre la particule de fluide était en \({M}_{0}\) et donc : \(\text{grad}{P}_{0}={\rho}_{f}\mathrm{g}\) dans \({V}_{f}\) .
On envisage des mouvements de faible amplitude autour de l’état d’équilibre (c’est l’hypothèse des petites perturbations) : alors \(M={M}_{0}+{x}_{f}({M}_{0},t)\)
Soient \(p\) la fluctuation de pression eulérienne et \({p}_{L}\) la fluctuation de pression lagrangienne, alors :
\(\begin{array}{}p(M,t)=P({M}_{0},t)-{P}_{0}({M}_{0})\\ {p}_{L}=P(M,t)-{P}_{0}({M}_{0})\end{array}\)
Compte tenu de l’hypothèse des petits déplacements :
\(\begin{array}{}{p}_{L}-p=\text{gradP}({M}_{o},t){x}_{f}({M}_{o},t)\\ \phantom{{p}_{L}-p}=-{\rho}_{f}g{x}_{f}({M}_{o},t)\end{array}\) éq 3.2.1-1
Si on considère le cas d’un fluide pesant présentant une surface libre en contact avec un milieu à pression constante \({P}_{\text{atm}}\) , on peut écrire, en négligeant les effets de tension superficielle :
\(P(M,t)={P}_{\text{atm}}\) sur la surface libre \(\text{SL}\) c’est-à-dire : \({p}_{L}=0\) . Soit, avec [éq 3.2.1-1], \(p={\rho}_{f}g({x}_{f}.z)\)
Compte tenu de l’hypothèse des petits mouvements, l’inclinaison instantanée du plan tangent est un infiniment petit du premier ordre. \({x}_{f}.z\) se confond donc au deuxième ordre près avec l’élévation verticale \(h\) .
Figure 3.2.1-a : approximation sur la surface libre
Ainsi, si on rajoute aux conditions aux limites la condition de pesanteur sur la surface libre, cela revient à considérer en \(z=H\) la condition linéarisée :
\(p={\rho}_{f}gz\) éq 3.2.1-2
Les équations du problème global deviennent :
\(\lbrace \begin{array}{cc}{C}_{ijkl}{x}_{{s}_{k,\text{lj}}}+{\omega}^{2}{\rho}_{s}{x}_{{s}_{i}}=0& \text{dans}{V}_{s}\\ \Delta p+\frac{{\omega}^{2}}{{c}^{2}}p=0& \text{dans}{V}_{f}\\ {\sigma}_{ij}{n}_{i}={C}_{ijkl}{x}_{{s}_{k,l}}{n}_{i}=-p{\delta}_{ij}{n}_{i}& \text{sur}\Sigma \\ \frac{\partial p}{\partial n}={\rho}_{f}{\omega}^{2}{x}_{{f}_{i}}{n}_{i}& \text{sur}\Sigma \text{et}\text{sur}\text{SL}\\ p={\rho}_{f}gz& \text{sur}\text{SL}\end{array}\)
Pour exprimer la fonctionnelle, on utilise la loi de comportement sur la surface libre. En considérant un champ de déplacement admissible \(\mathit{dz}\) on obtient [bib2] :
\(\underset{\text{SL}}{\int}\rho gz\delta z\text{ds}=\underset{\text{SL}}{\int}p\delta z\text{ds}\)
Soit, finalement, la fonctionnelle du système total fluide structure soumis à la pesanteur :
\(\begin{array}{}L({x}_{s},p,z)=\frac{1}{2}\underset{{V}_{s}}{\int}\left[{\sigma}_{ij}({x}_{s}){\varepsilon}_{ij}({x}_{s})-{\rho}_{s}{\omega}^{2}{x}_{s}^{2}\right]-\underset{\Sigma}{\int}p{x}_{s}.nd\Sigma +\frac{1}{2}{\rho}_{f}\underset{{V}_{f}}{\int}\left[\frac{1}{{\omega}^{2}}(\text{grad}p{)}^{2}-\frac{{p}^{2}}{{c}^{2}}\right]\text{dV}\\ \phantom{L({x}_{s},p,z)}+\frac{1}{2}\underset{\text{SL}}{\int}\rho g{z}^{2}\text{dS}-\underset{\text{SL}}{\int}pz\text{dS}\end{array}\)
Cette prise en compte de la pesanteur implique deux termes supplémentaires dans la fonctionnelle décrivant le fluide :
un terme d’énergie potentielle liée à la surface libre : \(\frac{1}{2}\underset{\text{SL}}{\int}\rho g{z}^{2}\text{ds}\)
un terme dû au travail de la pression hydrodynamique dans le déplacement de la surface libre : \(\underset{\text{SL}}{\int}pz\text{ds}\)
Cependant il faut noter que ce n’est pas l’unique effet de la pesanteur puisqu’en tout point de la paroi \(\Sigma\) s’exerce une pression permanente \(-\rho gz\) (où \(z\) est l’altitude du point \(M\) considéré : on suppose que \(z=0\) au niveau de la surface libre à l’équilibre). Le point \(M\) est animé d’un mouvement \({X}_{s}\) infinitésimal, l’élément de surface \(d\Sigma\) varie donc et l’effort dû à la pression permanente aussi. Cet effort est responsable d’un terme de rigidité supplémentaire s’ajoutant à la rigidité de structure dans le système. Il pourrait causer un flambage de la structure en annulant la rigidité structurale. Cet effet est négligeable sur les caractéristiques vibratoires ([bib2], [bib1]), on ne le prend donc pas en compte.
Discrétisation par éléments finis#
Pour obtenir la forme discrétisée de la fonctionnelle, on remplace chaque intégrale par une somme d’intégrales sur chaque élément \(i\) du système discrétisé, puis on utilise une approximation par éléments finis des fonctions inconnues de déplacement et de pression sur chaque élément \(i\) [bib18].
Les inconnues sont \({\mathrm{X}}_{\mathrm{s}}(u,v,w),p,z\) , on a alors en posant \(\text{Ni}\) les fonctions de formes (ou fonctions d’interpolation nodales sur l’élément \(i\) ) :
\(\lbrace \begin{array}{c}\sigma =D\varepsilon \\ \varepsilon =Bu\end{array}\) \(\lbrace \begin{array}{c}{x}_{s}={N}_{i}u\\ p(x,y,z)={N}_{i}p\\ \nabla p=\overline{{N}_{i}}p\end{array}\) \(\lbrace \begin{array}{c}{x}_{s}.n={N}_{\Sigma i}u\\ z={N}_{\mathrm{Si}}z\end{array}\)
où \(\delta ,p\) , sont les inconnues aux nœuds structures et aux nœuds fluides, et \(z\) les inconnues à la surface libre.
D’où l’expression discrétisée de la fonctionnelle associée au problème :
\(L={u}^{t}(K-{\omega}^{2}M)u+{p}^{t}(\frac{H}{{\rho}_{f}{\omega}^{2}}-\frac{Q}{{\rho}_{f}{c}^{2}})p+{z}^{t}{\rho}_{f}g{K}_{z}z-2{p}^{t}{M}_{z}z-2{p}^{t}Cu\)
avec
\(K=\sum_{i}\underset{\text{Vi}}{\int}{N}_{i}^{t}{B}_{i}^{t}{D}_{i}{B}_{i}{N}_{i}{\text{dV}}_{i}\) |
matrice raideur de la structure |
\(M=\sum_{i}\underset{\text{Si}}{\int}{N}_{i}^{t}{\rho}_{f}{N}_{i}{\text{dV}}_{i}\) |
matrice masse de la structure |
et
\(Q=\sum_{i}\underset{{\text{Vi}}_{\text{fl}}}{\int}{N}_{i}^{t}{N}_{i}{\text{dV}}_{i}\) |
\({M}_{z}=\sum_{i}\underset{\mathrm{Si}}{\int}{N}_{\Sigma i}^{t}{N}_{\Sigma i}{\text{dS}}_{i}\) |
\({K}_{z}=\sum_{i}\underset{\text{Si}}{\int}{N}_{\text{Si}}^{t}{N}_{\text{Si}}{\text{dS}}_{i}\) |
\({M}_{z}=\sum_{i}\underset{\text{Si}}{\int}{N}_{i}^{t}{N}_{\text{Si}}{\text{dS}}_{i}\) |
\(H=\sum_{i}\underset{{\text{Vi}}_{\text{fl}}}{\int}\overline{{N}_{i}^{t}}\overline{{N}_{i}}{\text{dV}}_{i}\) |
où \(c\) est la célérité du son dans le fluide, \({\rho}_{f}\) la masse volumique du fluide et où \({K}_{f}\) correspond à l’énergie potentielle du fluide, \({K}_{z}\) à l’énergie potentielle de la surface libre, \(H\) à l’énergie cinétique du fluide, \(M\) au couplage fluide-solide et \({M}_{z}\) au couplage \(p-z\) sur la surface libre.
L’approximation par éléments finis du problème complet conduit alors au système matriciel suivant :
\(\left[\begin{array}{ccc}K& -C& 0\\ 0& H& 0\\ 0& -{M}_{z}& {K}_{z}\end{array}\right]\left\lbrace \begin{array}{c}u\\ p\\ z\end{array}\right\rbrace -{\omega}^{2}\left[\begin{array}{ccc}M& 0& 0\\ {\rho}_{f}C& \frac{Q}{{c}^{2}}& {\rho}_{f}{M}_{z}\\ 0& 0& 0\end{array}\right]\left\lbrace \begin{array}{c}u\\ p\\ z\end{array}\right\rbrace =0\)
La première équation correspond au mouvement de la structure soumise aux forces de pression, la deuxième à celle du mouvement du fluide couplé à la structure et à la surface libre, la troisième est l’équation de surface libre.
Cependant le problème écrit de la sorte possède des matrices masse et rigidité non symétriques ce qui empêche l’utilisation des algorithmes de résolution classiques de Code_Aster .
Introduction d’une variable supplémentaire#
Pour rendre le problème symétrique et pouvoir utiliser les méthodes de résolution classiques, on introduit une variable supplémentaire : le potentiel des déplacements dans le fluide \(j\) [bib2].
\({\mathrm{X}}_{\mathrm{f}}=\text{grad}\varphi\) i.e. \(\rho {\omega}^{2}\varphi =p\)
Cette inconnue supplémentaire est liée aux inconnues du problème, ce qui conduit à une matrice de rigidité singulière.
On reformule le problème couplé structure-fluide pesant :
\(\lbrace \begin{array}{cc}{C}_{ijkl}{x}_{{s}_{k,\text{lj}}}+{\omega}^{2}{\rho}_{s}{x}_{{s}_{i}}=0& \text{dans}{V}_{s}\\ \Delta \varphi +\frac{{\omega}^{2}}{{\rho}_{f}{c}^{2}}p=0& \text{dans}{V}_{f}\\ p={\rho}_{f}{\omega}^{2}\varphi & \text{dans}{V}_{f}\\ {\sigma}_{ij}{n}_{i}={C}_{ijkl}{x}_{{s}_{k,l}}{n}_{i}=-{\rho}_{f}{\omega}^{2}\varphi {\delta}_{ij}{n}_{i}& \text{sur}S\\ \frac{\partial \varphi }{\partial n}={x}_{{f}_{i}}{n}_{i}& \text{sur}S\\ p={\rho}_{f}gz& \text{sur}\text{SL}\end{array}\)
Ce qui conduit à la fonctionnelle du système couplé :
\(\begin{array}{}L({x}_{s},p,j,z)=\frac{1}{2}\underset{{V}_{s}}{\int}\left[{\sigma}_{ij}({x}_{s}){\varepsilon}_{ij}({x}_{s})-{\rho}_{s}{\omega}^{2}{x}_{s}^{2}\right]+\frac{1}{2}\underset{{V}_{f}}{\int}\frac{{p}^{2}}{{c}^{2}}\text{dV}+\frac{1}{2}\underset{\text{SL}}{\int}\rho g{z}^{2}\text{ds}\\ \phantom{L({x}_{s},p,j,z)}-{\omega}^{2}\left[\underset{\Sigma}{\int}{\rho}_{f}\varphi {x}_{s}nd\Sigma +\underset{\text{SL}}{\int}{\rho}_{f}\varphi z\text{ds}+\underset{{V}_{f}}{\int}\left[\frac{{\rho}_{f}}{2}(\text{grad}\varphi {)}^{2}+\frac{p\varphi }{{c}^{2}}\right]\text{dV}\right]\end{array}\)
Soit en discrétisant :
\(\begin{array}{}L={\delta}^{t}(K-{\omega}^{2}M)\delta +\frac{1}{{\rho}_{f}{c}^{2}}{p}^{t}Qp+{\rho}_{f}g{z}^{t}{K}_{z}z\\ \phantom{\text{L=}}-2{\omega}^{2}\left[\frac{{\rho}_{f}}{2}{\varphi}^{t}H\varphi +{\rho}_{f}{\varphi}^{t}C\delta +{\rho}_{f}{\varphi}^{t}{M}_{z}z+\frac{1}{{c}^{2}}{p}^{t}Q\varphi \right]\end{array}\)
Ce qui s’écrit, sous forme matricielle :
\(\left[\begin{array}{cccc}K& 0& 0& 0\\ 0& \frac{Q}{{\rho}_{f}{c}^{2}}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& {\rho}_{f}g{K}_{z}\end{array}\right]\left[\begin{array}{c}\delta \\ p\\ \varphi \\ z\end{array}\right]-{\omega}^{2}\left[\begin{array}{cccc}M& 0& {\rho}_{f}C& 0\\ 0& 0& \frac{Q}{{c}^{2}}& 0\\ {\rho}_{f}{C}^{t}& \frac{{Q}^{t}}{{c}^{2}}& {\rho}_{f}H& {\rho}_{f}{M}_{z}\\ 0& 0& {\rho}_{f}{M}_{z}^{t}& 0\end{array}\right]\left[\begin{array}{c}\delta \\ p\\ \varphi \\ z\end{array}\right]=0\)
Implantation dans Code_Aster#
La bibliothèque des éléments finis de Code_Aster a été enrichie de cinq éléments surfaciques iso-paramétriques ayant comme degrés de liberté la déflexion de la surface libre et le potentiel des déplacements du fluide à la surface libre. Ils sont compatibles avec les éléments 3D qui traitent le problème de couplage fluide/structure [bib3]
On nomme :
MEFP_FACE3 et MEFP_FACE6 respectivement les triangles à 3 ou à 6 nœuds,
MEFP_FACE4, MEFP_FACE8 et MEFP_FACE9 respectivement les quadrangles à 4, 8 ou à 9 nœuds.
Ces éléments appartiennent à la modélisation 2D_FLUI_PESA du phénomène MECANIQUE.
Bibliographie#
J.TANI, J. TERAKI : Free vibration analysis of FBR vessels partially filled with liquid. SMIRT 1989
R.J. GIBERT : Vibration des sructures - interaction avec les fluides - sources d’excitation aléatoires. CEA/EDF/INRIA 1988
WAECKEL : Analyse modale en vibration acoustique dans ASTER. Note interne HP‑61/91160 EDF/DER
LEPOUTERE, F. WAECKEL : Effet de la pesanteur sur la surface libre d’un fluide couplé à une structure, Note interne HP - 61/93.139
Description des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
3 |
G. ROUSSEAU, Fe WAECKEL (EDF/EP/AMV) |
Texte initial |