v6.02.125 SSNL125 - Traction d’un barreau fragile : endommagement à gradient#
Résumé:
Ce test permet la vérification des lois d’endommagement fragile à gradient ENDO_SCALAIRE modélisation A) et ENDO_FISS_EXP (modélisation B) dans une situation unidimensionnelle non homogène. De par son caractère 1D, ce problème admet une solution analytique qui exhibe deux régimes de couches limites: l’une de longueur finie (existence d’une frontière libre entre la zone endommagée et la zone saine) et l’autre de longueur infinie (elle s’étend jusqu’à la frontière de la pièce).
Solution de référence#
Nous considérons une des lois de comportement à gradient d’endommagement [R5.04.01] dont la densité d’énergie libre se présente sous la forme suivante :
\(\Phi (\varepsilon ,a)=A(a)w(\varepsilon )+\omega (a)+c/2{(\nabla a)}^{2}\) , où \(0\le a\le 1\) désigne la variable d’endommagement et \(w(\varepsilon )\) l’énergie de déformation élastique. Pour la loi ENDO_SCALAIRE, \(\omega (a)=\mathrm{ka}\) et le préfacteur d’endommagement s’écrit comme \(A(a)={(\frac{1-a}{1+\gamma a})}^{2}\) . La solution dépend alors de trois paramètres matériau \(\gamma\) , \(c\) , et \(k\) (voir [R5.03.18]).
Dans le cas général, deux équations aux dérivées partielles doivent être résolues: l’équation d’équilibre (\(\delta \Phi (\varepsilon ,a)/\delta \varepsilon =0\) ) et l’équation de comportement \(\delta \Phi (\varepsilon ,a)/\delta a=0\) . Obtenir une solution analytique s’avère généralement délicat, même pour des structures unidimensionnelles. Pour valider néanmoins ce modèle, on s’attache à un problème plus simple pour lequel l’équation d’équilibre n’a pas besoin d’être résolue, c’est-à-dire que le champ de déplacement est fixé partout. L’équation de comportement est alors pilotée par l’énergie de déformation élastique \(w\) connue en tout point de l’espace. Dans ce cas simplifié, la solution semi-analytique est obtenue par une intégration numérique d’une équation différentielle avec les conditions initiales de Cauchy.
Quant à la loi ENDO_FISS_EXP, elle est totalement équivalente à la loi ENDO_SCALAIRE dans une situation de traction uniaxiale. La solution analytique est donc identique à celle pour ENDO_SCALAIRE.
Caractérisation de la solution#
Plus précisément, on considère une barre dont une partie est astreinte à demeurer sans déformation tandis que l’autre est soumise à une déformation homogène. On étudie alors la couche limite d’endommagement qui se développe à l’interface de ces deux zones. L’équation différentielle de comportement est la suivante dans les zones où le critère est atteint, c’est-à-dire là où l’endommagement évolue:
\(c\frac{{\partial}^{2}a}{\partial {x}^{2}}-\frac{\partial A}{\partial a}w-k=0\) , où \(\frac{\partial A}{\partial a}(a)=-2(1+\gamma )\frac{(1-a)}{{(1+\gamma a)}^{3}}\) , éq 2.1-1
où \(x\) désigne la variable d’espace.
L’énergie de déformation élastique \(w\) est connue dans chacune des deux parties de la barre : \(w=0\) à gauche et \(w({\varepsilon}_{0})=\mathrm{const}\) à droite. La solution est donnée par une fonction continue dérivable \(a\in {C}^{1}(-{L}_{0,}{L}_{1})\) .
On dispose de trois conditions aux limites suivantes :
\(\frac{da}{dx}(-{L}_{0})=\frac{da}{dx}({L}_{1})=0\) et \(a(-{L}_{0})=0\) éq 2.1-2
ainsi que deux conditions d’interface :
\({\left[\mid a\mid \right]}_{x=0}=0{\left[∣\frac{da}{dx}∣\right]}_{x=0}=0\) éq 2.1-3
Au total, nous avons donc cinq conditions aux limites pour cinq inconnues. Malgré ce fait, le problème est «mal posé» pour pouvoir être résolu numériquement, car les conditions aux limites sont définies aux différents points de l’espace. Nous procédons dans la suite à l’identification des inconnues à l’aide des conditions aux limites enfin d’obtenir une équation différentielle ordinaire de type Cauchy-Lipschitz.
Résolution du problème dans la zone déchargée#
On suppose que l’extrémité gauche de la barre \(x=-{L}_{0}\) se situe suffisamment loin de l’interface \(x=0\) pour que l’hypothèse d’une couche limite de longueur finie dans la zone déchargée reste valide. Notons alors \({b}_{0}>0\) la longueur ( a priori dépendante du niveau de déformation \({\varepsilon}_{0}\) dans la zone chargée) sur laquelle se développe cette couche limite.
Dans la partie \(x\in \left[{L}_{0},-{b}_{0}\right]\) l’endommagement n’évolue pas et reste nul :
\(a(x)=0\) sur \(\left[-{L}_{0},-{b}_{0}\right]\) éq 2.2-1
Dans la partie \(x\in \left[-{b}_{0},0\right]\) , le critère est atteint et le champ d’endommagement évolue selon l’équation [éq 2.1-1], qui se simplifie considérablement dans cette partie de la barre puisque l’énergie de déformation élastique y est nulle :
\(c\frac{{d}^{2}a}{{\mathrm{dx}}^{2}}=k\) sur \(\left[-{b}_{0},0\right]\) éq 2.2-2
Par définition de \({b}_{0}\) , on dispose des conditions suivantes :
\(a(-{b}_{0})=0\) et \(\frac{\mathrm{da}}{\mathrm{dx}}(-{b}_{0})=0\) éq 2.2-4
On peut alors exprimer analytiquement la variable d’endommagement en intégrant [éq 2.2-2] :
\(a(x)=\frac{k}{\mathrm{2c}}{(x+{b}_{0})}^{2}\) sur \(\left[-{b}_{0},0\right]\) éq 2.2-5
On connaît alors les valeurs de l’endommagement et de sa dérivée à l’interface \(x=0\) , en fonction de l’inconnue \({b}_{0}\) :
\(a({0}^{-})=\frac{{\mathrm{kb}}_{0}^{2}}{2c}\) et \(a'({0}^{-})=\frac{{\mathrm{kb}}_{0}}{c}\) éq 2.2-6
et l’endommagement à l’interface \({a}_{0}\) étant strictement compris entre \(0\) et \(1\) , on en déduit l’encadrement suivant de la longueur \({b}_{0}\) :
\(0<{a}_{0}<1\Rightarrow 0<{b}_{0}<\sqrt{\frac{\mathrm{2c}}{k}}\) éq 2.2-7
Il suffit donc de prendre la partie gauche de la barre plus longue que \(\sqrt{\mathrm{2c}/k}\) pour être certain d’avoir le profil d’endommagement confiné à gauche.
Résolution du problème dans la zone chargée#
Cette partie de la barre voit une déformation homogène \({\varepsilon}_{0}>0\) associée à une énergie de déformation élastique \({w}_{0}\) . Dans cette zone, la couche limite n’est plus bornée et s’étend asymptotiquement vers la réponse homogène. On supposera donc que l’extrémité droite de la barre \(x={L}_{1}\) se situe suffisamment loin de l’interface \(x=0\) pour pouvoir faire l’approximation d’une dérivée nulle pour le champ d’endommagement «à l’infini». Au voisinage de l’extrémité droite de la barre, l’équation de comportement [éq 2.1-1] se réduit alors à :
\(\frac{-\partial A}{\partial a}({a}_{\infty}){w}_{0}=k\) éq 2.3-1
où \({a}_{\infty}\) désigne la valeur asymptotique du champ d’endommagement. Compte tenu de l’expression de \(\partial A/\partial a\) ( cf. [éq 2.1-1]), on peut alors paramétrer le niveau de chargement par la seule valeur \({a}_{\infty}\) :
pour \({a}_{\infty}\in \left]0,1\right[,{w}_{0}=\frac{k}{2(1+\gamma )}\frac{{(1+\gamma {a}_{\infty})}^{3}}{(1-{a}_{\infty})}\) éq 2.3-2
D’autre part, l’équation différentielle non linéaire [éq 2.1-1] peut s’écrire sous la forme suivante :
\(\frac{c}{2}2a'a''=a'(A'(a)w-k)\) éq 2.3-3
elle admet donc une intégrale première :
\(\forall x\in \left[0,{L}_{1}\right],{\left[\frac{c}{2}{(\frac{\mathrm{da}}{\mathrm{dx}}(s))}^{2}\right]}_{0}^{x}={\left[A(a(s))+ka(s)\right]}_{0}^{x}\) éq 2.3-4
On rappelle que la condition d’interface [éq 2.1-3] impose un raccord \({C}^{1}\) en \(x=0\) , en tenant compte des expressions liant \(a({0}^{+})\) et \(a'({0}^{+})\) à la longueur \({b}_{0}\) de la couche limite on peut alors écrire :
\({a}_{0}=a({0}^{+})=\frac{{\mathrm{kb}}_{0}^{2}}{2c}\) et \(a{'}_{0}=a'({0}^{+})=\sqrt{\frac{2k}{c}{a}_{0}}\) éq 2.3-5
De plus, les conditions aux limites [éq 2.1-2] assurent la nullité de la dérivée de l’endommagement en \(x={L}_{1}\) , en faisant l’approximation \(a({L}_{1})={a}_{\infty}\) , l’intégrale première [éq 2.3-4] évaluée en \(x={L}_{1}\) revient à :
\(A({a}_{0}){w}_{0}+k{a}_{0}-\frac{c}{2}a{'}_{0}^{2}=A({a}_{\infty}){w}_{0}+k{a}_{\infty}\) éq 2.3-6
L’équation [éq 2.3-6] se simplifie d’après [éq 2.3-5], elle se réduit alors au trinôme du second degré suivant d’inconnue \({a}_{0}\) :
\((f({a}_{\infty}){\gamma}^{2}-1){a}_{0}^{2}+2(f({a}_{\infty})\gamma +1){a}_{0}+(f({a}_{\infty})-1)=0\) où \(f({a}_{\infty})=A({a}_{\infty})+\frac{k{a}_{\infty}}{{w}_{0}}\) éq 2.3-7
On dispose donc d’une relation simple entre le paramètre de chargement \({a}_{\infty}\) et la valeur \({a}_{0}\) de l’endommagement ainsi que de sa dérivée \(a{'}_{0}\) à l’interface.
La résolution de l’EDO non linéaire [éq 2.1-1] sur \(\left[0,{L}_{1}\right]\) a été réalisée en la ramenant à un système différentiel d’ordre 1, intégré numériquement avec un schéma de Runge-Kutta à l’ordre 4 (en espace). La mise en œuvre de ce schéma ne nécessite que la connaissance des conditions initiales \(({a}_{0,}a{'}_{0})\) associées à chaque paramètre de chargement \({a}_{\infty}\) , ce qui est assuré par la relation [éq 2.3-7]
Résolution du problème dans la barre entière#
Afin de calculer la solution de référence sur toute la barre, on adopte le schéma suivant :
choix d’un paramètre de chargement \({a}_{\infty}\in \left]0,1\right[\)
Résolution dans la partie chargée :
détermination avec [éq 2.3-2] de l’énergie de déformation élastique \({w}_{0}\) associée
détermination avec [éq 2.3-7] du jeu de conditions aux limites \(({a}_{0,}a{'}_{0})\)
intégration numérique avec un schéma de Runge-Kutta à l’ordre 4
Résolution dans la partie déchargée :
détermination avec [éq 2.3-7] et [éq 2.3-5] de la taille \({b}_{0}\) de la couche limite dans la zone déchargée
obtention avec [éq 2.2-5] du champ analytique \(a(x)\) dans cette zone
Application numérique#
Pour l’élasticité, l’écrouissage et le paramètre non local, on adopte les caractéristiques suivantes:
\(\begin{array}{cc}E={3.10}^{4}\mathrm{MPa}& {\sigma}^{y}=3.\mathrm{MPa}\\ \nu =0.& \gamma =4.\end{array}c=1.875N\) éq 2.5-1
Ces choix conduisent à une énergie de déformation élastique seuil :
\({w}_{y}=1.5{10}^{-4}\mathrm{MPa}\) éq 2.5-2
Résultats de référence#
La solution de référence est obtenue en prenant une barre de longueur \(-{L}_{0}=-125\mathrm{mm}\) et \({L}_{1}=250\mathrm{mm}\) . On examine la valeur du champ d’endommagement \(a\) pour trois niveaux de chargement et en deux endroits, l’un dans la zone déchargée, l’autre dans la zone chargée.
\({\varepsilon}_{0}\) |
\({w}_{0}\) (MPa) |
\(a(x=-7.5\mathrm{mm})\) |
\(a(x=+7.5\mathrm{mm})\) |
2,70000000000000E-04 |
1,09350000000000E-003 |
1,93274688119012E-02 |
1,41846675324338E-01 |
7,34846922834953E-04 |
8,10053999999999E-003 |
1,39107889370765E-01 |
3,80008828951670E-01 |
1,10464444958548E-02 |
1,83060308787201E+000 |
6,14240950943351E-01 |
9,77312427816067E-01 |
Table2.6-1 - Résultats de référence
Figure 2.6-1: Profil d’endommagement pour la solution de référence
Modélisations A et B#
Caractéristiques de la modélisation et du maillage#
Il s’agit d’une modélisation axisymétrique (AXIS_GRAD_VARI). La géométrie correspondante est un rectangle, c’est-à-dire que le barreau est disposé de manière verticale et sa section (sans influence) est circulaire.
Le maillage est constitué d’un seul élément selon le rayon. Selon l’axe, les éléments ont une taille de \(2.5\mathrm{mm}\) . Le maillage ainsi engendré est finalement constitué de 150 éléments quadrangulaires à 8 nœuds.
Grandeurs testées et résultats#
On valide la modélisation et l’algorithme d’intégration de lois non locales en examinant le niveau d’endommagement (variable interne \(\mathrm{V1}\) ) aux différents niveaux de chargement et aux différents lieux géométriques listés dans la [Table 2.6-1]. Les résultats sont réunis dans l’extrait du fichier de résultat ci‑dessous.
Identification |
Instant |
Référence |
\(\mathrm{V1}(x=-7.5)\) |
2,70000000000000E-04 |
1,93274688119012E-02 |
\(\mathrm{V1}(x=+7.5)\) |
2,70000000000000E-04 |
1,41846675324338E-01 |
\(\mathrm{V1}(x=-7.5)\) |
7,34846922834953E-04 |
1,39107889370765E-01 |
\(\mathrm{V1}(x=+7.5)\) |
7,34846922834953E-04 |
3,80008828951670E-01 |
\(\mathrm{V1}(x=-7.5)\) |
1,10464444958548E-02 |
6,14240950943351E-01 |
\(\mathrm{V1}(x=+7.5)\) |
1,10464444958548E-02 |
9,77312427816067E-01 |
Synthèse des résultats#
On note un très bon accord entre la modélisation et la solution analytique.