r3.03.04 Efforts extérieurs de pression en grands déplacements#

Résumé:

Un chargement de pression en grands déplacements est un chargement suiveur. En employant des éléments de peau, on est amené à calculer, d’une part, un second membre dont le calcul est proche de celui en petits déplacements, et d’autre part, un terme de rigidité supplémentaire qui n’est, en général, pas symétrique.

Écriture continue du problème#

Éléments de cinématique en grandes transformations#

On considère un solide \(\domain\) soumis à des grandes déformations (voir figure). Soit \(\gradTransfor\) le tenseur gradient de la transformation, \(\phi\) faisant passer la configuration initiale \(\domainRefe\) à la configuration actuelle déformée \(\domainCurr\). On note \(\posiRefe\) la position d’un point dans \(\domainRefe\) et \(\posiCurr\) la position de ce même point après déformation dans \(\domainCurr\). \(\disp\) est alors le déplacement entre les deux configurations. On a donc:

(286)#\[\posiCurr = \posiRefe + \disp\]

Le tenseur gradient de la transformation s’écrit:

(287)#\[\gradTransfor = \frac{\partial \posiCurr}{\partial \posiRefe} = \tensTwoUnit + \gradVector{\disp}\]
../../../../_images/Shape15.gif

Fig. 130 Solide en grandes transformations#

Travail virtuel des efforts extérieurs de pression#

On considère une pression \(\presRefe\) normale à la surface dans la configuration de référence. Cette pression s’écrit \(\presCurr\) dans la configuration actuelle.

../../../../_images/Shape24.gif

Fig. 131 Configuration de référence et configuration actuelle#

Dans la configuration actuelle, le travail virtuel des efforts extérieurs de pression \(\work{\pres}\) s’écrit simplement (voir Fig. 131):

(288)#\[\work{\pres}(\disp) = {\int}_{\boundNCurr(\disp)} {-\pres \, \normalCurr \cdot \dispVirt \measBound }\]

De plus, on suppose dorénavant que la valeur de la pression ne dépend pas explicitement du déplacement mais seulement du point matériel d’application:

(289)#\[\presCurr(\posiCurr) = \presCurr(\funcTransfor (\posiRefe)) = \presRefe(\posiRefe)\]

Le coté suiveur de la force vient de la dépendance de la normale au déplacement. Dans ce cas, on peut alors exprimer le travail virtuel des efforts de pression dans la configuration de référence (changement de variable dans l’intégrale):

(290)#\[\work{\pres}(\disp) = {\int}_{\boundNRefe(\disp)} {-\presRefe \, \det \gradTransfor \, \left[ \inverseTranspose{\gradTransfor} \cdot \normalRefe \right] \cdot \dispVirt \left(\funcTransfor \left(\posiRefe\right)\right)\measBound }\]

Sur le plan pratique, on utilisera la formule (5002) pour calculer le travail des efforts de pression. Toutefois, la formule (5004) est la mieux adaptée à une dérivation par rapport au déplacement, dont on va voir la nécessité au paragraphe suivant.

Variation du travail virtuel des efforts extérieurs de pression#

Dans l’optique d’une résolution du problème d’équilibre de la structure par une méthode de Newton, on est amené à exprimer la variation du travail virtuel des efforts extérieurs de pression par rapport au déplacement, de manière similaire à ce qui a été fait pour le travail virtuel des efforts intérieurs dans [R5.03.01]. Le domaine d’intégration étant fixe dans l’expression (5004), la dérivation sous le signe somme est licite (cf. [Bib2] ):

(291)#\[\frac{\partial \work{\pres}\left(\disp\right)}{\partial \disp} = {\int}_{\boundNRefe(\disp)} {-\presRefe \, \frac{\partial}{\partial \disp} \left[\det \gradTransfor \, \inverseTranspose{\gradTransfor} \right] \cdot \delta \disp \cdot \normalRefe \cdot \dispVirt\measBound}\]

Nous décidons de choisir comme configuration de référence la configuration actuelle, pour laquelle \(\gradTransfor=\tensTwoUnit\). Ce choix conduit à une expression simple de la dérivée du terme entre crochets:

(292)#\[\frac{\partial}{\partial \disp}\left[ \det \gradTransfor \, \inverseTranspose{\gradTransfor} \right] \cdot \delta \disp = \divVector{\delta \disp} \, \tensTwoUnit - \gradVector{ \delta \disp}\]

Finalement, la variation du travail virtuel des efforts extérieurs de pression s’écrit dans la configuration actuelle:

(293)#\[\frac{\partial \work{\pres}(\disp)}{\partial \disp}\cdot\delta \disp\cdot\dispVirt = {\int}_{\boundNCurr(\disp)} {-\presCurr \, \left[ \divVector{\delta \disp} \, \tensTwoUnit - \gradVector{ \delta \disp} \right] \cdot \normalCurr \cdot \dispVirt \measBound}\]

Dans l’expression (5007) subsiste une difficulté. En effet, on s’attend à obtenir une grandeur essentiellement surfacique alors que l’intégrande fait apparaître des termes de dérivation normale à la surface. Autrement dit, il faut connaître l’expression des déplacements virtuels non seulement sur la surface du domaine mais aussi à l’intérieur de celui-ci (dans un voisinage de la surface pour pouvoir exprimer les dérivées normales). Cet inconvénient n’est pas anodin puisque pour calculer les termes élémentaires dus aux efforts surfaciques, on emploie des éléments de peau pour lesquels une variation normale n’a pas de sens.

Adoption d’un paramétrage curviligne de la surface#

Pour remédier au problème mentionné précédemment, il faut chercher à exprimer la relation (5006) à l’aide de grandeurs surfaciques uniquement. Pour cela, on a recours à des éléments de géométrie différentielle, [Bib1], dont on adopte les notations (en particulier, on adopte la convention de sommation des indices répétés où les indices grecs prennent les valeurs \(1\) et \(2\) tandis que les indices latins prennent les valeurs \(1\) à \(3\) ).

../../../../_images/Shape32.gif

Fig. 132 Paramétrage curviligne du voisinage de la surface soumise à la pression#

Soit \(({\theta}^{1},{\theta}^{2})\) un paramétrage admissible de la surface. Pour décrire le volume constitué d’un voisinage de cette surface, on lui adjoint une troisième variable, \({\theta}^{3}\), qui mesure la progression suivant la normale unitaire \(n\) en \(({\theta}^{1},{\theta}^{2})\). On a ainsi (voir Fig. 132):

(294)#\[\vector{OM}\left({\theta}^{1},{\theta}^{2},{\theta}^{3}\right) = \vector{OS}\left({\theta}^{1},{\theta}^{2})+{\theta}^{3} \,\normalCurr({\theta}^{1},{\theta}^{2}\right)\]

Avec ce choix de paramétrage, la base naturelle covariante \((\vectorBaseCV{1},\vectorBaseCV{2},\vectorBaseCV{3})\) s’écrit:

(295)#\[\vectorBaseCV{i} = \frac{\partial \vector{OM}}{\partial {\theta}^{i}}\]

Tandis que le tenseur métrique \(\metric\) vaut:

(296)#\[\begin{split}\metric = \vectorBaseCV{i} \times \vectorBaseCV{j} = \left[ \begin{array}{ccc} \tensTwoCmpCO{\metric}{1}{1}& \tensTwoCmpCO{\metric}{1}{2}& 0\\ \tensTwoCmpCO{\metric}{2}{1}& \tensTwoCmpCO{\metric}{2}{2}& 0\\ 0& 0& 1\end{array}\right]\end{split}\]

Dans ce paramétrage curviligne, l’intégrande (5007) a pour expression:

(297)#\[ {-\pres \, \tensTwoCmpCO{\metric}{i}{j} \, \vectorCmpCT{\normal}{i} \, \left[ \derivee{\vectorCmpCT{\delta \disp}{k}}{k} \, \vectorCmpCT{\dispVirt}{j} - \derivee{\vectorCmpCT{\delta \disp}{j}}{k} \, \vectorCmpCT{\dispVirt}{k} \right] }\]

Ce terme se simplifie considérablement. En effet, on peut déjà noter que lorsque \(j=k\) , le terme entre crochet est nul. En outre, dans le système curviligne adopté, les composantes contravariantes de \(\normal\) sont:

(298)#\[\vectorCmpCT{\normal}{1}=0, \quad \vectorCmpCT{\normal}{2}=0 \quad \vectorCmpCT{\normal}{3}=1\]

Enfin, en tenant compte de la forme particulière de \(\metric\) (c’est-à-dire \(\tensTwoCmpCO{\metric}{1}{3}=0\) , \(\tensTwoCmpCO{\metric}{2}{3}=0\) et \(\tensTwoCmpCO{\metric}{3}{3}=1\) ), la variation du travail (5007) s’écrit simplement:

(299)#\[\frac{\partial \work{\pres}(\disp)}{\partial \disp}\cdot\delta \disp \cdot \dispVirt = {\int}_{\boundNCurr(\disp)} { -\presCurr \, \left[ \derivee{\vectorCmpCT{\delta \disp}{\alpha}}{\alpha} \, \vectorCmpCT{ \dispVirt}{3} - \derivee{\vectorCmpCT{\delta \disp}{3}}{\alpha} \, \vectorCmpCT{\dispVirt}{\alpha} \right] \measBound}\]

Sur cette expression, on constate que seuls interviennent des opérateurs différentiels surfaciques (dérivation covariante par rapport à \({\theta}^{1}\) et \({\theta}^{2}\) seulement), ce qui est bien le but recherché. En introduisant la base contravariante \((\vectorBaseCT{1},\vectorBaseCT{2},\vectorBaseCT{3}=\normal)\), appelée aussi base duale et qui s’exprime à partir de la base covariante par:

(300)#\[vectorBaseCT{i} = {\left[\metric^{ij}\right]}^{-1} \, vectorBaseCV{j}\]

On peut s’affranchir des composantes curvilignes:

(301)#\[\frac{\partial \work{\pres}(\disp)}{\partial \disp}\cdot\delta \disp\cdot\dispVirt = {\int}_{\boundNCurr(\disp)} { -\presCurr \, \left[ \left( \frac{\partial \delta \disp}{\partial {\theta}^{\alpha}} \cdot \vectorBaseCT{\alpha} \right) \, \left(\dispVirt \cdot \normal \right) - \left( \frac{\partial \delta \disp}{\partial {\theta}^{\alpha}} \cdot \normal \right) \, \left( \dispVirt \cdot \vectorBaseCT{\alpha} \right) \right] \measBound}\]

C’est dorénavant l’expression (5012) qui sera utilisée pour calculer la variation du travail virtuel des efforts de pression.

Cas particulier d’une structure soumise à une pression interne ou externe constante#

Dans le cas particulier d’une pression constante dans une cavité (voir Fig. 133), on montre que les efforts de pression dérivent d’un potentiel \(\Xi\) qui n’est autre que le produit de la pression par le volume de la cavité. Ce résultat s’étend au cas d’une structure plongée dans un fluide à pression constante.

../../../../_images/Shape41.gif

Fig. 133 Structure sous pression interne ou externe constante#

On écrit ce potentiel:

(302)#\[\Xi = \pres \, {\int}_{\boundNCurr(\disp)} { \measDomain } = \pres \, {\int}_{\boundNRefe(\disp)} { \det \gradTransfor \measDomain }\]

A nouveau, on choisit comme configuration de référence la configuration actuelle. La variation de \(\Xi\) conduit alors bien au travail virtuel des efforts extérieurs de pression:

(303)#\[\frac{\partial \Xi }{\partial \disp}\cdot\dispVirt = \pres \, {\int}_{\domainCurr(\disp)} { \divVector(\dispVirt) \, \measDomain } = - {\int}_{\boundNCurr(\disp)} { \pres \, \dispVirt \cdot \normal \measBound} = \work{\pres} \, \dispVirt\]

Dans ce cas particulier, la variation du travail virtuel est aussi la seconde variation du potentiel \(\Xi\) , c’est-à-dire une forme bilinéaire symétrique :

(304)#\[\frac{\partial {W}^{\pres}(\disp)}{\partial \disp}\cdot\delta \disp\cdot\dispVirt=\frac{{\partial}^{2}\Xi (\disp)}{\partial {\disp}^{2}}\cdot\delta \disp\cdot\dispVirt\]

Discrétisation#

Options calculées#

Des éléments finis de peau (éléments surfaciques plongés dans un espace tridimensionnel) sont employés pour discrétiser les déplacements réels et virtuels intervenant dans des expressions surfaciques telles que (5004) et (5012). Ces dernières permettent d’exprimer respectivement le vecteur second membre et la matrice de rigidité dus à la pression, dont l’emploi par l’algorithme de STAT_NON_LINE est précisé en [R5.03.01]

On développe donc quatre options:

  1. RIGI_MECA_PRSU_R : matrice de rigidité pour une pression suiveuse comme constante réelle

  2. RIGI_MECA_PRSU_F : matrice de rigidité pour une pression suiveuse comme fonction réelle

  3. CHAR_MECA_PRSU_R : vecteur second membre pour une pression suiveuse comme constante réelle

  4. CHAR_MECA_PRSU_F : vecteur second membre pour une pression suiveuse comme fonction réelle

Ces options sont développées pour les éléments 3D, D_PLAN et AXIS. On peut appliquer une pression suiveuse normale mais pas de cisaillement tangent suiveur. La pression peut être une fonction réelle du temps ou une constante réelle.

Discrétisation des termes de géométrie différentielle#

On note \(\discVect{G_{\alpha}}\) la version discrétisée des vecteurs de la base covariante avec \({\alpha={1,2}}\).` Les deux premiers vecteurs de la base covariante \(\discVect{G_{\alpha={1,2}}}\) se calcule à partir du vecteur des déplacements discrétisés aux noeuds \(\discVect{U}\) et des dérivées des fonctions de forme \(\discMatr{B_{\alpha}}\) :

(305)#\[\discVect{G_{\alpha}} = \discMatr{B_{\alpha}} \, \discVect{U}\]

La normale \(\normal\) est calculée comme le produit vectoriel de ces deux premiers vecteurs \(\vector{g_{\alpha}}\) :

(306)#\[\normal = \frac{{\vector{g}}_{1} {\wedge } {\vector{g}}_{2}}{\Vert {\vector{g}}_{1} {\wedge }{\vector{g}}_{2}\Vert}\]

On peut aussi calculer le tenseur métrique discrétisé en composantes covariantes \(\discMatr{G_{\alpha \beta}}\) :

(307)#\[\discMatr{G_{\alpha \beta}} = \discVect{G_\alpha} \discVectLigne{{G}_{\beta}}\]

Et son jacobien:

(308)#\[\jacobTransfor = \det \discMatr{G_{\alpha \beta}}\]

On peut calculer la matrice métrique contravariante qu’on notera \(\discMatr{G^{\delta \gamma}}\)

(309)#\[\discMatr{G^{\delta \gamma}} = \inverse{\discMatr{ G_{\alpha \beta}}}\]

Et finalement extraire la base contravariante \(\discVectCmp{G}{\delta}\) :

(310)#\[\discVect{G^\delta} = \discMatr{G^{\delta \gamma}} \, \discVect{G_\alpha}\]

Vecteur des forces suiveuses#

Le calcul du travail virtuel des efforts de pression (5004) est en fait identique à celui effectué en petits déplacements, moyennant une réactualisation préalable de la géométrie. On part de l’expression des travaux virtuels:

(311)#\[\work{\pres}(\disp) = {\int}_{\boundNCurr(\disp)} {-\pres \, \normalCurr \cdot \dispVirt \measBound }\]

Sous forme discrétisée:

(312)#\[\work{\pres}(\disp) = \discVectLigne{\delta V} \discVect{{F}^{\pres}\left(\disp\right)}\]

La variation des déplacements s’écrit à partir des fonctions de forme:

(313)#\[\dispVirt = \discMatr{ \shapeFunc } \, \left\lbrace \delta V\right\rbrace\]

On a discrétisé tous les termes de géométrie différentielle dans le paragraphe précédent, il ne nous reste plus qu’à discrétiser l’intégrale en utilisant un schéma de Gauss avec les poids \(\quadWeight{\quadPointIndex}\). Sur une quantité quelconque \(A\), on aura

(314)#\[{\int}_{\boundNCurr(\disp)} {A \measBound } = \sum_{\quadPointIndex} {\onQuadPoint{A}{\quadPointIndex} \, \quadWeight{\quadPointIndex} }\]

La pression \(\pres\) donnée par AFFE_CHAR_MECA est localisée aux nœuds, on doit donc interpoler la pression \(\pres\) à partir de ses valeurs sur les nœuds qu’on note \(\onNode{\pres}{\nodeIndex}\) vers le point de Gauss de coordonnées \(\onQuadPoint{\xi}{\quadPointIndex}\) :

(315)#\[\onQuadPoint{\pres}{\quadPointIndex} = \sum_{\nodeIndex}{\onNode{\shapeFunc(\onQuadPoint{\xi}{\quadPointIndex})}{\nodeIndex} \, \onNode{\pres}{\nodeIndex} }\]

On note \(\onNode{\shapeFunc(\onQuadPoint{\xi}{\quadPointIndex})}{\nodeIndex}\) la valeur de la fonction de forme du noeud \(\nodeIndex\) au point de gauss de coordonnées \(\onQuadPoint{\xi}{\quadPointIndex}\). Finalement:

(316)#\[\discVect{{F}^{\pres}\left(\disp\right)} = - \sum_{{\quadPointIndex}}{ \onQuadPoint{\pres}{\quadPointIndex} \, \quadWeight{\quadPointIndex} \, \onQuadPoint{\jacobTransfor}{\quadPointIndex} \, \transpose{\discMatr{\onQuadPoint{ \shapeFunc} {\quadPointIndex} }} \, \discVect{\onQuadPoint{\normal} {\quadPointIndex} } }\]

Matrice des forces suiveuses#

La variation du travail virtuel des efforts de pression est exprimée par l’expression (5012)

(317)#\[\frac{\partial \work{\pres}(\disp)}{\partial \disp}\cdot\delta \disp\cdot\dispVirt = {\int}_{\boundNCurr(\disp)} { -\presCurr \, \left[ \left( \frac{\partial \delta \disp}{\partial {\theta}^{\alpha}} \cdot \vectorBaseCT{\alpha} \right) \, \left(\dispVirt \cdot \normal \right) - \left( \frac{\partial \delta \disp}{\partial {\theta}^{\alpha}} \cdot \normal \right) \, \left( \dispVirt \cdot \vectorBaseCT{\alpha} \right) \right] \measBound}\]

Il permet d’extraire la matrice tangente \(\discMatr{{K}^{\pres}\left(\disp\right)}\) des forces suiveuses:

(318)#\[\frac{\partial \work{\pres}(\disp)}{\partial \disp}\cdot\delta \disp\cdot\dispVirt = \discVectLigne{\delta V} \discMatr{{-K}^{\pres}\left(\disp\right)} \discVect{\delta U}\]

Le signe moins vient du fait que la contribution de la matrice est au premier membre.

Finalement:

(319)#\[\discMatr{{K}^{\pres}\left(\disp\right)} = \sum_{{\quadPointIndex}}{ \onQuadPoint{\pres}{\quadPointIndex} \, \quadWeight{\quadPointIndex} \, \onQuadPoint{\jacobTransfor}{\quadPointIndex} \, \transpose{\discMatr{\onQuadPoint{ \shapeDFunc} {\quadPointIndex} }} \left( \discVect{\onQuadPoint{G^\delta}{\quadPointIndex} } \, \discVectLigne{\onQuadPoint{\normal} {\quadPointIndex} } - \discVect{\onQuadPoint{\normal} {\quadPointIndex} } \, \discVectLigne{\onQuadPoint{G^\delta}{\quadPointIndex} } \right ) \, {\discMatr{\onQuadPoint{ \shapeFunc} {\quadPointIndex} }} }\]

avec \(\discMatr{\onQuadPoint{ \shapeDFunc} {\quadPointIndex} }\) la matrice des dérivées des fonctions de forme.

Choix de la matrice#

En général, la matrice \(\discMatr{{K}^{\pres}\left(\disp\right)}\) n’est pas symétrique (sauf cas particulier d’une structure soumise à une pression interne ou externe constante). On constate également en pratique, que pour de fortes variations de la géométrie (en utilisant un comportement hyper-élastique en grandes déformations comme ELAS_HYPER), le fait de symétriser cette matrice n’est pas une bonne stratégie (échec de la convergence). On décide donc de garder cette matrice non-symétrique, malgré le (léger) surcoût induit par la factorisation d’une telle matrice.

Bibliographie#

[Bib1]

Foundations of solid mechanics, Fung Y. C., Prentice Hall. 1965, pp 31-57.

[Bib2]

Calcul de la dérivée d’une grandeur par rapport à un fond de fissure par la méthode thêta, Mialon P., EDF - Bulletin de la Direction des Etudes et Recherches - SérieC - n° 3. 1988, pp1-28.