r4.02.05 Éléments de frontière absorbante#
Résumé
Ce document décrit l’implantation dans code_aster des éléments de frontière absorbante. Ces éléments de type paraxiaux, dont on décrit ici la théorie, sont affectés à des frontières de domaines élastiques ou fluides pour traiter des problèmes 2D ou 3D d’interaction sol-structure ou sol-fluide-structure. Ils permettent de satisfaire la condition de Sommerfeld vérifiant l’hypothèse d’anéchoïcité : l’élimination des ondes planes élastiques ou acoustiques diffractées et non physiques venant de l’infini.
Théorie des éléments paraxiaux#
On présente dans cette partie le principe de l’approximation paraxiale dans le cas de l’élastodynamique linéaire. Deux approches théoriques permettent de cerner l’esprit et la mise en pratique des éléments paraxiaux élastiques : on doit la première à Cohen et Jennings [bib2] et la seconde à Modaressi [bib3]. L’application de la théorie des éléments paraxiaux au cas fluide sera faite dans la partie suivante.
Dans toute la suite, comme présenté sur la r4.02.05-domaine_isfs, on suppose que la frontière du maillage du sol est située dans un domaine au comportement élastique. L’approche de Modaressi implantée dans code_aster permet à la fois de construire des frontières absorbantes et d’introduire la champ sismique incident.
Impédance spectrale de la frontière#
Pour obtenir l’équation paraxiale, il nous faut d’abord déterminer la forme du champ de déplacement diffracté au voisinage de la frontière. Pour cela, on part des équations de l’élastodynamique 3D:
Avec le déplacement \(\disp\) qui se décompose en une partie normale et tangentielle :
et
La constante \(\soundSpeed\), homogène à une vitesse, est introduite pour rendre certaines quantités adimensionnelles. Les équations et leurs solutions sont bien entendu indépendantes de cette constante.
On appelle \(\posiTang\) et \(\dispTang\) les directions et les composantes du déplacement dans le plan tangent et \(\posiNorm\) et \(\dispNorm\) selon \(\basisVector{3}\), la direction normale à la frontière.
On procède à deux transformées de Fourier, l’une par rapport au temps, l’autre par rapport aux variables d’espace dans le plan à la frontière. On se limite au cas d’une frontière plane et sans coin:
Les équations s’écrivent alors. Selon la composante tangentielle:
et selon la composante normale
où \(\widehat{\disp}\) et \(\widehat{\dispNorm}\) désignent les transformées de Fourier et \(\vector{\xi}_{\tau}\) le vecteur d’onde associé à \(\posiTang\).
Il s’agit d’un système différentiel en \(\posiNorm\) que l’on sait résoudre en le diagonalisant. On en déduit:
Avec \(\vector{\xi}_P=\sqrt{\frac{{\omega}^{2}}{\soundSpeedComp^{2}}-{\normEucl{\vector{\xi}_{\tau}}}^{2}}\) et \(\vector{\xi}_S=\sqrt{\frac{{\omega}^{2}}{\soundSpeedCisa^{2}}-{\normEucl{\vector{\xi}_{\tau}}}^{2}}\).
Pour déterminer les constantes \(A\), \({A}_{S}\) et \({A}_{P}\), on suppose connu \(\dispTang \left(\vector{\xi}, 0\right)\) sur la frontière du domaine éléments finis. On les exprime en fonction de \(\widehat{\dispTang}\left(\vector{\xi}_{\tau}, 0\right)={\widehat{\dispTang}}_{,0}\) et \(\widehat{\dispNorm}\left(\vector{\xi}_{\tau},0\right)={\widehat{\dispNorm}}_{,0}\).
On va maintenant évaluer le vecteur contrainte \(\sigmVector\) sur une facette de normale \(\basisVector{3}\) en \(\posiNorm=0\), ce qui nous donnera l’impédance de la frontière. On fait subir à \(\sigmVector(\posiTang, \posiNorm)\) la même transformée de Fourier en espace que pour les équations de l’élastodynamique, si bien que:
\(\lameLambda\) et \(\lameMu\) sont les coefficients élastiques de Lamé.
On souhaite s’affranchir en \(\posiNorm=0\) des termes contenant des dérivées en \(\posiNorm\) . Le système obtenu précédemment nous le permet en fonction de \({\widehat{\dispTang}}_{,0}\) et \({\widehat{\dispNorm}}_{,0}\).
On obtient ainsi l’impédance spectrale de la frontière:
où \({a}^{0}\), \({b}^{0}\) et \({c}^{0}\) sont des fonctions de \(\normEucl{\vector{\xi}_{\tau}}\) et de \(\omega\) qui dépendent linéairement de \(\widehat{\dispTang}_{,0}\) et \({\widehat{\dispNorm}_{,0}}\) .
On peut alors écrire :
où \(A\) désigne l’opérateur impédance spectrale globale. On revient à l’espace physique par deux transformées de Fourier inverses.
Approximation paraxiale de l’impédance#
L’impédance spectrale calculée précédemment n’est locale ni en espace ni en temps puisqu’elle fait intervenir \(\left( \vector{\xi}_{\tau},\omega\right)\), la transformée de Fourier de \({\disp}_{0} \left(\posiTang, t\right)\) pour tout \(\posiTang\) et tout \(t\) .
L’idée est alors de développer \({\vector{\xi}}_{P}\) et \({\vector{\xi}}_{S}\) selon les puissances de \(\frac{\normEucl{\vector{\xi}_{\tau}}}{\omega}\). Cette approximation sera bonne soit à haute fréquence, soit pour \(\normEucl{\vector{\xi}_{\tau}}\) petit.
Examinons la dépendance en \(\posiNorm\), par exemple de \(\widehat{\dispNorm}\). On aura, pour \(\dispNorm\left(\posiTang, \posiNorm,t\right)\), des termes de la forme \(\exp\left[i\left(\vector{\xi}_{\tau}\posiTang+\omega t-{\vector{\xi}}_{P}\posiNorm\right)\right]\).
Avec le développement de \(\vector{\xi}_P\):
On montre que pour \(\normEucl{\vector{\xi}_{\tau}}\) petit, on aura des ondes se propageant selon des directions voisines de la normale \(\basisVector{3}\) à la frontière, car l’exponentielle s’écrit:
Dès lors, avec un développement asymptotique de \({\vector{\xi}}_{P}\) et \({\vector{\xi}}_{S}\), en multipliant par une puissance convenable de \(\omega\) pour supprimer cette quantité au dénominateur, on obtient :
où \({A}_{0}\) et \({A}_{1}\) sont des fonctions polynomiales en \(\vector{\xi}_{\tau}\) et \(\omega\).
Soit, après les deux transformées de Fourier inverses:
On obtient ainsi la forme finale de l’impédance locale transitoire approchée en fonction du dernier terme en \(\frac{\normEucl{\vector{\xi}_{\tau}}}{\omega}\) retenu. On peut trouver le calcul détaillé des \({A}_{i}\) dans [bib4].
Par exemple, pour l’ordre 0 :
Ceci correspond à des amortisseurs visqueux distribués le long de la frontière du domaine éléments finis, avec le coefficient correctif \(C^{a}\) renseigné derrière le mot clé COEF_AMOR affecté par [DEFI_MATERIAU] au matériau des éléments de frontière absorbante.
Par analogie, pour l’ordre 0, on peut ajouter un terme en déplacement correspondant à des rigidités distribuées le long de la frontière du domaine éléments finis :
Cette formulation supplémentaire comprend les coefficients de Lamé \(\lameLambda\) et \(\lameMu\) , ainsi qu’une dimension caractéristique \(L\) .
A l’ordre 1:
On voit apparaître la dérivée par rapport au temps du vecteur contrainte \(\sigmVector\). Dans le traitement numérique, il faudra avoir recours à une intégration de ce terme sur les éléments de la frontière.
Pour conclure, on retiendra que l’approximation paraxiale conduit à une impédance locale transitoire ne faisant intervenir essentiellement que des dérivées en temps et dans le plan tangent à la frontière, ainsi qu’éventuellement un terme supplémentaire sans dérivée en temps.
De façon symbolique, on écrit à l’ordre 0:
et à l’ordre 1:
Prise en compte du champ sismique incident#
On rappelle que le comportement du sol est supposé élastique au moins au voisinage de la frontière. On introduit donc le champ diffracté \({\disp}_r\) tel que:
A l’infini, le champ total \(\disp\) doit être égal au champ incident \({\disp}_{i}\) (une des conséquences de la condition de radiation de Sommerfeld), c’est-à-dire:
A la frontière du maillage éléments finis, on écrit la condition d’absorption pour le champ diffracté à l’ordre 0:
On en déduit le vecteur contrainte total sur la frontière du maillage éléments finis:
On obtient ainsi la formulation variationnelle du problème au voisinage de la frontière pour l’ordre 0:
Pour tout champ \(\dispVirt\) cinématiquement admissible.
La condition d’absorption pour le champ diffracté à l’ordre 1 s’écrit:
Pour l’ordre 1, on conserve la formulation classique:
où \(\sigmVector \left(\disp\right)\) suit la loi d’évolution suivante:
La sollicitation due au champ incident apparaît explicitement dans le cas de l’ordre 0, mais elle est contenue dans la loi d’évolution de \(\sigmVector \left(\disp\right)\) pour l’ordre 1.
Éléments fluides anéchoïques en transitoire#
Cette partie présente l’essentiel des contraintes générales d’implémentation d’éléments fluides anéchoïques de frontière absorbants avec l’approximation paraxiale d’ordre 0 dans code_aster. Pour des raisons de simplicité liées à la manipulation de grandeurs scalaires telles que la pression ou le potentiel de déplacement, en opposition aux grandeurs vectorielles comme le déplacement, on s’intéresse d’abord aux éléments fluides.
Formulation standard#
On reprend ici le raisonnement de Modaressi en l’adaptant à un domaine fluide acoustique. Dans un premier temps, on s’intéresse à la seule donnée de la grandeur pression dans ce fluide. On reviendra ensuite sur cette modélisation pour s’adapter aux contraintes de code_aster, en soulignant les ajustements à faire.
Soit donc la configuration suivante, en reprenant les conventions de la partie précédente au voisinage de la frontière:
La définition d’un repère local au niveau de l’élément permet de nous ramener systématiquement dans une telle situation.
Formulation éléments finis#
La pression \(p\) vérifie l’équation d’Helmholtz dans tout le domaine \(\Omega\) modélisé aux éléments finis, ce qui donne, pour tout champ de pression virtuel \(q\) :
\(-\underset{\Omega}{\int}\nabla p.\nabla q-\frac{1}{{c}^{2}}\frac{{\partial}^{2}}{\partial {t}^{2}}\underset{\Omega}{\int}\text{pq}+\underset{\Sigma}{\int}\frac{\partial p}{\partial n}q=0\)
\(\Sigma\) représente la frontière du domaine \(O\) .
La grandeur à estimer sur \(\Sigma\) grâce à l’approximation paraxiale est ici \(\frac{\partial p}{\partial n}\) .
Approximation paraxiale#
Dans la configuration proposée, le terme \(\frac{\partial p}{\partial n}\) correspond à \(\frac{\partial p}{\partial {x}_{3}}\) .
Considérons dès lors une onde plane harmonique se propageant dans le fluide:
\(p=A\exp\left[i({k}_{1}{x}_{1}+{k}_{2}{x}_{2}+{k}_{3}{x}_{3}-\omega t)\right]\)
En remplaçant dans l’équation d’Helmholtz, on obtient:
\({k}_{3}=\frac{\omega}{c}\sqrt{1-\frac{{c}^{2}}{{\omega}^{2}}({k}_{1}^{2}+{k}_{2}^{2})}\)
On obtient alors le développement suivant, pour des hautes fréquences (\(\omega\) grand) ou au voisinage de la frontière (\({k}_{1}\) et \({k}_{2}\) petits):
\({k}_{3}=\frac{\omega}{c}(1-\frac{{c}^{2}}{2{\omega}^{2}}({k}_{1}^{2}+{k}_{2}^{2}))\)
Soit, en multipliant par \(\omega\) pour faire disparaître cette quantité au dénominateur et après une transformée de Fourier inverse en espace et en temps:
\(\frac{{\partial}^{2}p}{\partial {x}_{3}\partial t}=-\frac{1}{c}\frac{{\partial}^{2}p}{\partial {t}^{2}}+\frac{1}{2}c(\frac{{\partial}^{2}p}{\partial {x}_{1}^{2}}+\frac{{\partial}^{2}p}{\partial {x}_{2}^{2}})\)
Comme l’avait présenté Modaressi, cette équation fait intervenir la dérivée par rapport au temps du terme surfacique. Dans le cadre de cette partie, on ne s’intéresse qu’au terme d’ordre 0, soit, après une intégration en temps, ce qui fait disparaître la dérivée gênante:
\(\frac{\partial p}{\partial {x}_{3}}=-\frac{1}{c}\frac{\partial p}{\partial t}\) ou plus généralement: \(\frac{\partial p}{\partial n}=-\frac{1}{c}\frac{\partial p}{\partial t}\)
C’est cette relation d’impédance que nous allons discrétiser sur la frontière du domaine éléments finis.
Remarque :
Compte tenu de la disparition du terme d’ordre 1 dans le développement de la racine carrée, l’ordre minimal d’approximation pour les paraxiaux fluides est en fait 1 et non 0. Nous conserverons l’appellation d’éléments d’ordre 0 pour la cohérence avec le solide. Toutefois, on parle d’éléments fluides d’ordre 2 au moment de considérer des éléments d’ordre strictement positif.
Impédance des éléments vibro-acoustiques dans code_aster#
code_aster dispose d’éléments vibro-acoustiques. On rappelle dans ce paragraphe les choix de formulation faits au moment de leur implémentation. On s’inspire pour présenter l’existant de la documentation de référence r4.02.02.
Limites de la formulation en p#
Dans le cadre de l’interaction fluide-structure en harmonique, la formulation en pression uniquement du fluide acoustique conduit à des matrices non symétriques. En effet, le système global s’exprime, sous forme variationnelle, de la façon suivante:
\(\underset{{\Omega}_{s}}{\int}{C}_{ijkl}.{u}_{k,l}{v}_{i,j}-{\omega}^{2}\underset{{\Omega}_{s}}{\int}{\rho}_{s}{u}_{i}{v}_{i}-\underset{\Sigma}{\int}p{v}_{i}.{n}_{i}=0\) pour la structure
\(\frac{1}{{\rho}_{f}{\omega}^{2}}\underset{{\Omega}_{f}}{\int}\nabla p.\nabla q-{k}^{2}\underset{{\Omega}_{f}}{\int}\mathrm{pq}-\underset{\Sigma}{\int}{u}_{i}.{n}_{i}q=0\) pour le fluide
avec \(k=\frac{\omega}{c}\) , nombre d’onde pour le fluide, \(v\) et \(q\) deux champs virtuels dans la structure et dans le fluide respectivement.
Après discrétisation par éléments finis, on obtient le système matriciel suivant:
\(\left[\begin{array}{cc}K& -C\\ 0& H\end{array}\right]\left[\begin{array}{c}u\\ p\end{array}\right]-{\omega}^{2}\left[\begin{array}{cc}M& 0\\ {\rho}_{f}{C}^{T}& \frac{Q}{{c}^{2}}\end{array}\right]\left[\begin{array}{c}u\\ p\end{array}\right]=0\)
où: |
\(K\) et \(M\) sont les matrices de rigidité et de masse de la structure |
\(H\) et \(Q\) sont les matrices fluides obtenues respectivement à partir des formes bilinéaires: \(\underset{{\Omega}_{f}}{\int}\nabla p.\nabla q\) et \(\underset{{O}_{f}}{\int}\mathrm{pq}\) |
\(C\) est la matrice de couplage obtenue à partir de la forme bilinéaire: \(\underset{\Sigma}{\int}p{u}_{i}.{n}_{i}\)
Le caractère non symétrique de ce système ne permet pas d’utiliser les algorithme de résolution classiques de code_aster. Ceci motive l’introduction d’une variable supplémentaire dans la description du fluide.
Formulation symétrique en p et phi#
La nouvelle grandeur introduite est le potentiel des déplacements \(\phi\) , tel que \(u=\nabla \varphi\) . D’après , on obtient la nouvelle forme variationnelle du système couplé fluide-structure:
\(\underset{{\Omega}_{s}}{\int}{C}_{ijkl}.{u}_{k,l}{v}_{i,j}-{\omega}^{2}\underset{{\Omega}_{s}}{\int}{\rho}_{s}{u}_{i}{v}_{i}-{\rho}_{f}{\omega}^{2}\underset{\Sigma}{\int}\phi p{v}_{i}.{n}_{i}=0\) pour la structure
\(\frac{1}{{\rho}_{f}{c}^{2}}\underset{{\Omega}_{f}}{\int}\text{pq}-{\rho}_{f}{\omega}^{2}\left[\frac{1}{{\rho}_{f}{c}^{2}}\underset{{\Omega}_{f}}{\int}\left(\varphi q+p\chi \right)-\underset{{\Omega}_{f}}{\int}\nabla \varphi .\nabla \chi +\underset{S}{\int}\chi {\mathrm{u}}_{i}{\mathrm{n}}_{i}\right]=0\) pour le fluide
Avec: \(p={\rho}_{f}{\omega}^{2}f\) dans le fluide et \(\chi\) un champ de potentiel de déplacement virtuel
Ceci nous conduit au système matriciel symétrique:
\(\left[\begin{array}{ccc}K& 0& 0\\ 0& \frac{{M}_{f}}{{\rho}_{f}{c}^{2}}& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}u\\ p\\ \phi \end{array}\right]-{\omega}^{2}\left[\begin{array}{ccc}M& 0& {\rho}_{f}{M}_{\Sigma}\\ 0& 0& \frac{{M}_{\text{fl}}}{{c}^{2}}\\ {\rho}_{f}{M}_{\Sigma}^{T}& \frac{{M}_{\text{fl}}^{T}}{{c}^{2}}& {\rho}_{f}H\end{array}\right]\left[\begin{array}{c}u\\ p\\ \phi \end{array}\right]=0\)
où: \(K\) et \(\mathrm{M}\) sont les matrices de rigidité et de masse de la structure
\({\mathrm{M}}_{\Sigma}\) est la matrice de couplage obtenue à partir de la forme bilinéaire \(\underset{\Sigma}{\int}\varphi {\mathrm{u}}_{i}{\mathrm{n}}_{i}\)
\({\mathrm{M}}_{\mathrm{f}},{\mathrm{M}}_{\text{fl}}\) et \(\mathrm{H}\) sont les matrices fluides obtenues à partir des formes bilinéaires: \(\underset{{\Omega}_{f}}{\int}\text{pq}\) , \(\underset{{\Omega}_{f}}{\int}pq\) (ou \(\underset{{\Omega}_{f}}{\int}\varphi q\) ) et \(\underset{{\Omega}_{f}}{\int}\nabla \varphi .\nabla \chi\)
Formulation en psi#
Elle existe une troisième formulation définie en \(\psi\) , c’est-à-dire un potentiel de vitesse. Pour plus de détails on renvoie à la documentation . Cette formulation est également disponible dans code_aster. Son point négatif est de ne pas avoir un accès direct au post-traitement en pression du fluide.
Imposition d’une impédance d’ordre 0#
D’une manière générale, une relation d’impédance à la frontière du fluide s’exprime ainsi:
\(p=Z\mathrm{v}.\mathrm{n}\)
où:
\(Z=\rho c{q}_{\alpha}\) est l’impédance imposée (\(\rho\) est la densité du fluide, \(c\) est la célérité des ondes fluides et \({q}_{\alpha}\) est le coefficient de réflexion);
\(\mathrm{v}.\mathrm{n}\) est la vitesse normale sortante des particules fluides.
On en déduit, d’après la loi de comportement du fluide, qui relie la pression au déplacement des particules fluides pour un fluide acoustique \(\nabla p-{\rho}_{f}\frac{{\partial}^{2}u}{\partial {t}^{2}}=0\) :
\(\frac{{\rho}_{f}}{Z}\frac{\partial p}{\partial t}=\frac{\partial p}{\partial n}\)
La discrétisation d’un telle équation conduit à un terme non symétrique dans une formulation en \(p\) et \(\varphi\) . On préfère formuler la condition par rapport au potentiel de déplacement, soit:
\(\nabla \varphi +\frac{{\rho}_{f}}{Z}\frac{\partial \varphi }{\partial t}=0\)
On obtient alors comme expression pour le terme de bord associé à la relation d’impédance:
\({\rho}_{f}\frac{{\partial}^{2}}{\partial {t}^{2}}\underset{\Sigma}{\int}\varphi \frac{\partial \psi }{\partial n}=\frac{{\partial}^{3}}{\partial {t}^{3}}\underset{\Sigma}{\int}\frac{{\rho}_{f}^{2}}{Z}\varphi \psi\)
On constate alors l’apparition (quelque peu artificielle) d’un terme en dérivée troisième par rapport au temps. En harmonique, cela ne pose pas de problème. On traite un terme en \({\omega}^{3}\) sans difficulté. Ce terme de troisième ordre sera calculé à l’aide de l’option IMPE_MECA dans CALC_MATR_ELEM (ou ASSEMBLAGE) et utilisait comme entrée de MATR_IMPE_PHI dans DYNA_VIBRA (TYPE_CALCUL = ‘HARM’). Pour le calcul transitoire, plutôt que d’introduire une approximation d’une dérivée troisième dans le schéma de Newmark implémenté dans les opérateurs d’intégration directe en dynamique dans code_aster DYNA_VIBRA (ou DYNA_LINE) et DYNA_NON_LINE, on peut transformer le terme de troisième ordre dans un terme de premier ordre et traiter cela comme un terme d’amortissement classique (pour plus de détails faire référence à ). Ce terme sera calculé à l’aide de l’option AMOR_MECA dans CALC_MATR_ELEM (ou ASSEMBLAGE). Cette option est disponible pour les trois formulations fluides (pour les détails sur la formulation \(u-\psi\) faire référence à ).
Imposition d’une impédance d’ordre 1#
La condition de rayonnement est exacte pour une onde plane. Cependant, pour approcher asymptotiquement le comportement à l’infini d’une onde sphérique, pour un rayon \(R\) donné, on utilise une approximation d’ordre 1. Cette condition est rappelée dans , ainsi que la description complète des termes imposés par cette condition. En fonction de la formulation choisie, on va impacte la masse ou la raideur du modèle:
Si on choisit la formulation \(u-p-\varphi\) , on rajoute un terme de masse calculé à l’aide de l’option MASS_MECA dans CALC_MATR_ELEM (ou ASSEMBLAGE);
Si on choisit les formulations \(u-p\) ou \(u-\psi\) , on rajoute un terme de raideur calculé à l’aide de l’option RIGI_MECA dans CALC_MATR_ELEM (ou ASSEMBLAGE).
Le rayon \(R\) sera donné via l’option LONG_CARA dans DEFI_MATERIAU(FLUIDE…). Ce paramètre est nul par défaut (dans ce cas l’impédance d’ordre 1 n’est pas prise en compte).
Le taux de réflexion#
Si on veut représenter un taux de réflexion \(\alpha\) , par exemple dans le cas de la présence d’alluvions en fond de retenue, il faut alors corriger la valeur de la célérité \(c\) par une nouvelle valeur \(c'\) telle que \(c'=c\frac{1+\alpha }{1-\alpha }\) . Le paramètre \(\alpha\) sera défini via l’option COEF_AMOR dans DEFI_MATERIAU (FLUIDE…). Ce paramètre est nul par défaut (dans ce cas la frontière est parfaitement absorbante). Ce paramètre influence à la fois l’impédance d’ordre 0 et d’ordre 1.
Formulation détaillée#
On propose ici la formulation précise pour un fluide acoustique modélisé sur un domaine \(\Omega ` avec une condition anéchoïque sur une partie :math:`{\Sigma}_{a}\) de la frontière \(\Sigma ` du domaine. En dehors de cela, on décompose la frontière en une surface libre et une partie en contact avec un solide rigide. L’introduction de sollicitations extérieures ou la présence d’une structure élastique se modélise aisément par les méthodes courantes. Les éléments de volume et de surface sont formulés en :math:`p\) et \(\varphi\) .
Les équations dans le fluide sont:
\({\rho}_{f}\Delta \varphi +\frac{1}{{c}^{2}}p=0\) dans le volume :math:`Omega ` éq 3.2.4-1
\(p={\rho}_{f}\frac{{\partial}^{2}\varphi }{\partial {t}^{2}}\) dans le volume :math:`Omega ` éq 3.2.4-2
\(p=0\) sur la surface libre éq 3.2.4-3
\(\frac{\partial \varphi }{\partial n}=0\) sur la paroi rigide éq 3.2.4-4
\(\frac{\partial p}{\partial n}=-\frac{1}{c}\frac{\partial p}{\partial t}\) sur la partie de la frontière avec condition anéchoïque éq 3.2.4-5
On multiplie l’équation [éq 3.2.4-1] par un champ de potentiel virtuel \(y\) et on intègre dans :math:`Omega ` :
\(\underset{{\Omega}_{f}}{\int}\left[\frac{1}{{c}^{2}}p\psi +{\rho}_{f}\frac{{\partial}^{2}}{\partial {t}^{2}}\left(\nabla \varphi .\nabla \psi \right)\right]+\underset{\Sigma}{\int}\psi {\rho}_{f}\frac{{\partial}^{2}}{\partial {t}^{2}}\left(\frac{\partial \varphi }{\partial n}\right)=0\) d’après la formule de Green
Soit, avec les conditions aux limites sur \(\Sigma ` et l’équation [:ref:`éq 3.2.4-2 <éq 3.2.4-2>\)]:
\(\underset{{\Omega}_{f}}{\int}\left[\frac{1}{{c}^{2}}p\psi +{\rho}_{f}\frac{{\partial}^{2}}{\partial {t}^{2}}\left(\nabla \varphi .\nabla \psi \right)\right]+\underset{\Sigma}{\int}\psi {\rho}_{f}\frac{\partial p}{\partial n}=0\)
On peut dès lors appliquer la condition d’impédance formulée en pression:
\(\underset{{\Sigma}_{a}}{\int}\psi {\rho}_{f}\frac{\partial p}{\partial n}=-\frac{1}{c}\underset{{\Sigma}_{a}}{\int}\psi {\rho}_{f}\frac{\partial p}{\partial t}\)
De plus, pour parvenir à une formulation symétrique des termes de volume, on multiplie l’équation [éq3.2.4-2] par un champ de pression virtuel \(q\) et on intègre dans :math:`Omega ` :
\(\underset{{\Omega}_{f}}{\int}\frac{\text{pq}}{{\rho}_{f}{c}^{2}}-\frac{{\partial}^{2}}{\partial {t}^{2}}\underset{{\Omega}_{f}}{\int}\frac{\varphi q}{{c}^{2}}=0\)
En sommant les deux équations variationnelles, on obtient:
\(\frac{1}{{\rho}_{f}{c}^{2}}\underset{{\Omega}_{f}}{\int}\text{pq}+{\rho}_{f}\frac{{\partial}^{2}}{\partial {t}^{2}}\left[\frac{1}{{\rho}_{f}{c}^{2}}\underset{{\Omega}_{f}}{\int}\left(\varphi q+p\psi \right)-\underset{{\Omega}_{f}}{\int}\nabla \varphi .\nabla \psi \right]-\frac{1}{c}\underset{{\Sigma}_{a}}{\int}\psi {\rho}_{f}\frac{\partial p}{\partial t}=0\)
Matriciellement:
\(\left[\begin{array}{cc}{\mathrm{M}}_{\mathrm{f}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}p\\ \varphi \end{array}\right]-\frac{1}{c}\left[\begin{array}{cc}0& 0\\ \mathrm{A}& 0\end{array}\right]\left[\begin{array}{c}\dot{p}\\ \dot{\varphi}\end{array}\right]+\left[\begin{array}{cc}0& \frac{{\mathrm{M}}_{\text{fl}}}{{c}^{2}}\\ \frac{{\mathrm{M}}_{\text{fl}}^{\mathrm{T}}}{{c}^{2}}& {\rho}_{f}\mathrm{H}\end{array}\right]\left[\begin{array}{c}\ddot{p}\\ \ddot{\varphi}\end{array}\right]=0\)
où les sous-matrices \({\mathrm{M}}_{\mathrm{f}},{\mathrm{M}}_{\text{fl}}\) et \(\mathrm{H}\) discrétisent les mêmes formes bilinéaires précédentes.
La sous-matrice \(\mathrm{A}\) discrétise le terme \(\underset{{\Sigma}_{a}}{\int}\psi {\rho}_{f}\frac{\partial p}{\partial t}\) . On remarque que la matrice d’amortissement obtenue n’est pas symétrique.
Utilisation dans code_aster#
La prise en compte d’éléments fluides anéchoïques et le calcul de leur impédance nécessite une modélisation spécifique sur les frontières absorbantes:
en 2D avec la modélisation “2D_FLUI_ABSO’ou “AXIS_FLUI_ABSO” sur les arêtes absorbantes à \(n\) nœuds;
en 3D avec la modélisation “3D_FLUI_ABSO” sur les faces absorbantes à \(n\) nœuds.
En analyse harmonique pour la formulation \(u-p-\varphi\) avec l’opérateur DYNA_VIBRA, on calcule au préalable une impédance mécanique par l’option IMPE_MECA de l’opérateur CALC_MATR_ELEM et on la renseigne dans DYNA_LINE_HARM (mot clé MATR_IMPE_PHI). Dans tous les autres cas, l’impédance d’ordre 0 sera prise en compte via l’option AMOR_MECA de l’opérateur CALC_MATR_ELEM. L’impédance d’ordre 1 sera prise en compte via des termes de rigidité ou masse (en fonction de la formulation) comme précédemment expliqué. Les caractéristiques matériaux des frontières absorbantes fluides seront définies dans DEFI_MATERIAU (FLUIDE…).
Éléments absorbants élastiques dans code_aster#
Cette partie présente l’essentiel des contraintes générales d’implémentation d’éléments élastiques de frontière absorbants avec l’approximation paraxiale d’ordre 0 dans code_aster. On rappelle la relation d’impédance paraxiale d’ordre 0 telle qu’elle a été établie par Modaressi pour un domaine élastique linéaire:
\(\mathrm{t}\left(\mathrm{u}\right)=\rho \left({c}_{p}\frac{\partial {\mathrm{u}}_{\phantom{.}\perp \phantom{.}}}{\partial t}+{c}_{s}\frac{\partial {\mathrm{u}}_{\text{//}}}{\partial t}\right)\)
\({\mathrm{u}}_{\mathrm{T}}\) devient \({\mathrm{u}}_{3}\) et \({\mathrm{u}}_{\text{//}}\) devient \({\mathrm{u}}^{\text{/}}\)
Adaptation du chargement sismique aux éléments paraxiaux#
On a présenté dans la première partie le principe de prise en compte du champ incident grâce aux éléments paraxiaux. Il convient ici de présenter les méthodes de modélisation du chargement sismique dans code_aster pour pouvoir adapter les données aux exigences des éléments paraxiaux.
L’équation fondamentale de la dynamique associée à un modèle quelconque 2D ou 3D discrétisé en éléments finis de milieu continu ou de structure et en l’absence de chargement extérieur s’écrit dans le repère absolu:
\(\mathrm{M}{\ddot{\mathrm{X}}}_{\mathrm{a}}+\mathrm{C}{\dot{\mathrm{X}}}_{\mathrm{a}}+{\text{KX}}_{\mathrm{a}}=0\)
On décompose le mouvement des structures en un mouvement d’entraînement \({\mathrm{X}}_{\mathrm{e}}\) et un mouvement relatif \({\mathrm{X}}_{\mathrm{r}}\) .
Figure 4.1-a: Décomposition du mouvement des structures
Ainsi, \({\mathrm{X}}_{\mathrm{a}}={\mathrm{X}}_{\mathrm{r}}+{\mathrm{X}}_{\mathrm{e}}\)
\({\mathrm{X}}_{\mathrm{a}}\) est le vecteur des déplacements dans le repère absolu,
\({\mathrm{X}}_{\mathrm{r}}\) est le vecteur des déplacements relatifs, c’est-à-dire le vecteur des déplacements de la structure par rapport à la déformée qu’elle aurait sous l’action statique des déplacements imposés au niveau des supports \({\mathrm{X}}_{\mathrm{S}}\) . \({\mathrm{X}}_{\mathrm{r}}\) est donc nul aux points d’ancrage,
\({\mathrm{X}}_{\mathrm{e}}\) est le vecteur des déplacements d’entraînement de la structure produit statiquement par le déplacement imposé des supports \({\mathrm{X}}_{\mathrm{S}}\) : \({\mathrm{X}}_{\mathrm{e}}=\Psi {\mathrm{X}}_{\mathrm{s}}\) ,
\(\Psi ` est la matrice des modes statiques. Les modes statiques représentent la réponse de la structure à un déplacement unitaire imposé à chaque degré de liberté de liaison (les autres étant bloqués), en l’absence de forces extérieures. Ainsi, :math:\)mathrm{K}Psi =0` , c’est-à-dire, \({\text{KX}}_{e}=0\) .
Dans le cas du mono-appui (tous les appuis subissent le même mouvement imposé), :math:`Psi ` est un mode de corps rigide.
Hypothèse dans code_aster:
On suppose que l’amortissement dissipé par la structure est de type visqueux c’est-à-dire que la force d’amortissement est proportionnelle à la vitesse relative de la structure. Ainsi, \(\mathrm{C}{\dot{\mathrm{X}}}_{e}=0\) .
L’équation fondamentale de la dynamique dans le repère relatif s’écrit alors:
\(\mathrm{M}{\ddot{\mathrm{X}}}_{\mathrm{r}}+\mathrm{C}{\dot{\mathrm{X}}}_{\mathrm{r}}+{\text{KX}}_{\mathrm{r}}=-\mathrm{M}\Psi \ddot{{\mathrm{X}}_{\mathrm{S}}}\)
L’opérateur CALC_CHAR_SEISME [U4.63.01] calcule le terme \(-\mathrm{M}\Psi ` , ou plus exactement , :math:\)-mathrm{M}Psi mathrm{d}` où \(\mathrm{d}\) est un vecteur unitaire tel que \({\mathrm{X}}_{\mathrm{s}}=\mathrm{d}.f\left(t\right)\) avec \(f\) une fonction scalaire du temps.
On distingue deux types de chargements sismiques introduits dans code_aster grâce à l’opérateur CALC_CHAR_SEISME:
Le chargement de type MONO_APPUI, pour lequel :math:`Psi ` est la matrice identité (les modes statiques sont des modes de corps rigide),
Le chargement de type MULTI_APPUI, pour lequel :math:`Psi ` est quelconque.
D’après la méthode de prise en compte du champ incident avec les éléments paraxiaux présentée dans la première partie, il nous faut connaître sur la frontière le déplacement et les contraintes dus au champ incident. Pour le chargement de type MULTI_APPUI, seul le déplacement est directement accessible à tout instant. Il semble donc difficile de permettre l’utilisation d’un tel mode de chargement avec des éléments paraxiaux dans le sol. D’ailleurs, si un tel chargement modélise les déplacements imposés des appuis, il ne requiert pas une modélisation du sol puisque toute l’influence est prise en compte par ces déplacements.
Le cas MONO_APPUI peut être perçu différemment. Il représente une accélération d’ensemble appliquée au modèle. Dès lors, la propagation d’onde dans le sol peut avoir un rôle à jouer dans le comportement de la structure, puisque les mouvements de l’interface sol-structure ne sont pas imposés. De plus, les éléments paraxiaux sont utilisables avec ce type de chargement car il ne crée pas de contraintes à la frontière du maillage (un mode de corps rigide ne crée pas de déformations). Dès lors, on dispose de toutes les données nécessaires au calcul de l’impédance absorbante sur la frontière.
Remarque1 :
Dans le cas d’une sollicitation sismique MONO_APPUI, le calcul dynamique se fait dans le repère relatif. Si on revient sur le terme à discrétiser sur les éléments paraxiaux (voir première partie), on remarque que \({u}_{i}\) correspond exactement au déplacement d’entraînement \({\mathrm{X}}_{\mathrm{e}}\) présenté plus haut. Ainsi, \(\mathrm{u}-{\mathrm{u}}_{i}\) correspond au déplacement relatif calculé pendant le calcul. Dès lors, la relation à prendre en compte sur les éléments paraxiaux dans une telle configuration est simplement:
\(t(u)={A}_{0}(\frac{\partial u}{\partial t})\)
Remarque 2:
Dans le cas d’un calcul d’interaction sol-fluide-structure avec fluide infini, la pression à prendre en compte pour le calcul de l’impédance anéchoïque dans le fluide est bien la pression absolue, si on n’a pas de champ incident dans le fluide (ce qui est souvent le cas). La correction que l’on pouvait se dispenser de faire pour le sol doit alors être faite pour les éléments paraxiaux fluides.
Implémentation des éléments en transitoire et en harmonique#
Implémentation en transitoire#
Le mode d’implémentation des éléments paraxiaux élastiques en transitoire est très voisin de celui présenté pour les éléments fluides. La différence vient essentiellement de la nécessité de décomposer le déplacement en une composante selon la normale à l’élément, correspondant à une onde \(P\) , et une composante dans le plan de l’élément, correspondant à une onde \(S\) . On est alors à même de discrétiser la relation d’impédance introduite dans la première partie:
\(\mathrm{t}\left(\mathrm{u}\right)=\rho {C}_{p}\frac{\partial {\mathrm{u}}_{3}}{\partial t}+\rho {C}_{s}\frac{\partial {\mathrm{u}}^{\text{/}}}{\partial t}\)
On ne revient pas sur le schéma d’intégration temporelle que l’on a déjà décrit dans la partie précédente, sachant qu’on considère la relation d’impédance par un opérateur d’amortissement ajouté au premier membre.
Pour la prise en compte du terme supplémentaire:
\({\mathrm{t}}_{1}\left(\mathrm{u}\right)=\frac{\lambda +2\mu }{L}{u}_{3}{\mathrm{e}}_{3}+\frac{\mu}{L}{\mathrm{u}}^{'}\)
on utilise cette fois-ci un opérateur de rigidité ajoutée au premier membre.
Implémentation en harmonique#
Les éléments acoustiques fluides de code_aster proposent déjà la possibilité de prendre en compte une impédance imposée à la frontière du maillage en harmonique. Cela correspond au traitement d’un terme en \({\omega}^{3}\) dans les équations, comme cité plus haut. Il est tentant d’introduire la possibilité d’imposer une impédance absorbante pour un problème élastique en harmonique.
Pour un calcul de réponse harmonique d’une structure infinie, la prise en compte de l’impédance absorbante comme une correction du second membre n’est évidemment pas applicable. Cependant, la relation d’impédance à l’ordre 0 exprime les termes surfaciques en fonction de la vitesse des nœuds de l’élément. On peut donc construire une pseudo-matrice d’amortissement visqueux traduisant la présence du domaine infini.
De la même façon, on construit (cf. 4.1.1) une pseudo-matrice de rigidité complétant le rôle de la matrice d’amortissement définie précédemment.
La décomposition de la relation d’impédance selon les composantes normale ou tangentielle du déplacement sur l’élément nous contraint à construire la matrice d’impédance dans un repère local sur l’élément. On définit ce repère local dans la routine élémentaire ainsi que la matrice de passage qui permet le retour à la base globale.
Mode de chargement sismique par onde plane#
En complément des méthodes de prise en compte du chargement sismique déjà disponible et en raison de l’inadéquation du mode MULTI_APPUI avec les éléments paraxiaux, il semble intéressant d’introduire un principe de chargement par onde plane. Cela correspond aux chargements classiquement rencontrés lors des calculs d’interaction sol-structure par les équations intégrales.
Caractérisation d’une onde plane en transitoire#
En harmonique, une onde plane élastique est caractérisée par sa direction, sa pulsation et son type (onde \(P\) pour les ondes de compression, ondes \(\mathrm{SV}\) ou \(\mathrm{SH}\) pour les ondes de cisaillement). En transitoire, la donnée de la pulsation, correspondant à une onde stationnaire en temps, doit être remplacée par la donnée d’un profil de déplacement dont on va prendre en compte la propagation au cours du temps dans la direction de l’onde.
Les directions des ondes :math:`P`, :math:`\mathit{SV}`et :math:`\mathit{SH}`sont déterminées à partir du vecteur :math:`V`renseigné par le paramètre DIRECTION. A savoir :
* :math:`P`est colinéaire à :math:`V`et normé à 1,
* :math:`\mathit{SH}`est l'intersection du plan horizontal et du plan normal à :math:`V`, et normé à 1,
* :math:`\mathit{SV}`est le produit vectoriel de :math:`\mathit{SH}`et de :math:`P`. Il existe un cas d'indétermination avec cette règle quand le plan horizontal et le plan normal sont confondus. Dans ce cas, si :math:`V=Z`purement vertical, on impose :math:`\mathit{SH}=Y`, et :math:`\mathit{SV}=X`.
Plus précisément, on va considérer une onde plane sous la forme:
\(u(x,t)=f(k.x-{C}_{p}t)k\) pour une onde \(P\) (avec \(k\) unitaire)
\(\mathrm{u}\left(\mathrm{x},t\right)=f\left(k.\mathrm{x}-{C}_{S}t\right)\wedge k\) pour une onde \(S\) (avec toujours \(k\) unitaire)
\(f\) représente alors le profil de l’onde donné selon la direction \(k\) .
\(H\) est la distance de l’origine au front d’onde principal.
Données utilisateur pour le chargement par onde plane#
Conformément à la théorie exposée en première partie, il nous faut calculer la contrainte à la frontière du maillage due à l’onde incidente et le terme d’impédance correspondant au déplacement incident, soit:
\(\mathrm{t}\left({\mathrm{u}}_{i}\right)\) et \({A}_{0}\left(\frac{\partial {\mathrm{u}}_{i}}{\partial t}\right)+{A}_{1}\left({\mathrm{u}}_{i}\right)\)
Pour exprimer les contraintes, il nous faut disposer des déformations dues à l’onde incidente, la loi de comportement du matériau nous permettant de passer des unes aux autres.
Sur les éléments de frontière, on peut exprimer le tenseur des déformations linéarisées en chaque nœud par la formule classique:
\(\epsilon \left(x,t\right)=\frac{1}{2}\left[\nabla \mathrm{u}\left(x,t\right)+{}^{t}\text{}\nabla \mathrm{u}\left(x,t\right)\right]\)
Finalement, pour estimer les contraintes dues au champ incident, il nous faut donc déterminer les dérivées \(\frac{\partial {\left({\mathrm{u}}_{i}\right)}_{j}}{\partial {x}_{k}}\) pour \(j\) et \(k\) parcourant les trois directions de l’espace. On obtient ces quantités à partir de la définition de l’onde plane incidente:
\(\frac{\partial {\left({\mathrm{u}}_{i}\right)}_{j}}{\partial {x}_{k}}={k}_{k}f'\left(k.\mathrm{x}-{C}_{m}t\right){k}_{j}\) avec \(m=S\) ou \(P\)
En ce qui concerne le terme d’impédance, il nous faut \(\frac{\partial {\mathrm{u}}_{i}}{\partial t}=-{C}_{m}\dot{f}\left(k.\mathrm{x}-{C}_{m}t\right)k\) , toujours avec \(m=S\) ou \(P\) .
On voit alors que la fonction importante pour un chargement par onde plane avec des éléments paraxiaux d’ordre 0 n’est pas le profil de l’onde \(f\) , mais sa dérivée, soit \(f'\) soit \(\dot{f}\) . L’onde étant plane, le front d’onde est caractérisé par les plans \(k.\mathrm{x}-{C}_{m}t=\mathit{cte}\) , d’où la relation: \(k.d\mathrm{x}={C}_{m}\mathit{dt}\) . Il y a donc l’équivalence suivante entre les deux dérivées de \(f\) : \(f'=\frac{1}{{C}_{m}}\dot{f}\) . On choisit de demander la fonction \(\dot{f}\) à l’utilisateur comme donnée du calcul.
Toutefois, dans le cas optionnel où on ajoute des rigidités distribuées le long de la frontière du domaine éléments finis , l’utilisateur doit en plus fournir le profil de l’onde \(f\) .
D’autre part, pour bien calculer le déphasage en temps dû au passage de l’onde, il faut également indiquer le point d’entrée du chargement par onde plane. Pour la prise en compte de l’onde réfléchie, celle-ci est activée si on donne la position du point de sortie. Pour plus de détails sur l’utilisation d’un chargement onde plane dans code_aster, il faut consulter la section dédiée de la commande AFFE_CHAR_MECA.
Le chargement onde plane inclinée#
Dans la prise en compte des ondes sismiques à partir d’une source profonde, les ondes résultantes peuvent prendre une direction inclinée. Or, dans la plupart des modélisations, seules les ondes planes de pression ou cisaillement dans la direction verticale sont prises en compte. Généralement, dans les études d’interaction sol-structure (ISS), la propagation d’une onde sismique en direction verticale est considérée comme l’approche standard. En effet, dans la plupart des cas, on s’attend à une verticalisation des ondes dans le bassin sédimentaire en lien avec la variation de la vitesse de propagation des ondes de cisaillement (\(v_S\)). Dans le cas où les ouvrages sont fondés sur du rocher (sans avoir donc l’effet de verticalisation par la loi de Snell dû à des couches superficielles à \(v_S\) plus faible) et pour des scénarios sismiques peu profonds, l’hypothèse de propagation à direction verticale est peu réaliste. La modélisation à propagation verticale permet de représenter assez simplement la propagation des ondes à partir de solutions analytiques ou modèles simplifiés. Toutefois, cette approche ne permet pas de prendre en compte certains déphasages dus au passage d’onde inclinée. C’est pourquoi, il est important d’adapter les modes de calcul de la propagation d’ondes verticales aux ondes inclinées. Dans cette partie, nous présentons :
Les éléments théoriques sur les ondes élastiques dans les solides isotropes
La mise en place d’une modélisation exclusivement pour des milieux homogènes non amortis. Pour ce cas spécifique, la formulation analytique est connue et permet d’imposer directement le chargement onde plane oblique aux frontières [bib5]
La mise en place d’une méthodologie généralisant l’application des ondes planes inclinées aux frontières des modèles dans les milieux hétérogènes et/ou amortis
Ondes planes élastiques dans les solides isotropes : réflexion et transmission aux interfaces#
Le phénomène de réflexion et transmission d’ondes aux interfaces est traité dans ce chapitre [bib6]. Par la suite, les ondes étudiées sont harmoniques de pulsation \(\omega\) et ,par souci de simplification, la dépendance temporelle \(e^{j\omega t}\) est omise dans les expressions des champs de déplacements. On étudie la réflexion d’une onde plane harmonique de pulsation \(\omega\) à la surface libre d’un demi-espace solide isotrope et homogène situé en \(\vector{x}_3<0\) (Fig. 152).
Fig. 152 Schéma d’une onde élastique se propageant dans un milieu homogène.#
Le solide est caractérisé par la masse volumique \(\rho\) et par les coefficients de Lamé \(\lambda\) et \(\mu\). Dans un premier temps, la nature de l’onde incidente n’est pas définie. Cette onde se propage dans le plan (\(\vector{x}_1\), \(\vector{x}_3\)) avec un angle d’incidence \(\theta_I\) et une vitesse \(v_I\). Sa direction de propagation est représentée par le vecteur \(n_I=(\sin{\theta_I}, 0, \cos{\theta_I})^T\). Le champ de déplacement en un point de coordonnées \(\vector{r}=(x_1, x_2, x_3)^T\) représentant l’onde incidente a pour expression
avec \(\vector{Q}_I=(Q_1,Q_2,Q_3 )^T\) et \(k_I=\frac{\omega}{v_I}\). On peut démontrer que trois ondes peuvent exister dans les solides isotropes : une onde longitudinale et deux ondes transverses. L’hypothèse est donc faite que ces trois ondes peuvent être générées lors de l’interaction de l’onde incidente avec la surface lorsque \(x_3=0\). L’onde longitudinale réfléchie se propage dans la direction \(\vector{n}_L=(\sin{\theta_L}, 0, -\cos{\theta_L})^T\) et les deux ondes transverses réfléchies se propagent dans la direction \(\vector{n}_T=(\sin{\theta_T},0,-\cos{\theta_T})^T\). Ainsi, le champ de déplacement des ondes réfléchies a pour expression :
avec:
où \(r_{IL}\), \(r_{IT_1}\) et \(r_{IT_2}\) sont des coefficients de réflexion en amplitude de déplacement. Le champ de déplacement longitudinal est irrotationnel, ce qui signifie que le vecteur polarisation \(\vector{Q}_L\) est colinéaire au vecteur \(\vector{n}_L\). De même, les deux champs de déplacement transverses sont à divergence nulle, imposant ainsi que les vecteurs polarisation \(\vector{Q}_{T_1}\) et \(\vector{Q}_{T_2}\) sont orthogonaux au vecteur \(\vector{n}_T\). Les vecteurs polarisation ont donc pour expression :
La première onde transverse est appelée par la suite onde « transverse verticale » et la seconde est appelé onde « transverse horizontale ». Compte tenu des lois de Snell-Descartes, on obtient :
Le déplacement \(u=u_I+u_T\) dans le demi-espace peut être écrit :
où \(f(x_1)=e^{-jk_I \sin{\theta_I}x_1 }=e^{-jk_L \sin{\theta_L}x_1 }=e^{-jk_T \sin{\theta_T}x_1 }\). Afin de déterminer les trois coefficients de réflexion en amplitude, il est nécessaire d’appliquer les conditions aux frontières. Ces conditions sont des conditions de contraintes nulles à la surface (x_3=0) du matériau :
L’application des conditions aux frontières correspond alors aux trois équations scalaires suivantes :
Il est donc nécessaire de calculer les expressions en \(x_3=0\) des contraintes \(\sigma_{13}\), \(\sigma_{23}\) et \(\sigma_{33}\) associées à chacune des quatre ondes en présence. Pour cela, on utilisera la loi de Hooke :
Donc, en utilisant les équations (976) et (973) on est en capacité de pouvoir déterminer les contraintes \(\sigma_{13}\), \(\sigma_{23}\) et \(\sigma_{33}\). Le système d’équations peut être donc réécrit comme un système de trois équations en trois inconnues :
avec \(Z_L=\rho v_L\) impédance associée à l’onde longitudinale et \(Z_T=\rho v_T\) impédance associée à l’onde transverse. Il vient immédiatement que la seconde équation de ce système est indépendante des deux autres équations. Les coefficients de réflexion \(r_IL\) et \(r_{IT_1}\) sont indépendants de la composante \(Q_2\) du vecteur polarisation de l’onde incidente et le coefficient de réflexion \(r_{IT_2}\) est indépendant des composantes \(Q_1\) et \(Q_3\). Il s’ensuit que si l’onde incidente est une onde longitudinale ou transverse verticale, il n’existe pas d’onde transverse horizontale réfléchie (voir Fig. 153 (a) et Fig. 153 (b)) et, de la même manière, si l’onde incidente est une onde transverse horizontale, il n’y a pas d’ondes longitudinale et transverse verticale réfléchies (voir Fig. 153 (c)). Dans ce second cas, on déduit de la seconde équation du système (977) que le coefficient de réflexion \(r_{IT_2}\) est égal à 1. Dans la suite, les cas d’une onde incidente longitudinale et d’une onde incidente transverse verticale sont traités séparément.
Fig. 153 Schéma de la réflexion d’ondes dans le cas d’une onde incidente (a) longitudinale (L), (b) transverse verticale (TV) et (c) transverse horizontale (TH).#
Onde incidente longitudinale#
Dans le cas des solides isotropes, les vitesses des ondes longitudinales et transverses sont les mêmes dans toutes les directions. Il s’en suit que les lenteurs, c’est-à-dire l’inverse des vitesses, sont également les mêmes dans toutes les directions. Les lois de Snell-Descartes peuvent alors être facilement interprétées en observant les surfaces des lenteurs présentées sur la Fig. 154. La lecture de ces surfaces démontre rapidement que l’angle \(\theta_T\) est toujours inférieur à l’angle d’incidence \(\theta_L\) et donc que les ondes transverses verticales réfléchies sont toujours des ondes planes propagatrices.
Fig. 154 Courbes des surfaces de lenteur pour une surface libre et une onde longitudinale incidente.#
Si l’onde incidente est une onde longitudinale, le vecteur polarisation est colinéaire au vecteur \(\vector{n}_L\). Le vecteur polarisation est donc du type \(\vector{Q}_L=(sin{\theta_L}, 0, cos{\theta_L})^T\). La première et la troisième équation du système (977) peuvent alors être mises sous la forme matricielle suivante :
avec \(\kappa=v_T/v_L\). La résolution de ce système conduit à l’expression des coefficients de réflexion suivante :
avec \(N = \kappa^2 \sin{2\theta_L} \sin{2\theta_T} + \cos^2{2\theta_T}\).
Onde incidente transverse verticale#
Les surfaces des lenteurs correspondant à une onde transverse verticale incidente sont présentées sur la Fig. 155. La lecture de ces surfaces démontre rapidement que l’angle \(\theta_T\) est toujours inférieur à l’angle d’incidence \(\theta_L\). Néanmoins, lorsque l’angle \(\theta_L=90°\), l’angle d’incidence est égal à un angle critique \(\theta_L^c=arcsin{(v_T/v_L)}\) traduisant la limite entre le domaine des ondes longitudinales réfléchies propagatives (\(\theta_T \leq \theta_L^c\)) et le domaine des ondes longitudinales réfléchies évanescentes (\(\theta_T \geq \theta_L^c\)).
Fig. 155 Courbes des surfaces de lenteur pour une surface libre et une onde transverse verticale incidente.#
Si l’onde incidente est une onde transverse verticale, le vecteur polarisation est perpendiculaire au vecteur \(\vector{n}_T\). Le vecteur polarisation est donc du type \(\vector{Q}_I=(\cos(\theta_T),0,-\sin(\theta_T))^T\). La première et la troisième équation du système (977) peuvent alors être mises sous la forme matricielle suivante :
La résolution de ce système conduit à l’expression des coefficients de réflexion suivante :
Utilisation dans code_aster#
La prise en compte d’éléments élastiques absorbants et le calcul de leur impédance nécessite une modélisation spécifique sur les frontières absorbantes :
en 2D avec la modélisation “D_PLAN_ABSO” sur les arêtes absorbantes.
en 3D avec la modélisation “3D_ABSO” sur les faces absorbantes.
La formulation de ces éléments étant assez rudimentaire pour pouvoir justement qu’ils soient assimilés à des discrets amortisseurs et donc être utilisés dans des analyses harmoniques, en conséquence, d’une part, on ne doit pas les bloquer lors d’analyses dynamiques, et d’autre part, une contrepartie est que la qualité de leur utilisation dépend de la qualité de la forme de la frontière. Un bon test de cette qualité s’inspire des cas tests SDLV120 et SDLV121 et peut être basé sur l’absorption plus ou moins totale au niveau de cette frontière d’un passage d’onde de déplacement imposé en haut de structure. Comme pour ces tests, on peut s’assurer que l’onde ne revient pas en regardant un degré de liberté de vitesse sur un nœud à proximité de la frontière.
En analyse harmonique avec l’opérateur DYNA_LINE_HARM [U4.53.11], on calcule au préalable un amortissement mécanique par l’option AMOR_MECA de l’opérateur CALC_MATR_ELEM [U4.61.01] et on le renseigne dans DYNA_LINE_HARM (mot clé MATR_AMOR).
En analyse transitoire, le calcul de l’amortissement mécanique est automatique avec les modélisations d’éléments absorbants dans l’opérateur DYNA_NON_LINE [U4.53.01]. Avec l’opérateur DYNA_LINE_TRAN [U4.53.02], on calcule au préalable cet amortissement mécanique de manière explicite par l’option AMOR_MECA de l’opérateur CALC_MATR_ELEM [U4.61.01] et on le renseigne par le mot clé MATR_AMOR.
Pour ces deux analyses, le calcul de la rigidité mécanique ajoutée est automatique avec les modélisations d’éléments absorbants quand on calcule systématiquement l’option RIGI_MECAquel que soit l’opérateur de calcul.
Bibliographie#
Crepel, M. (1983). Modélisation tridimensionnelle de l’interaction sol-structure par des éléments finis et infinis. Ecole Centrale de Paris. Phd thesis.
Cohen, P. & Jennings, C. (1983). Silent boundary methods for transient analysis. Comput. Methods for Transient Analysis. 301–357.
Modaressi, H. (1987). Modélisation numérique de la propagation des ondes dans les milieux poreux élastiques. Ecole Centrale de Paris. Phd thesis.
Engquist, B. & Majda, A. (1977). Absorbing boundary conditions for the numerical simulation of waves. Mathematics of Computation.
Humbert , N., & Péton, P. (2016). Comportement dynamique de barrage voûte. PFE Centrale Lyon.
Ewing, W. (1957). Elastic Waves in Layered Media. McGraw Hill Text.