v6.03.152 SSNP152 – Inclusion de deux couronnes

Résumé :

Ce test sert à évaluer les performances en ce qui concerne le traitement du contact entre deux structures avec compatibilité géométrique entre les surfaces maître et esclave à l’instant initial et différentes positions d’incompatibilité géométrique au cours du temps. Il permet de valider le traitement du contact avec la méthode continue en prenant en compte de grandes rotations.

On considère une structure composée de deux couronnes concentriques. La couronne externe est soumise à une pression uniforme sur son bord libre alors qu’on impose une rotation finie à la couronne intérieure. Leur rigidité, représentée par leurs modules de Young joue un rôle important dans l’évaluation de la valeur des déformations et des fluctuations de la pression de contact. On cherche aussi à savoir quels sont les effets de l’utilisation d’un maillage d’ordre supérieur.

Une solution analytique a été développée pour ce problème afin de valider les résultats numériques calculés. La validation de ce test porte sur les valeurs de la pression de contact.

Solution de référence

Nous développons ici une solution analytique pour le problème présenté ci-dessus. Cette solution sera développée sur l’hypothèse de petites déformations en considérant que les matériaux des couronnes sont isotropes, régis par une loi élastique linéaire sans variation de température.

Du fait des symétries du problème, la solution en déplacement du problème a la forme générique suivante :

\(u={u}_{r}(r,z).\underline{{e}_{r}}+{u}_{z}(r,z).\underline{{e}_{z}}\)

Cas 1 : déformations planes

En utilisant les symétries du problème et l’hypothèse d’invariance selon z des contraintes planes, la solution du problème prend la forme suivante:

\({u}_{r}={u}_{r}(r)\)

\({u}_{\theta}=0\) éq 2.1

\({u}_{z}=0\)

En utilisant l’équation de Lamé-Navier:

\(\begin{array}{c}(\lambda +\mu )\underline{\mathrm{grad}}(d(\underline{u}))+\mu \Delta \underline{u}+\underline{\mathrm{fd}}=\underline{0}\end{array}\) éq 2.2

où \(\underline{\mathrm{fd}}=\underline{0}\) sont les efforts volumiques nuls ici, et la formule du Laplacien:

\(\begin{array}{c}\Delta \underline{u}=\underline{\mathit{grad}}(div(\underline{u}))+\underline{\mathit{rot}}\underline{\mathit{rot}}(\underline{u})\end{array}\) éq 2.3

On peut écrire éq 2.2 sous la forme:

\(\begin{array}{c}(\lambda +2\mu )\underline{\mathit{grad}}(div(\underline{u}))+\mu \underline{\mathit{rot}}\underline{\mathit{rot}}(\underline{u})+\underline{\mathit{fd}}=\underline{0}\end{array}\) éq 2.4

soit encore en utilisant \(\underline{\mathrm{rot}}(\underline{u})=\overrightarrow{0}\text{et}\underline{\mathrm{fd}}=\overrightarrow{0}\text{et}u={u}_{r}(r).\underline{{e}_{r}}\) :

\(\begin{array}{c}d(\underline{u})=\frac{d}{\mathrm{dr}}{u}_{r}(r)+\frac{1}{r}{u}_{r}(r)\\ \underline{\mathrm{grad}}(d\underline{u})=\frac{d}{{d}_{r}}[\frac{1}{r}\frac{d}{{d}_{r}}(r{u}_{r}(r))].{\underline{e}}_{r}\\ \text{soit encore}(\lambda +2\mu )\frac{d}{{d}_{r}}[\frac{1}{r}\frac{d}{{d}_{r}}(r{u}_{r}(r))]=0\end{array}\) éq 2.5

en intégrant l’équation, on obtient pour les couronnes (extérieure et intérieure) la forme suivante du champ de déplacement:

\({u}_{r}={C}_{i}r+\frac{{D}_{i}}{r}\text{}{u}_{\theta}=0\text{}{u}_{z}=0\) éq 2.6

Pour déterminer \({C}_{i}\) et \({D}_{i}\) , il nous reste à imposer les conditions limites en pression et en déplacement. Pour cela, il faut d’abord calculer les déformations puis les contraintes associées au champ de déplacement.

Les déformations sont la partie symétrique du gradient des déplacements. On obtient:

\(\begin{array}{c}{\epsilon}_{\mathrm{rr}}={C}_{i}-\frac{{D}_{i}}{\mathrm{r²}}\\ {\epsilon}_{\theta \theta }={C}_{i}+\frac{{D}_{i}}{\mathrm{r²}}\\ {\epsilon}_{zz}={\epsilon}_{r\theta }={\epsilon}_{\mathrm{rz}}={\epsilon}_{\theta z}=0\end{array}\) éq 2.7

En appliquant la loi de Hooke:

\(\underline{\underline{\sigma}}=\lambda \mathrm{tr}(\underline{\underline{\epsilon}})\underline{\underline{1}}+2\mu \underline{\underline{\epsilon}}\) éq 2.8

on obtient la forme générale suivante pour les contraintes:

\(\begin{array}{c}{\sigma}_{\mathrm{rr}}=\frac{E}{1+\nu }(\frac{{C}_{i}}{1-2\nu }-\frac{{D}_{i}}{\mathrm{r²}})\\ {\sigma}_{\theta \theta }=\frac{E}{1+\nu }(\frac{{C}_{i}}{1-2\nu }+\frac{{D}_{i}}{\mathrm{r²}})\\ {\sigma}_{zz}=\frac{2\nu {\mathrm{EC}}_{i}}{(1+\nu )(1-2\nu )}\\ {\sigma}_{r\theta }={\sigma}_{\mathrm{rz}}={\sigma}_{\theta z}=\underline{0}\end{array}\) éq 2.9

En posant:

\(\begin{array}{c}{A}_{i}=\frac{{E}_{i}}{(1+{\nu}_{i})(1-2{\nu}_{i})}{C}_{i}\text{}{B}_{i}=\frac{{E}_{i}}{1+{\nu}_{i}}{D}_{i}\end{array}\) éq 2.10

les contraintes non nulles deviennent:

\(\begin{array}{c}{\sigma}_{\mathrm{rr}}={A}_{i}-\frac{{B}_{i}}{\mathrm{r²}}\\ {\sigma}_{\theta \theta }={A}_{i}+\frac{{B}_{i}}{\mathrm{r²}}\\ {\sigma}_{zz}=2\nu {A}_{i}\end{array}\) éq 2.11

Il ne nous reste plus qu’à calculer les valeurs de \({A}_{i}\) et \({B}_{i}\) pour chacune des couronnes. On notera \({\lambda}_{n}\) la pression de contact entre les deux couronnes telle que:

\(\begin{array}{}{\underline{\underline{\sigma}}}_{\mathrm{1rr}}({R}_{2}).(-{\underline{e}}_{r})={\lambda}_{n}{\underline{e}}_{r}\\ {\underline{\underline{\sigma}}}_{\mathrm{2rr}}({R}_{2}).{\underline{e}}_{r}=-{\lambda}_{n}{\underline{e}}_{r}\end{array}\) éq 2.12

avec les conditions aux limites:

\(\begin{array}{}{\underline{\underline{\sigma}}}_{\mathrm{1rr}}({R}_{1}).{\underline{e}}_{r}=-\mathrm{p.}{\underline{e}}_{r}\\ {\underline{\underline{\sigma}}}_{\mathrm{2rr}}({R}_{3}).(-{\underline{e}}_{r})=\underline{0}\end{array}\) éq 2.13

La condition de continuité sur le déplacement à l’interface entre les deux couronnes en contact donne de plus:

\(\begin{array}{}{u}_{r;1}(\mathrm{R2})={u}_{r;2}(\mathrm{R2})\end{array}\) éq 2.14

Nous avons donc 5 équations pour les 5 inconnues \(\begin{array}{}{A}_{1},{B}_{1},{A}_{2},{B}_{2,}{\lambda}_{n}\end{array}\) .

Le système des 4 premières équations nous permet d’obtenir:

\(\begin{array}{c}{A}_{1}=\frac{-p{R}_{1}^{2}+{\lambda}_{n}{R}_{2}^{2}}{{R}_{1}^{2}-{R}_{2}^{2}};{B}_{1}=(-p+{\lambda}_{n})\frac{{R}_{1}^{2}{R}_{2}^{2}}{{R}_{1}^{2}-{R}_{2}^{2}}\\ {A}_{2}=-{\lambda}_{n}\frac{{R}_{2}^{2}}{{R}_{2}^{2}-{R}_{3}^{2}};{B}_{2}=-{\lambda}_{n}\frac{{R}_{2}^{2}{R}_{3}^{2}}{{R}_{2}^{2}-{R}_{3}^{2}}\end{array}\) éq 2.15

et l’équation de continuité sur le déplacement permet enfin d’avoir la pression de contact:

\({\lambda}_{n}=\frac{2p{R}_{1}^{2}(1-{\nu}_{1})}{{R}_{1}^{2}+{R}_{2}^{2}(1-2{\nu}_{1})+\frac{{E}_{1}}{{E}_{2}}\frac{1+{\nu}_{2}}{1+{\nu}_{1}}\frac{{R}_{1}^{2}-{R}_{2}^{2}}{{R}_{2}^{2}-{R}_{3}^{2}}({R}_{2}^{2}(1-2{\nu}_{2})+{R}_{3}^{2})}\) éq 2.16

Cas 2 : contraintes planes

Nous supposons d’abord qu’il n’y a pas de contraintes suivant la direction perpendiculaire au plan des couronnes (\({\underline{\underline{\sigma}}}_{i}.{\underline{e}}_{z}=\underline{0}\) ). Les symétries du problème nous amènent à un champ de contraintes qui peut être écrit sous la forme:

\(\underline{\underline{{\sigma}_{i}}}(r,z)={\sigma}_{\mathrm{rr};i}{\underline{e}}_{r}\otimes {\underline{e}}_{r}+{\sigma}_{\theta \theta ;i}{\underline{e}}_{\theta}\otimes {\underline{e}}_{\theta}+{\sigma}_{r\theta ;i}{\underline{e}}_{r}\otimes {\underline{e}}_{\theta}\)

où l’indice \(i\) vaut 1 pour la couronne extérieure et 2 pour la couronne intérieure. En l’absence de forces volumiques et en considérant le problème quasi-statique (on néglige les effets de l’accélération), on a:

\(div\underline{\underline{\sigma}}=\underline{0}\)

En utilisant la solution en déplacement du problème générique:

\(u={u}_{r}(r,z).\underline{{e}_{r}}+{u}_{z}(r,z).\underline{{e}_{z}}\)

Le champ de déformation s’écrit alors:

\(\begin{array}{}{\epsilon}_{\mathrm{rr}}=\frac{\partial {u}_{r}(r,z)}{\partial r}\\ {\epsilon}_{\theta \theta }=\frac{{u}_{r}}{r}\\ {\epsilon}_{zz}=\frac{\partial {u}_{z}(r,z)}{\partial z}\\ {\epsilon}_{\theta z}={\epsilon}_{r\theta }=0\\ {\epsilon}_{\mathrm{rz}}=\frac{1}{2}(\frac{\partial {u}_{r}}{\partial z}+\frac{\partial {u}_{z}}{\partial r})\end{array}\) éq 2.17

Comme:

\(\begin{array}{}d\underline{\underline{\sigma}}=\underline{0}\text{devient}\frac{\partial {\sigma}_{\mathrm{rr}}}{\partial r}+\frac{{\sigma}_{\mathrm{rr}}-{\sigma}_{\theta \theta }}{r}=0\\ \text{avec}\\ {\sigma}_{\mathrm{rr}}=\frac{E}{1+\nu }({\epsilon}_{\mathrm{rr}}+\frac{\nu}{1-2\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta }+{\epsilon}_{zz}))\\ {\sigma}_{\theta \theta }=\frac{E}{1+\nu }({\epsilon}_{\theta \theta }+\frac{\nu}{1-2\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta }+{\epsilon}_{zz}))\\ {\sigma}_{zz}=\frac{E}{1+\nu }({\epsilon}_{zz}+\frac{\nu}{1-2\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta }+{\epsilon}_{zz}))\\ {\sigma}_{zz}={\sigma}_{z\theta }={\sigma}_{\mathrm{zr}}=0\\ \text{impliquant}\\ {\epsilon}_{\mathrm{zr}}=0=\frac{\partial {u}_{r}}{\partial z}+\frac{\partial {u}_{z}}{\partial r}\end{array}\) éq 2.18

Le fait que \({\sigma}_{zz}=0\) nous donne:

\(\begin{array}{}\frac{E}{1+\nu }({\epsilon}_{zz}+\frac{\nu}{1-2\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta }+{\epsilon}_{zz}))=0\\ {\epsilon}_{zz}(1-\nu )=-\nu ({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta })\\ {\epsilon}_{zz}=\frac{-\nu }{(1-\nu )}({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta })\end{array}\) éq 2.19

ce qui implique:

\(\begin{array}{}{\sigma}_{\mathrm{rr}}-{\sigma}_{\theta \theta }=\frac{E}{1+\nu }({\epsilon}_{\mathrm{rr}}-{\epsilon}_{\theta \theta })\\ \mathrm{tr}\underline{\underline{\epsilon}}=\frac{\nu -1}{\nu}{\epsilon}_{zz}+{\epsilon}_{zz}\\ \mathrm{tr}\underline{\underline{\epsilon}}=\frac{2\nu -1}{\nu}{\epsilon}_{zz}\end{array}\) éq 2.20

et permet d’écrire:

\(\begin{array}{}\mathrm{tr}\underline{\underline{\epsilon}}=\frac{1-2\nu }{1-\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta })\end{array}\) éq 2.21

On remplace alors \(\mathrm{tr}\underline{\underline{\epsilon}}\) par sa valeur dans \({\sigma}_{\mathrm{rr}}\) et \({\sigma}_{\theta \theta }\) pour obtenir:

\(\begin{array}{}{\sigma}_{\mathrm{rr}}=\frac{E}{1+\nu }({\epsilon}_{\mathrm{rr}}+\frac{\nu}{1-2\nu }\frac{1-2\nu }{1-\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta }))\\ {\sigma}_{\theta \theta }=\frac{E}{1+\nu }({\epsilon}_{\theta \theta }+\frac{\nu}{1-2\nu }\frac{1-2\nu }{1-\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta }))\\ \text{}\\ {\sigma}_{\mathrm{rr}}=\frac{E}{1+\nu }({\epsilon}_{\mathrm{rr}}+\frac{\nu}{1-\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta }))\\ {\sigma}_{\theta \theta }=\frac{E}{1+\nu }({\epsilon}_{\theta \theta }+\frac{\nu}{1-\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta }))\end{array}\) éq 2.22

L’équation d’équilibre devient ainsi :

\(\begin{array}{c}\frac{\partial {\sigma}_{\mathit{rr}}}{\partial r}+\frac{{\sigma}_{\mathit{rr}}-{\sigma}_{\theta \theta }}{r}=0\text{}\to \text{}\frac{\partial {\epsilon}_{\mathit{rr}}}{\partial r}+\frac{\nu}{1-\nu }(\frac{\partial}{\partial r}({\epsilon}_{\mathit{rr}}+{\epsilon}_{\theta \theta }))+\frac{{\epsilon}_{\mathit{rr}}-{\epsilon}_{\theta \theta }}{r}=0\end{array}\) éq 2.23

Comme nous avons:

\(\begin{array}{}{\epsilon}_{\mathrm{rr}}=\frac{\partial {u}_{r}(r,z)}{\partial r}\\ {\epsilon}_{\theta \theta }=\frac{{u}_{r}}{r}\end{array}\) éq 2.24

la substitution des contraintes par les déformations dans l’équation d’équilibre permet d’écrire finalement:

\(\begin{array}{}\frac{\partial \mathrm{²}{u}_{r}(r,z)}{\partial \mathrm{r²}}+\frac{\nu}{1-\nu }(\frac{\partial \mathrm{²}{u}_{r}(r,z)}{\partial \mathrm{²}r}+\frac{\partial {u}_{r}(r,z)}{\partial r})+\frac{1}{r}(\frac{\partial {u}_{r}(r,z)}{\partial r}-\frac{{u}_{r}(r,z)}{r})=0\\ \frac{\partial \mathrm{²}{u}_{r}(r,z)}{\partial \mathrm{r²}}+\frac{\nu}{1-\nu }(\frac{\partial \mathrm{²}{u}_{r}(r,z)}{\partial \mathrm{r²}}+\frac{1}{r}\frac{\partial {u}_{r}(r,z)}{\partial r}-\frac{1}{\mathrm{r²}}{u}_{r}(r,z))+\frac{1}{r}(\frac{\partial {u}_{r}(r,z)}{\partial r}-\frac{{u}_{r}(r,z)}{r})\\ \frac{\partial \mathrm{²}{u}_{r}(r,z)}{\partial \mathrm{r²}}+\frac{1}{r}\frac{\partial {u}_{r}(r,z)}{\partial r}-\frac{1}{\mathrm{r²}}{u}_{r}(r,z)=0\\ \frac{\partial}{\partial r}[\frac{1}{r}\frac{\partial}{\partial r}({\mathrm{r.u}}_{r}(r,z))]=0\end{array}\) éq 2.25

Par intégrations successives de éq 2.25 nous obtenons la forme suivante du champ \({u}_{r}(r,z)\) :

\(\begin{array}{}\frac{1}{r}\frac{\partial}{\partial r}({\mathrm{r.u}}_{r})=f(z)\\ \frac{\partial}{\partial r}({\mathrm{r.u}}_{r})=\mathrm{r.f}(z)\\ {\mathrm{r.u}}_{r}=\frac{\mathrm{r²}}{2}f(z)+g(z)\\ {u}_{r}=\frac{r}{2}f(z)+\frac{g(z)}{r}\end{array}\) éq 2.26

Soit:

\(\begin{array}{}{u}_{r}(r,z)=C(z)r+\frac{D(z)}{r}\end{array}\) éq 2.27

Comme nous avons:

\(\begin{array}{}\frac{\partial {u}_{z}}{\partial r}+\frac{\partial {u}_{r}}{\partial z}=0\\ {\epsilon}_{\mathrm{rr}}=\frac{{u}_{r}(r,z)}{\partial r}=C(z)-\frac{D(z)}{\mathrm{r²}}\\ {\epsilon}_{\theta \theta }=\frac{{u}_{r}}{r}=C(z)+\frac{D(z)}{\mathrm{r²}}\\ {\epsilon}_{zz}=\frac{\partial {u}_{z}}{\partial z}=\frac{\nu}{1-\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta })\\ {\epsilon}_{zz}=\frac{-2\nu }{1-\nu }C(z)\end{array}\) éq 2.28

On aboutit ainsi à:

\(\begin{array}{}\frac{\partial {u}_{z}}{\partial z}=\frac{-2\nu }{1-\nu }{C}_{1}z+{C}_{2}\end{array}\) éq 2.29

On peut donc écrire en intégrant \({\epsilon}_{zz}\) que \(\begin{array}{}{u}_{z}(r,z)=f(z)+g(r)\end{array}\) . En utilisant la première relation de éq 2.28 on obtient alors:

\(\begin{array}{}\frac{\partial {u}_{z}}{\partial r}=g'(r)=\frac{-\partial {u}_{r}}{\partial z}\\ g'(r)=-C'(z)r-\frac{D'(z)}{r}\\ C'(z)=\mathrm{cte}\to C(z)={C}_{1}z+{C}_{2}\\ D'(z)=\mathrm{cte}\to D(z)={D}_{1}z+{D}_{2}\\ g(r)=-{C}_{1}\frac{\mathrm{r²}}{2}-{D}_{1}\ln(r)+{C}_{0}\end{array}\) éq 2.30

On aboutit ainsi à:

\(\begin{array}{}\frac{\partial {u}_{z}}{\partial r}=g'(r)\end{array}\) éq 2.31

En intégrant les deux équations aux dérivées partielles éq 2.29 et éq 2.31 on obtient pour chacune des deux couronnes i=1 extérieure et i=2 intérieure:

\(\begin{array}{}{u}_{z}=\frac{-2\nu }{1-\nu }({C}_{1}\frac{\mathrm{z²}}{2}+{C}_{2}z)-{C}_{1}\frac{\mathrm{r²}}{2}-{D}_{1}\ln(r)+{C}_{0}\\ {u}_{z}(r,0)=g(r)+f(0)=0\to {C}_{0}={C}_{1}={D}_{1}=0\\ {u}_{z}(r,z)=\frac{-2\nu }{1-\nu }{C}_{2}z\end{array}\) éq 2.32

d’où l’on obtient:

\(\begin{array}{}{\sigma}_{\mathrm{rr}}=\frac{E}{1+\nu }({\epsilon}_{\mathrm{rr}}+\frac{\nu}{1-\nu }({\epsilon}_{\mathrm{rr}}+{\epsilon}_{\theta \theta }))\\ {\sigma}_{\mathrm{rr}}=\frac{E}{1+\nu }(\frac{\partial {u}_{r}(r,z)}{\partial r}+\frac{\nu}{1-\nu }(\frac{\partial {u}_{r}(r,z)}{\partial r}+\frac{{u}_{r}}{r}))\\ {\sigma}_{\mathrm{rr}}=\frac{E}{1+\nu }(\frac{1+\nu }{1-\nu }{C}_{i}+\frac{{D}_{i}}{\mathrm{r²}})\end{array}\) éq 2.33

En posant:

\({A}_{i}=\frac{E}{1+\nu }\frac{1+\nu }{1-\nu }{C}_{i}\text{}{B}_{i}=\frac{E}{1+\nu }{D}_{i}\) éq 2.34

le champ des contraintes s’écrit:

\(\begin{array}{c}{\sigma}_{\mathit{rr};i}(r)={A}_{i}-\frac{{B}_{i}}{{r}^{2}}\\ {\sigma}_{\theta \theta ;i}(r)={A}_{i}+\frac{{B}_{i}}{{r}^{2}}\\ {\underline{\underline{\sigma}}}_{i}(r)={A}_{i}\underline{\underline{1}}-\frac{{B}_{i}}{{r}^{2}}({\underline{e}}_{r}\otimes {\underline{e}}_{r}-{\underline{e}}_{\theta}\otimes {\underline{e}}_{\theta})\end{array}\) éq 2.35

Il ne nous reste plus qu’à calculer les valeurs de \({A}_{i}\) et \({B}_{i}\) pour chacune des couronnes. On notera \({\lambda}_{n}\) la pression de contact entre les deux couronnes telle que:

\(\begin{array}{c}{\underline{\underline{\sigma}}}_{1}({R}_{2}).(-{\underline{e}}_{r})={\lambda}_{n}{\underline{e}}_{r}\\ {\underline{\underline{\sigma}}}_{2}({R}_{2}).{\underline{e}}_{r}=-{\lambda}_{n}{\underline{e}}_{r}\end{array}\) éq 2.36

avec les conditions aux limites:

\(\begin{array}{c}{\underline{\underline{\sigma}}}_{1}({R}_{1}).{\underline{e}}_{r}=-\mathit{p.}{\underline{e}}_{r}\\ {\underline{\underline{\sigma}}}_{2}({R}_{3}).(-{\underline{e}}_{r})=\underline{0}\end{array}\) éq 2.37

on a:

\(\begin{array}{c}{A}_{1}=\frac{-p{R}_{1}^{2}+{\lambda}_{n}{R}_{2}^{2}}{{R}_{1}^{2}-{R}_{2}^{2}};{B}_{1}=(-p+{\lambda}_{n})\frac{{R}_{1}^{2}{R}_{2}^{2}}{{R}_{1}^{2}-{R}_{2}^{2}}\\ {A}_{2}=-{\lambda}_{n}\frac{{R}_{2}^{2}}{{R}_{2}^{2}-{R}_{3}^{2}};{B}_{2}=-{\lambda}_{n}\frac{{R}_{2}^{2}{R}_{3}^{2}}{{R}_{2}^{2}-{R}_{3}^{2}}\end{array}\) éq 2.38

par ailleurs la loi de comportement des matériaux permet d’écrire:

\(\begin{array}{c}{\underline{\underline{\epsilon}}}_{i}=\frac{1+{\nu}_{i}}{{E}_{i}}{\underline{\underline{\sigma}}}_{i}-\frac{{\nu}_{i}}{{E}_{i}}\mathit{Tr}({\underline{\underline{\sigma}}}_{i})\underline{\underline{1}}=\frac{1}{{E}_{i}}[{A}_{i}(1-{\nu}_{i})\underline{\underline{1}}-{B}_{i}(1+{\nu}_{i})\frac{1}{{r}^{2}}({\underline{e}}_{r}\otimes {\underline{e}}_{r}-{\underline{e}}_{\theta}\otimes {\underline{e}}_{\theta})]\\ \text{}=\frac{\partial {u}_{r;i}}{\partial r}{\underline{e}}_{r}\otimes {\underline{e}}_{r}+\frac{1}{2}(\frac{1}{r}\frac{\partial {u}_{r;i}}{\partial \theta }+\frac{\partial {u}_{\theta ;i}}{\partial r}-\frac{{u}_{\theta ;i}}{r})({\underline{e}}_{r}\otimes {\underline{e}}_{\theta}+{\underline{e}}_{\theta}\otimes {\underline{e}}_{r})+(\frac{1}{r}\frac{\partial {u}_{\theta ;i}}{\partial \theta }+\frac{{u}_{r;i}}{r}){\underline{e}}_{\theta}\otimes {\underline{e}}_{\theta}\end{array}\) éq 2.39

ce qui permet d’obtenir:

\({\underline{u}}_{i}=\frac{1}{{E}_{i}}[{A}_{i}(1-{\nu}_{i})r+{B}_{i}(1+{\nu}_{i})\frac{1}{r}]{\underline{e}}_{r}={f}_{\mathit{CP}}(r){\underline{e}}_{r}\) éq 2.40

La fonction \(f(r)\) du déplacement radial est donnée en fonction des propriétés des matériaux et de la pression \(p\) . Dans le cas de contraintes planes (MODELISATION = “C_PLAN”) on a:

\(\begin{array}{c}f({R}_{1})=\frac{1}{{E}_{1}}({A}_{1}(1-{\nu}_{1}){R}_{1}+\frac{{B}_{1}(1+{\nu}_{1})}{{R}_{1}})\\ f({R}_{3})=\frac{1}{{E}_{2}}({A}_{2}(1-{\nu}_{2}){R}_{3}+\frac{{B}_{2}(1+{\nu}_{2})}{{R}_{3}})\end{array}\) éq 2.41

Pour obtenir la valeur de \({\lambda}_{n}\) , on impose la continuité du vecteur déplacement en \(r={R}_{2}\) :

\(\begin{array}{c}{\underline{u}}_{1}({R}_{2})={\underline{u}}_{2}({R}_{2})\\ \frac{1}{{E}_{1}}[{A}_{1}(1-{\nu}_{1}){R}_{2}+{B}_{1}(1+{\nu}_{1})\frac{1}{{R}_{2}}]{\underline{e}}_{r}=\frac{1}{{E}_{2}}[{A}_{2}(1-{\nu}_{2}){R}_{2}+{B}_{2}(1+{\nu}_{2})\frac{1}{{R}_{2}}]{\underline{e}}_{r}\\ \frac{1}{{E}_{1}}[\frac{-p{R}_{1}^{2}+{\lambda}_{n}{R}_{2}^{2}}{{R}_{1}^{2}-{R}_{2}^{2}}(1-{\nu}_{1}){R}_{2}+(-p+{\lambda}_{n})\frac{{R}_{1}^{2}{R}_{2}^{2}}{{R}_{1}^{2}-{R}_{2}^{2}}(1+{\nu}_{1})\frac{1}{{R}_{2}}]{\underline{e}}_{r}\\ \text{}=\frac{1}{{E}_{2}}[-{\lambda}_{n}\frac{{R}_{2}^{2}}{{R}_{2}^{2}-{R}_{3}^{2}}(1-{\nu}_{2}){R}_{2}-{\lambda}_{n}\frac{{R}_{2}^{2}{R}_{3}^{2}}{{R}_{2}^{2}-{R}_{3}^{2}}(1+{\nu}_{2})\frac{1}{{R}_{2}}]{\underline{e}}_{r}\end{array}\) éq 2.42

soit encore:

\({\lambda}_{n}=\frac{2p{R}_{1}^{2}}{{R}_{1}^{2}(1+{\nu}_{1})+{R}_{2}^{2}(1-{\nu}_{1})+\frac{{E}_{1}}{{E}_{2}}\frac{{R}_{1}^{2}-{R}_{2}^{2}}{{R}_{2}^{2}-{R}_{3}^{2}}({R}_{2}^{2}(1-{\nu}_{2})+{R}_{3}^{2}(1+{\nu}_{2}))}\) éq 2.43

Remarque

Une fois que l’on a calculé la valeur des déplacements et de la pression de contact pour le cas en contraintes planes, on peut aisément calculer ces valeurs en déformations planes en remplaçant les valeurs du module d’Young et le coefficient de Poisson:

\({E}_{\mathit{DP}}=\frac{{E}_{\mathit{CP}}}{1-{\nu}_{\mathit{CP}}^{2}};{\nu}_{\mathit{DP}}=\frac{{\nu}_{\mathit{CP}}}{1-{\nu}_{\mathit{CP}}}\) éq 2.44

Où \({E}_{\mathit{CP}}\) et \({\nu}_{\mathit{CP}}\) prennent les valeurs \({E}_{i}\) et \({\nu}_{i}\) du §2.1.2 .

On obtient ainsi les valeurs des déplacements:

\({\underline{u}}_{i}=\frac{1+{\nu}_{i}}{{E}_{i}}[{A}_{i}(1-2{\nu}_{i})r+\frac{{B}_{i}}{r}]{\underline{e}}_{r}={f}_{\mathit{DP}}(r){\underline{e}}_{r}\) éq 2.45

et de la pression de contact:

\({\lambda}_{n}=\frac{2p{R}_{1}^{2}(1-{\nu}_{1})}{{R}_{1}^{2}+{R}_{2}^{2}(1-2{\nu}_{1})+\frac{{E}_{1}}{{E}_{2}}\frac{1+{\nu}_{2}}{1+{\nu}_{1}}\frac{{R}_{1}^{2}-{R}_{2}^{2}}{{R}_{2}^{2}-{R}_{3}^{2}}({R}_{2}^{2}(1-2{\nu}_{2})+{R}_{3}^{2})}\) éq 2.46

On remarque que les valeurs de \({A}_{i}\) et \({B}_{i}\) demeurent inchangées.

Calcul de la norme L2 de la pression de contact

Soit \({\Gamma}_{c}\) la surface de contact. La norme \({L}^{2}\) de la pression de contact est définie par:

\({\parallel \lambda \parallel }_{{L}^{2}}^{2}={\int}_{{\Gamma}_{c}}{\lambda}^{2}\mathit{dS}.\)

Dans le cas présent, \({\Gamma}_{c}\) l est le cercle de centre \(O\) et de rayon \({R}_{2}\) et la pression de contact est contact est uniforme.

On a donc:

\({\parallel \lambda \parallel }_{{L}^{2}}^{2}={\lambda}^{2}\mid {\Gamma}_{c}\mid =2\pi {\lambda}^{2}{R}_{2}.\)

Et finalement:

\({\parallel \lambda \parallel }_{{L}^{2}}=\mid \lambda \mid \sqrt{2\pi {R}_{2}}.\)

Valeurs testées

On teste la pression de contact à l’interface entre les 2 couronnes. Pour un nœud d’interface, la solution analytique est obtenue par les équations 1.3 et 1.5 en déformations planes et contraintes planes respectivement.

La valeur de la pression appliquée au bord à chaque pas de temps est donnée par la formule:

\(p(t)={p}_{0}{10}^{\frac{t}{10}-1,1},t\in \left\lbrace n\in \mathrm{\mathbb{Z}},1\le n\le 21\right\rbrace ,{p}_{0}=1,0\mathit{MPa}\) éq 2.47

pour \(t=1\) , la pression vaut \(0,1\mathit{MPa}\) et pour \(t=21\) elle vaut \(\mathrm{10MPa}\) . On remarque que les simulations ont été faites en négligeant les effets de l’accélération (commande STAT_NON_LINE) ce qui implique que les pas de temps possèdent des unités arbitraires.

Modélisation A

Il s’agit d’une modélisation en contraintes planes (C_PLAN), sans rotation de la couronne intérieure et où on impose un chargement augmentant exponentiellement avec le temps, le but étant de mesurer l’écart entre les valeurs calculées et celles obtenues analytiquement afin de mieux connaître le domaine de validité de la solution.

Les modules de Young \({E}_{1}={E}_{2}\) et les coefficients de Poisson \({\nu}_{1}={\nu}_{2}\) sont respectivement \(1.0E9\mathit{Pa}\) et \(0.2\) . La pression appliquée sur le bord de la couronne extérieure vaut \(1.0E6\mathit{Pa}\) et elle varie de 10% à 1000% de sa valeur au cours du temps.

La couronne extérieure définit la surface maîtresse.

Le maillage (Figure 3.2-a) comporte:

• 160 mailles de type SEG2;

• 240 mailles de type QUAD4.

Figure 3.2-1: Le maillage de la modélisation A

On calcule la pression de contact (LAGS_C) pour le nœud \(A\) de coordonnées \((\mathrm{0.6,0}.0)\) , celui qui à l’instant initial se trouve le plus à droite de l’interface entre les deux couronnes.

Identification	Référence
LAGS_C au nœud \(A\)	\(p(t)\)	ANALYTIQUE

Pour quelques pas de charge, on compare la valeur calculée à celle donnée par l’équation 1.5. La précisionest inférieure à 1%.

On teste également la norme \({L}^{2}\) de la pression de contact pour l’instant \(t=1\) , La valeur analytique est:

\({\parallel \lambda \parallel }_{{L}^{2}}=p(t=1)\sqrt{2\pi {R}_{2}}\approx 179780,177088\text{Mpa}\times {\text{m}}^{\frac{1}{2}}.\)

Identification	Type de référence	Valeur de référence	tolérance
Norme L2	Analytique	179780,177088	0, 052%

Puisque la simulation est faite sur l’hypothèse d’un comportement élastique, la configuration à un instant \(t\) ne dépend nullement des instants précédents: tout se passe comme si l’on effectuait plusieurs simulations indépendantes, chacune avec une valeur de chargement différent.

Modélisation B

La modélisation est identique à la modélisation A, mais dans ce cas on travaillera en déformations planes (D_PLAN) et, comme déjà mentionné ci-dessus, les coefficients de Poisson n’interviennent plus dans la solution.

La couronne extérieure définit la surface maîtresse.

Idem modélisation A.

On calcule la pression de contact (LAGS_C) pour le nœud \(A\) de coordonnées \((0.6,0.0)\) , celui qui à l’instant initial se trouve le plus à droite de l’interface entre les deux couronnes.

Identification	Référence
LAGS_C au nœud \(A\)	\(p(t)\)	ANALYTIQUE

Pour quelques pas de charge, on compare la valeur calculée à celle donnée par l’équation 1.5. La précisionest inférieure à 1%.

Figure 4.4-1: évolution de l’écart entre la solution analytique et la solution donnée par Code_Aster pour une intégration aux nœuds et aux points de Gauss.

Idem modélisation A.

Modélisation C

Par la suite, on ne travaillera que avec la modélisation en contraintes planes.

Les modules de Young et les coefficients de Poisson restent les mêmes. On fixe la valeur de la pression sur le bord de la couronne extérieure à (\(p=1.0E7\) ) et on impose une rotation de la couronne intérieure conformément à l’équation 1.6 pour un nombre \(N=10\) de pas de temps.

La couronne extérieure définit la surface maîtresse.

Idem modélisation A.

On calcule la pression de contact (LAGS_C) pour le nœud \(A\) de coordonnées \((0.6,0.0)\) , celui qui à l’instant initial se trouve le plus à droite de l’interface entre les deux couronnes. Les valeurs calculées sont comparées à la valeur obtenue d’après l’équation 1.5 pour une pression externe de \(p=1.0E7\) . La rotation de la couronne intérieure est appliquée. On regarde les variations de la pression extérieure au cours de cette rotation. On teste lorsque les maillages sont de nouveau en vis-à-vis.

Identification	Référence	Aster	tolérance
LAGS_C au nœud \(A\) Schéma d’intégration AUTO	\(\lambda =9.26E6\)	Analytique	\(0,6\text{\%}\)
LAGS_C au nœud \(A\) Schéma d’intégration GAUSS	\(\lambda =9.26E6\)	Analytique	\(0,8\text{\%}\)

Outre le calcul de la pression de contact, on est intéressé aussi par sa variation au cours du temps. L’incompatibilité des mailles maîtres et esclaves induit des fluctuations sur la valeur cette pression. La figure 5.4-1 montre (d’une façon un peu grossière) cet effet. Une façon de l’atténuer consiste à utiliser soit un maillage plus fin, soit des éléments d’ordre supérieur.

Figure 5.4-1: Fluctuation de la pression de contact due à l’incompatibilité des mailles

Pour estimer le plus correctement possible la pression de contact, lors de la rotation, lorsque les surfaces de contact maître et esclave ne sont plus compatibles, il faut utiliser les schémas d’intégration d’ordre le plus élevé possible avec un raffinement substantiel du maillage par rapport à la situation avec maillages compatibles des modélisations A et B. Le résultat est très nettement amélioré par le recours aux éléments quadratiques (cf. modélisation G).

Modélisation D

Idem modélisation C mais on joue avec les modules de Young et les coefficients de Poisson:

Couronne extérieure : \({E}_{1}=1.0E9\) , \({\nu}_{1}=0.3\)

Couronne intérieure : \({E}_{2}=1.0E8\) , \({\nu}_{2}=0.2\)

La couronne extérieure définit la surface maîtresse.

Idem modélisation A.

On calcule la pression de contact (LAGS_C) pour le nœud \(A\) de coordonnées \((0.6,0.0)\) , celui qui à l’instant initial se trouve le plus à droite de l’interface entre les deux couronnes. Les valeurs calculées sont comparées à la valeur obtenue d’après l’équation 1.5 pour une pression externe de \(p=1.0E7\) . La rotation de la couronne intérieure est appliquée. On regarde les variations de la pression extérieure au cours de cette rotation.

Identification	Référence	Aster	tolérance
LAGS_C au nœud \(A\) Schéma d’intégration AUTO	\(\lambda =2.418E6\)	Analytique	\(1,7\text{\%}\)
LAGS_C au nœud \(A\) Schéma d’intégration GAUSS	\(\lambda =2.418E6\)	Analytique	\(1,9\text{\%}\)

Lorsque les déplacements deviennent trop importants dû à une faible rigidité des couronnes, la valeur calculée s’écarte de celle calculée analytiquement, une fois que cette solution a été développée sur l’hypothèse de petits déformations et que la simulation a été faite en grands déplacements.

De l’autre coté, lorsque la rigidité des couronnes est trop importante, la fluctuation de la pression de contact augmente sensiblement. Ceci vient du fait que lorsqu’on impose un déplacement sur une structure, les contraintes auxquelles celle-ci est soumise peuvent être trop élevées (voir infinies) afin d’être compatibles avec les lois de la mécanique. Pour une structure très rigide, un petit déplacement est possible seulement avec des contraintes considérables.

Modélisation E

Idem modélisation C mais la couronne intérieure définit la surface maîtresse de façon à satisfaire la condition LBB (\(\mathit{P1}\) en contact et \(\mathit{P2}\) en déplacement).

Le maillage comporte:

• 80 mailles de type SEG2;

• 80 mailles de type SEG3;

• 120 mailles de type QUAD4 sur la couronne extérieure;

• 120 mailles de typeQUAD8sur la couronne intérieure.

On calcule la pression de contact (LAGS_C) pour le nœud \(A\) de coordonnées \((0.6,0.0)\) , celui qui à l’instant initial se trouve le plus à droite de l’interface entre les deux couronnes. Les valeurs calculées sont comparées à la valeur obtenue d’après l’équation 1.5 pour une pression externe de \(p=1.0E7\) . La rotation de la couronne intérieure est appliquée. On regarde les variations de la pression extérieure au cours de cette rotation. On teste lorsque les maillages ont un décalage maximal d’un demi-élément surfacique.

Identification	Référence	Aster	tolérance
LAGS_C au nœud \(A\) Schéma d’intégration AUTO	\(\lambda =9.26E6\)	Analytique	\(0,4\text{\%}\)

Modélisation F

Idem modélisation C.

La couronne extérieure définit la surface maîtresse de façon à satisfaire la condition LBB (\(\mathit{P1}\) en contact et \(\mathit{P2}\) en déplacement).

Nous travaillerons avec des éléments d’ordre 2 afin d’observer les effets de leurs utilisations lors de calculs. Pour la modélisation E, la couronne intérieure comportera des éléments d’ordre 2 et l’intérieure des éléments d’ordre 1. Pour la modélisation F, on inverse. Pour la modélisation G, tout le maillage aura des éléments d’ordre 2.

Le maillage comporte:

• 80 mailles de type SEG2;

• 80 mailles de type SEG3;

• 120 mailles de type QUAD4 sur la couronne intérieure;

• 120 mailles de typeQUAD8sur la couronne extérieure.

On calcule la pression de contact (LAGS_C) pour le nœud \(A\) de coordonnées \((0.6,0.0)\) , celui qui à l’instant initial se trouve le plus à droite de l’interface entre les deux couronnes. Les valeurs calculées sont comparées à la valeur obtenue d’après l’équation 1.5 pour une pression externe de \(p=1.0E7\) . La rotation de la couronne intérieure est appliquée. On regarde les variations de la pression extérieure au cours de cette rotation. On teste lorsque les maillages ont un décalage maximal d’un demi-élément surfacique.

Identification	Référence	Aster	tolérance
LAGS_C au nœud \(A\) Schéma d’intégration AUTO	\(\lambda =9.26E6\)	Analytique	\(0,9\text{\%}\)

Modélisation G

Idem modélisation C.

La couronne extérieure définit la surface maîtresse.

Le maillage quadratique comporte:

• 160 mailles de type SEG3;

• 240 mailles de typeQUAD8.

On calcule la pression de contact (LAGS_C) pour le nœud \(A\) de coordonnées \((0.6,0.0)\) , celui qui à l’instant initial se trouve le plus à droite de l’interface entre les deux couronnes. Les valeurs calculées sont comparées à la valeur obtenue d’après l’équation 1.5 pour une pression externe de \(p=1.0E7\) . La rotation de la couronne intérieure est appliquée. On regarde les variations de la pression extérieure au cours de cette rotation. On teste lorsque les maillages ont un décalage maximal d’un demi-élément surfacique.

Pour quelques pas de charge, on compare la valeur calculée à celle donnée par l’équation 1.5. La précisionest inférieure à 1,1%.

Identification	Référence	Aster	tolérance
LAGS_C au nœud \(A\) Schéma d’intégration GAUSS	\(\lambda =9.26E6\)	Analytique	\(1,5\text{\%}\)

Modélisation H

Idem modélisation G sauf qu’on utilise une modélisation sous-intégrée, “C_PLAN_SI”.

La couronne extérieure définit la surface maîtresse.

Le maillage quadratique comporte:

• 160 mailles de type SEG3;

• 240 mailles de typeQUAD8.

On calcule la pression de contact (LAGS_C) pour le nœud \(A\) de coordonnées \((0.6,0.0)\) , celui qui à l’instant initial se trouve le plus à droite de l’interface entre les deux couronnes. Les valeurs calculées sont comparées à la valeur obtenue d’après l’équation 1.5 pour une pression externe de \(p=1.0E7\) . La rotation de la couronne intérieure est appliquée. On regarde les variations de la pression extérieure au cours de cette rotation. On teste lorsque les maillages sont de nouveau en vis-à-vis.

Identification	Référence	Aster	tolérance
LAGS_C au nœud \(A\) Schéma d’intégration AUTO	\(\lambda =9.26E6\)	Analytique	\(1,5\text{\%}\)