Gestion de l’incompressibilité

État de l’art

Depuis les années 60, différentes stratégies ont été proposées afin de réduire ou éliminer le phénomène de verrouillage volumique et les oscillations de pression dans les solutions éléments finis. La technique d’intégration réduite, introduite par Zienkiewicz [Zienkiewicz et al.(1971)], est une des méthodes connues les plus anciennes. Elle consiste à intégrer la matrice de rigidité en utilisant une règle d’intégration d’ordre plus faible que celui nécessaire pour avoir une intégration exacte. Il s’est cependant avéré que cette méthode peut conduire à une matrice de rigidité singulière et introduire des modes à énergie nulle. La méthode d’intégration sélective, a été mise en place par Doherty et al [Doherty et al.(1969)]. 61 Gestion de l’incompressibilité Elle consiste à ne sous-intégrer que certains termes de la matrice de rigidité. Il s’agit de ceux qui sont liés à la contrainte d’incompressibilité (la partie dilatationnelle de la matrice). Les autres termes sont, quant à eux, intégrés de façon exacte. La sous-intégration sélective donne des résultats plus satisfaisants que la technique d’intégration réduite, bien que pouvant présenter de légères instabilités dans certaines situations. De plus, cette technique n’est pas suffisante pour éviter le phénomène de verrouillage [Ponthot(1995)]. Notons de plus que cette technique n’est valable que dans le cas isotrope, et que l’extension aux cas anisotropes ou orthotropes est compliquée à mettre en œuvre[Hughes(2012)]. De plus, elle ne fonctionne que pour des modèles pour lesquels il est possible de séparer le tenseur des déformations en partie déviatorique et volumétrique, comme c’est le cas pour des matériaux élastiques isotropes par exemple. Ainsi, ces techniques ont l’avantage d’être certes peu coûteuses en temps de calcul, mais dans le but de remédier aux problèmes d’apparition de modes parasites, des efforts ont été faits afin de les combiner avec des méthodes de stabilisation. Notons enfin que Malkus et Hughes [Malkus et Hughes(1978)] ont montré l’équivalence entre la méthode d’intégration réduite et d’intégration sélective et les formulations mixtes, sous certaines conditions, excluant notamment les cas axisymétriques [Ponthot(1995)]. C’est d’ailleurs grâce à eux que ces techniques sont passées du statut de simples « astuces » à celui de méthodes reconnues. Afin de palier aux limitations de la méthode de sous-intégration, la méthode B-bar a été mise en place par Hughes [Hughes(1980), Hughes(2012)] au début des années 80, dans le cadre des petites déformations. Il s’agit d’une méthode qui a l’avantage, tout comme les méthodes de sous-intégration, d’être basée sur une formulation en déplacements seulement, et qui fonctionne également dans le cas anisotrope. De plus, elle fonctionne avec tous types de lois de comportement, indépendamment de la décomposition en partie déviatorique et sphérique des équations de la loi de comportement, ce qui n’est pas le cas pour la méthode de sous-intégration sélective. Le principe de la méthode B-bar consiste à séparer la matrice de dérivation classique B (reliant les déformations aux déplacements dans la discrétisation éléments finis) en partie déviatorique et en partie dilatationnelle, et de modifier cette dernière, afin d’en réduire la contribution.

Formulations mixtes

Nous aurons besoin de définir les espaces des champs de déplacements, du champ de pression et éventuellement du champ de gonflement. Sauf mention contraire, on considérera les espaces : V = {u ∈ (H 1 (Ω0))ndim, u|Γu = u0 sur Ω0} avec u0 le déplacement imposé sur la frontière Γu (cf Chapitre 3). V 0 = {u ∈ (H 1 (Ω0))ndim, u|Γu = 0 sur Ω0} où ndim est la dimension de l’espace. P = G = L 2 (Ω0) où L 2 (Ω0) est l’ensemble des fonctions de carré intégrable, et H1 (Ω0) lest l’espace défini par : H 1 (Ω0) = {u ∈ L 2 (Ω0), ∀i ∈ {1, …, dim} ∂u ∂xi ∈ L 2 (Ω0) ndim} (4.8) Nous aurons besoin de discrétiser le problème. On considère un maillage Th du domaine Ω0. On note Ωh le domaine discret composé d’éléments Ωe ∈ Th, dont on précisera le type dans les applications. On s’intéressera aux espaces Vh et V 0 h définis ci-dessous ainsi qu’aux espaces Ph et Gh. Vh = {uh ∈ (C 0 (Ωh))ndim, uh|Ωe ∈ (P k (Ωe))ndim, uh|Γu = u0, ∀Ωe ∈ Th} 65 V 0 h = {uh ∈ (C 0 (Ωh))ndim, uh|Ωe ∈ (P k (Ωe))ndim, uh|Γu = 0, ∀Ωe ∈ Th} Ph = {ph ∈ C 0 (Ωh), ph|Ωe ∈ P m(Ωe), ∀Ωe ∈ Th} Gh = {θh ∈ C 0 (Ωh), θh|Ωe ∈ P n (Ωe), ∀Ωe ∈ Th} Notons que bien que les espaces Ph et Gh sont en apparence identiques, nous choisissons de les distinguer car l’ordre d’interpolation choisi pour le champ de pression ne sera pas nécessairement le même que celui choisi pour le champ de gonflement. 

Définition des fonctions de forme

On se propose de représenter les différents types d’éléments que nous considèrerons, avec la répartition des nœuds, ainsi que de donner l’expression des fonctions de forme pour les éléments linéaires. Pour les éléments quadratiques, nous nous contentons de rappeler l’expression des fonctions de forme pour les cas 2D, au vu du nombre important de nœuds en 3D (10 pour le tétraèdre P2, et 20 pour l’hexaèdre P2) et renvoyons à un livre traitant de la méthode des éléments finis. • Pour un élément triangle linéaire (P1):    N1 = 1 − ξ − η N2 = ξ N3 = η (4.9) • Pour un élément tétraèdre linéaire (P1): 66    N1 = 1 − ξ − η − ζ N2 = ξ N3 = η N4 = ζ (4.10) • Pour un élément quadrangle linéaire (Q1):    N1 = 1 4 (1 − ξ) (1 − η) N2 = 1 4 (1 + ξ) (1 − η) N3 = 1 4 (1 + ξ) (1 + η) N4 = 1 4 (1 − η) (1 + ξ) (4.11) • Pour un élément hexaèdre linéaire (Q1):    N1 = 1 8 (1 − ξ)(1 − η)(1 − ζ) N2 = 1 8 (1 + ξ)(1 − η)(1 − ζ) N3 = 1 8 (1 + ξ)(1 + η)(1 − ζ) N4 = 1 8 (1 − ξ)(1 + η)(1 − ζ) N5 = 1 8 (1 − ξ)(1 − η)(1 + ζ) N6 = 1 8 (1 + ξ)(1 − η)(1 + ζ) N7 = 1 8 (1 + ξ)(1 + η)(1 + ζ) N8 = 1 8 (1 − ξ)(1 + η)(1 + ζ) (4.12) • Pour un élément triangle quadratique (P2): 67    N1 = (1 − ξ − η) (1 − 2ξ − 2η) N2 = ξ (2ξ − 1) N3 = η (2η − 1) N4 = 4ξ (1 − ξ − η) N5 = 4ξη N6 = 4η (1 − ξ − η) (4.13) • Pour un élément tétraèdre quadratique (P2): • Pour un élément quadrangle quadratique (Q2):    N1 = 1 4 (1 − ξ) (1 − η) (−1 − η − ξ) N2 = 1 4 (1 + ξ) (1 − η) (−1 − η + ξ) N3 = 1 4 (1 + ξ) (1 + η) (−1 + η + ξ) N4 = 1 4 (1 − ξ) (1 + η) (−1 + η − ξ) N5 = 1 2 (1 − ξ 2 ) (1 − η) N6 = 1 2 (1 + ξ) (1 − η 2 ) N7 = 1 2 (1 − ξ 2 ) (1 + η) N8 = 1 2 (1 − ξ) (1 − η 2 ) (4.14) • Pour un élément hexaèdre quadratique (Q2):