r7.02.01 Taux de restitution de l’énergie en thermo-élasticité linéaire

Résumé :

On présente le calcul du taux de restitution de l’énergie \(G\) par la méthode \(\theta\) en 2D ou en 3D pour un problème thermo-élastique linéaire. On explique comment le champ \(\theta\) est introduit et comment le taux de restitution d’énergie est calculé avec l’opérateur CALC_G.

Des études mécano-fiabilistes d’évaluation de probabilité d’amorçage de la rupture requièrent, en plus, sa dérivée par rapport à une variation de domaine pilotée par un autre champ. On détaille l’implantation de cette option dans le code.

Le calcul de \(G\) est disponible aussi bien pour une fissure maillée (éléments finis classiques) que pour une fissure non maillée (éléments finis enrichis: méthode X-FEM). Le calcul de sa dérivée n’est quant à lui disponible qu’en 2D, pour une fissure maillée.

Discrétisation du taux de restitution d’énergie

Méthode \(\theta\) en dimension 2

On rappelle que le taux de restitution d’énergie \(G\) est solution de l’équation variationnelle :

\[\int_{{\Gamma}_{0}} G(s) \theta (s) \cdot m(s) \, ds = G(\theta), \quad \forall \theta \in \Theta\]

où :

• \(m\) est la normale unitaire au fond de fissure \({\Gamma}_{0}\), située dans le plan tangent à \(\partial \Omega\) et orientée vers l’intérieur de \(\Omega\),

• \(\Theta = \lbrace \mu \text{ tels que } \mu \cdot n = 0 \text{ sur } \partial \Omega \rbrace\).

En dimension 2, le fond de fissure \({\Gamma}_{0}\) se ramène à un point \(M_{0}\), et on peut choisir un champ \(\theta\) unitaire au voisinage de ce point, de telle sorte que : \(G(M_{0}) = G(\theta)\)

Fig. 300 Figure 2.1-a : Fond de fissure en 2D

Méthode \(\theta\) en dimension 3

La dépendance de \(G(\theta )\) vis-à-vis du champ \(\theta\) sur le fond de fissure est plus complexe. Le champ scalaire \(G(s)\) peut être discrétisé sur une base que nous noterons \({({p}_{j}(s))}_{1\le j\le N}\).

Fig. 301 **Figure 2.2-a : Discrétisation du fond de fissure en 3D (abscisse curviligne)

Soit \({G}_{j}\) les composantes de \(G(s)\) dans cette base :

\[G(s) = \sum_{j=1}^{N} G_{j} p_{j}(s).\]

On se donne également une base de fonctions tests pour les champs \(\theta\), de longueur \(P\) :

\[\tilde{\Theta}=\lbrace {\theta}^{i} \in \Theta, \quad i=1, \dots, P \rbrace.\]

\(\tilde{\Theta}\) est un sous-ensemble de longueur finie de l’ensemble \(\Theta\).

\(G(s)\) étant solution de l’équation variationnelle :

\[\int_{{\Gamma}_{0}} G(s) \theta (s) \cdot m(s) \, ds = G(\theta), \quad \forall \theta \in \Theta,\]

les coefficients \({G}_{j}\) vérifient alors en particulier :

\[\int_{{\Gamma}_{0}} \left( \sum_{j=1}^{N} G_{j} p_{j}(s) \right) \bar{\theta}^{i} \cdot m(s) \, ds = G(\theta^{i}), \quad \forall i \in [1,P],\]

soit :

\[\sum_{j=1}^{N} \left( \int_{{\Gamma}_{0}} p_{j}(s) \bar{\theta}^{i} \cdot m(s) \, ds \right) G_{j} = G(\theta^{i}), \quad \forall i \in [1,P].\]

Les coefficients \({G}_{j}\) peuvent donc être déterminés en résolvant le système linéaire à \(P\) équations et \(N\) inconnues :

\[\begin{split}\begin{cases} \sum_{j=1}^{N} a_{ij} G_{j} = b_{i}, \quad i = 1, P, \\ a_{ij} = \int_{{\Gamma}_{0}} p_{j}(s) \bar{\theta}^{i}(s) \cdot m(s) \, ds, \\ b_{i} = G(\theta^{i}). \end{cases}\end{split}\]

Ce système a une solution si on choisit \(P\) champs \(\theta^{i}\) indépendants tels que \(P \geq N\). Il peut comporter plus d’équations que d’inconnues, auquel cas il est résolu au sens des moindres carrés.

Remarque :

On se donne en réalité une base de fonctions tests pour la trace des champs \(\theta\) sur le fond de fissure \({\Gamma}_{0}\), notée \(\bar{\theta}\) et définie par \(\bar{\theta}(s) = {\theta}_{\mid {\Gamma}_{0}}(s)\). La donnée de la valeur de \(\bar{\theta}^{i}\) sur le fond de fissure \({\Gamma}_{0}\) suffit alors pour identifier le champ \(\theta^{i}\) sur tout le domaine.

Choix de la discrétisation de G en dimension 3

Description des espaces d’approximations

Remarque :

• En dimension 2, il n’y a pas de problème car en choisissant un champ \(\theta\) unitaire au voisinage du fond de fissure, on obtient la relation \(G = G(\theta)\). Le taux de restitution de l’énergie est donc indépendant du champ \(\theta\).

En dimension 3, la dépendance de \(G(\theta )\) vis-à-vis du champ \(\theta\) sur le fond de fissure est plus complexe. On calcule le taux de restitution d’énergie local \(G(s)\), solution de l’équation variationnelle :

\[\int_{{\Gamma}_{0}} G(s) \theta (s) \cdot m(s) \, ds = G(\theta ), \quad \forall \theta \in \Theta.\]

On a choisi deux familles de bases de fonctions tests pour \(\bar{\theta}\) et de décomposition pour \(G(s)\) :

• Les polynômes de Legendre \({\gamma}_{j}(s)\), de degré \(j\) avec \(0 \leq j \leq \mathrm{Deg}_{\max}\).

• Les fonctions de forme associées au nœud \(k\) de \({\Gamma}_{0}\), notées \({\phi }_{k}(s)\), avec \(1\leq k\leq \mathrm{NNO}\) (nombre de nœuds de \({\Gamma}_{0}\)), de degré un pour les éléments linéaires et deux pour les éléments quadratiques.

Rappelons que les polynômes de Legendre constituent une famille orthogonale non normée. Ils sont obtenus par la relation de récurrence :

\[(n+1) P_{n+1}(t) - (2n+1) t P_{n}(t) + n P_{n-1}(t) = 0.\]

En particulier :

\[\begin{split}P_{0}(t) = 1, \\ P_{1}(t) = t, \\ P_{2}(t) = \frac{3t^2 - 1}{2}, \\ P_{3}(t) = \frac{5t^3 - 3t}{2}.\end{split}\]

On les normalise sous la forme :

\[{\gamma}_{j}(s) = \sqrt{\frac{2j+1}{l}} P_{j} \left(\frac{2s}{l} - 1\right),\]

où :

• \(s\) est l’abscisse curviligne de \({\Gamma}_{0}\),

• \(l\) est la longueur du fond de fissure \({\Gamma}_{0}\).

Fig. 302 Figure 2.3-b : Polynômes de Legendre

On se limite à \(\mathrm{Deg}_{\max} = 7\) comme degré maximal.

Les fonctions de forme \(\varphi_{k}(s)\) sont associées à la discrétisation de \({\Gamma}_{0}\).

Fig. 303 Figure 2.3-c : Fonctions de forme du fond de fissure (éléments linéaires)

Rappelons que l’on est amené à discrétiser \(G(s)\) le long du fond de fissure :

\[G(s) = \sum_{j=1}^{N} G_{j} p_{j}(s).\]

et à se donner une base de fonctions tests \({(\bar{\theta}^{i}(s))}_{1\leq i\leq N}\) pour la trace du champ \(\theta\) sur le fond de fissure \({\Gamma}_{0}\).

Il existe donc plusieurs choix possibles de discrétisations, résumés dans le tableau ci-dessous :

Tableau 76 Tableau 2.3-1 : Choix de la discrétisation

	Polynômes de LEGENDRE	Fonctions de forme
\(G(s)\)	\(G(s)=\sum_{j=0}^{\mathrm{NDEG}}{G}_{j}{\gamma}_{j}(s)\)	\(\sum_{j=1}^{\mathrm{NNO}}{G}_{j}{\varphi}_{j}(s)\)
\(\bar{\theta}(s)\)	\({({\gamma}_{i}(s))}_{1\leq i\leq \mathrm{NDEG}}\)	\({({\varphi}_{i}(s))}_{1\leq i\leq \mathrm{NNO}}\)

Premier cas : Discrétisation LEGENDRE

\(G(s)\) est décomposé suivant les polynômes de Legendre :

\[G(s) = \sum_{j=0}^{\mathrm{NDEG}} G_{j} \gamma_{j}(s).\]

La base des fonctions tests pour \(\theta (s)\) est définie à partir des polynômes de Legendre :

\[\tilde{\Theta} = \lbrace {\theta}^{i} \in \Theta \mid \bar{\theta}^{i}(s) \cdot m(s) = {\varphi}_{i}(s), \quad i = 1, \dots, \mathrm{NDEG} \rbrace.\]

Les \(\mathrm{NDEG}\) composantes \(G_{j}\) sont déterminées en résolvant le système linéaire à \(\mathrm{NDEG}\) équations :

\[\begin{split}\begin{cases} \sum_{j=0}^{\mathrm{NDEG}} a_{ij} G_{j} = b_{i}, \quad i = 1, \dots, \mathrm{NDEG}, \\ \text{avec} \quad \begin{cases} a_{ij} = \int_{{\Gamma}_{0}} \gamma_{i}(s) \gamma_{j}(s) \, ds, \\ b_{i} = G(\theta^{i}). \end{cases} \end{cases}\end{split}\]

Les polynômes de Legendre formant une base orthonormée sur \({\Gamma}_{0}\), on a \(a_{ij} = \delta_{ij}\) et le système linéaire se résume donc à :

\[G_{j} = G(\theta^{j}),\]

et donc :

\[G(s) = \sum_{j=0}^{\mathrm{NDEG}} G(\theta^{j}) \gamma_{j}(s).\]

Deuxième cas : Discrétisation LINEAIRE

\(G(s)\) est défini par les fonctions de forme des \(\mathrm{NNO}\) nœuds du fond de fissure :

\[G(s) = \sum_{j=0}^{\mathrm{NNO}} G_{j} \varphi_{j}(s).\]

La base des fonctions tests pour \(\theta (s)\) est définie à partir des fonctions de forme des \(\mathrm{NNO}\) nœuds du fond de fissure :

\[\tilde{\Theta} = \lbrace {\theta}^{i} \in \Theta \mid \bar{\theta}^{i}(s) \cdot m(s) = {\varphi}_{i}(s), \quad i = 1, \dots, \mathrm{NNO} \rbrace,\]

où \({\varphi}_{i}(s)\) est la fonction de forme du nœud \(i\) du fond de fissure.

Système à résoudre :

\[\begin{split}\begin{cases} \sum_{j=0}^{\mathrm{NNO}} a_{ij} G_{j} = G(\theta^{i}), \quad (i=1,\dots,\mathrm{NNO}), \\ \text{avec} \quad a_{ij} = \int_{{\Gamma}_{0}} \varphi_{i}(s) \varphi_{j}(s) \, ds. \end{cases}\end{split}\]

Éléments linéaires : Si on a des éléments linéaires, les fonctions de forme sont définies comme suit :

\[\varphi_{1}(x) = \frac{1}{2} (1 - x), \quad \varphi_{2}(x) = \frac{1}{2} (1 + x).\]

Les coefficients de la matrice \(A_{ij}\) sont donnés par :

\[a_{i(i-j)} = a_{i(i+j)} = 0, \quad \text{si } j \geq 2.\]

\[a_{i(i-1)} = \int_{{\Gamma}_{0}} \varphi_{i}(s) \varphi_{i-1}(s) \, ds = \int_{s_{i-1}}^{s_{i}} \varphi_{i}(s) \varphi_{i-1}(s) \, ds.\]

\[a_{i(i-1)} = \frac{s_{i} - s_{i-1}}{2} \int_{-1}^{+1} \varphi_{i}(x) \varphi_{i-1}(x) \, dx = \frac{s_{i} - s_{i-1}}{2} \int_{-1}^{+1} \frac{1}{4} (1 - x^2) \, dx = \frac{1}{6} (s_{i} - s_{i-1}).\]

\[a_{ii} = \frac{s_{i} - s_{i-1}}{2} \int_{-1}^{+1} \varphi_{2}^{2}(x) \, dx + \frac{s_{i+1} - s_{i}}{2} \int_{-1}^{+1} \varphi_{1}^{2}(x) \, dx.\]

\[a_{ii} = \frac{s_{i} - s_{i-1}}{2} \int_{-1}^{+1} \frac{1}{4} (1 + x)^2 \, dx + \frac{s_{i+1} - s_{i}}{2} \int_{-1}^{+1} \frac{1}{4} (1 - x)^2 \, dx.\]

\[a_{ii} = \frac{1}{3} \left[ (s_{i+1} - s_{i}) + (s_{i} - s_{i-1}) \right].\]

Dans le cas d’un fond ouvert, la matrice \(A_{ij}\) est :

\[\begin{split}\frac{1}{6} \begin{bmatrix} 2(s_2 - s_1) & (s_2 - s_1) & 0 & 0 & \dots & 0 & 0 \\ (s_2 - s_1) & 2(s_3 - s_1) & (s_3 - s_2) & 0 & \dots & 0 & 0 \\ 0 & (s_3 - s_2) & 2(s_4 - s_2) & (s_4 - s_3) & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \dots & 2(s_n - s_{n-2}) & (s_n - s_{n-1}) \\ 0 & 0 & 0 & 0 & \dots & (s_n - s_{n-1}) & 2(s_n - s_{n-1}) \end{bmatrix}\end{split}\]

Dans le cas d’un fond fermé, la condition \(b_1 = b_{\mathit{NNO}}\) implique que \(G_1 = G_{\mathit{NNO}}\). Il faut donc rajouter des termes à la matrice du système linéaire \(A_{ij}\) sur la première et la dernière ligne :

\[\begin{split}\frac{1}{6} \begin{bmatrix} 2(s_2 - s_1) + 2(s_n - s_{n-1}) & (s_2 - s_1) & 0 & \dots & (s_n - s_{n-1}) & 0 \\ (s_2 - s_1) & 2(s_3 - s_1) & (s_3 - s_2) & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 2(s_n - s_{n-2}) & (s_n - s_{n-1}) \\ 0 & (s_2 - s_1) & 0 & \dots & (s_n - s_{n-1}) & 2(s_2 - s_1) + 2(s_n - s_{n-1}) \end{bmatrix}\end{split}\]

Éléments quadratiques: Si on a des éléments quadratiques, les fonctions de forme sont définies comme suit :

\[\varphi_{1}(x) = \frac{1}{2} x (x - 1)\]

\[\varphi_{2}(x) = \frac{1}{2} x (x + 1)\]

\[\varphi_{3}(x) = (1 - x)(1 + x)\]

Il faut distinguer nœud sommet et nœud milieu.

1. Nœud sommet (\(i\) = nœud sommet)

\[a_{i(i-j)} = a_{i(i+j)} = 0, \quad \text{si } j \geq 3\]

\[a_{i(i-2)} = \int_{{\Gamma}_{0}} \varphi_{i}(s) \varphi_{i-2}(s) \, ds\]

\[a_{i(i-2)} = \int_{s_{i-2}}^{s_{i}} \varphi_{i}(s) \varphi_{i-2}(s) \, ds\]

\[a_{i(i-2)} = \frac{s_{i} - s_{i-2}}{2} \int_{-1}^{+1} \varphi_{1}(x) \varphi_{2}(x) \, dx\]

\[a_{i(i-2)} = \frac{s_{i} - s_{i-2}}{2} \int_{-1}^{+1} \frac{1}{4} x^2 (x^2 - 1) \, dx\]

\[a_{i(i-2)} = -\frac{s_{i} - s_{i-2}}{30}\]

\[a_{i(i-1)} = \int_{{\Gamma}_{0}} \varphi_{i}(s) \varphi_{i-1}(s) \, ds\]

\[a_{i(i-1)} = \int_{s_{i-2}}^{s_{i}} \varphi_{i}(s) \varphi_{i-1}(s) \, ds\]

\[a_{i(i-1)} = \frac{s_{i} - s_{i-2}}{2} \int_{-1}^{+1} \varphi_{2}(x) \varphi_{3}(x) \, dx\]

\[a_{i(i-1)} = \frac{s_{i} - s_{i-2}}{2} \int_{-1}^{+1} \frac{1}{2} x (x+1)^2 (1-x) \, dx\]

\[a_{i(i-1)} = \frac{s_{i} - s_{i-2}}{15}\]

\[a_{ii} = \frac{s_{i} - s_{i-2}}{2} \int_{-1}^{+1} \varphi_{2}^{2}(x) \, dx + \frac{s_{i+2} - s_{i}}{2} \int_{-1}^{+1} \varphi_{1}^{2}(x) \, dx\]

\[a_{ii} = \frac{2}{15} (s_{i+2} - s_{i-2})\]

Fig. 304 Figure 2.3-f : Nœud sommet

2. Nœud milieu (\(i\) = nœud milieu)

\[a_{i(i-j)} = a_{i(i+j)} = 0, \quad \text{si } j \geq 2\]

\[a_{i(i-1)} = \frac{s_{i+1} - s_{i-1}}{2} \int_{-1}^{+1} \varphi_{3}(x) \varphi_{1}(x) \, dx\]

\[a_{i(i-1)} = \frac{s_{i+1} - s_{i-1}}{2} \int_{-1}^{+1} \frac{1}{2} x(x-1)(1-x)(1+x) \, dx\]

\[a_{i(i-1)} = \frac{s_{i+1} - s_{i-1}}{15}\]

\[a_{ii} = \int_{s_{i-1}}^{s_{i+1}} \varphi_{i}^{2}(s) \, ds\]

\[a_{ii} = \frac{s_{i+1} - s_{i}}{2} \int_{-1}^{+1} (1-x)^2 (1+x)^2 \, dx\]

\[a_{ii} = \frac{8}{15} (s_{i+1} - s_{i-1})\]

Fig. 305 Figure 2.3-g : Nœud milieu

La matrice \(A_{ij}\) s’écrit :

\[\begin{split}\frac{1}{30} \begin{bmatrix} 4(s_3 - s_1) & 2(s_3 - s_1) & -(s_3 - s_1) & 0 & 0 & \dots & 0 \\ 2(s_3 - s_1) & 16(s_3 - s_1) & 2(s_3 - s_1) & 0 & 0 & \dots & 0 \\ -(s_3 - s_1) & 2(s_3 - s_1) & 4(s_5 - s_1) & 2(s_5 - s_3) & -(s_5 - s_3) & \dots & 0 \\ 0 & 0 & 2(s_5 - s_3) & 16(s_5 - s_3) & 2(s_5 - s_3) & \dots & 0 \\ 0 & 0 & -(s_5 - s_3) & 2(s_5 - s_3) & 4(s_7 - s_3) & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & 0 & 0 & 0 & \dots & -(s_n - s_{n-2}) \\ 0 & 0 & 0 & 0 & 0 & \dots & 2(s_n - s_{n-2}) \\ 0 & 0 & 0 & 0 & 0 & \dots & 4(s_n - s_{n-2}) \end{bmatrix}\end{split}\]

Cas d’un fond fermé Dans ce cas, on rajoute des termes à la matrice pour respecter la condition \(G_1 = G_{\mathit{NNO}}\). La matrice \(A_{ij}\) devient :

\[\begin{split}\frac{1}{30} \begin{bmatrix} 4(s_3 - s_1) + 4(s_n - s_{n-2}) & 2(s_3 - s_1) & -(s_3 - s_1) & 0 & \dots & -(s_n - s_{n-1}) & 2(s_n - s_{n-2}) & 0 \\ 2(s_3 - s_1) & 16(s_3 - s_1) & 2(s_3 - s_1) & 0 & 0 & \dots & 0 & 0 \\ -(s_3 - s_1) & 2(s_3 - s_1) & 4(s_5 - s_1) & 2(s_5 - s_3) & -(s_5 - s_3) & 0 & \dots & 0 \\ 0 & 0 & 2(s_5 - s_3) & 16(s_5 - s_3) & 2(s_5 - s_3) & 0 & \dots & 0 \\ 0 & 0 & -(s_5 - s_3) & 2(s_5 - s_3) & 4(s_7 - s_3) & 0 & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & 0 & 0 & 0 & 0 & \dots & -(s_n - s_{n-2}) \\ 0 & 0 & 0 & 0 & 0 & 0 & \dots & 2(s_n - s_{n-2}) \\ 0 & 2(s_3 - s_1) & -(s_3 - s_1) & 0 & 0 & 0 & \dots & 4(s_3 - s_1) + 4(s_n - s_{n-2}) \end{bmatrix}\end{split}\]

Correction des valeurs aux extrémités du fond de fissure pour la discrétisation LINEAIRE

Lorsque le fond de fissure débouche, un traitement supplémentaire particulier est réalisé.

Il a été observé que, pour la discrétisation linéaire, les valeurs numériques des \(G({\theta}^{i})\) pour \(i=1\) et \(i=n\) sont particulièrement «fausses». Cela est probablement dû à l’ordre de la singularité qui n’est plus en \(-1/2\) à proximité des surfaces libres.

En effet, l’ordre de la singularité aux bords dépend de \(\nu\) (et de l’angle entre le fond de fissure et la surface libre, ainsi que de l’angle entre la normale à la surface de la fissure et la normale à la surface libre). Même pour une fissure orthogonale à la surface libre, la singularité n’est pas en \(-1/2\) si \(\nu \neq 0\).

Correction aux points de bord Afin d’améliorer les résultats aux points du bord, on modifie les valeurs de \(G({\theta}^{i})\) pour \(i=1\) et \(i=n\) dans le cas de la discrétisation linéaire. On modifie également les coefficients \(K_{j}({\theta}^{i})\), pour \(j=I, II, III\).

Ce traitement est appliqué même en présence de conditions de symétrie sur les bords du fond, bien que dans ce cas, l’ordre de la singularité soit effectivement en \(-1/2\) et aucun traitement ne soit théoriquement nécessaire.

Cas 1 : Discrétisation linéaire et maillage linéaire On fait l’hypothèse que :

\[G(s_1) = G(s_2) = G(s_3) = G_a, \quad G(s_n) = G(s_{n-1}) = G(s_{n-2}) = G_b.\]

Alors, en multipliant la matrice \(A\) par le vecteur \(G(s)\), on obtient :

\[G({\theta}^{1}) = \frac{1}{2} (s_2 - s_1) G_a\]

\[G({\theta}^{2}) = \frac{1}{2} (s_3 - s_1) G_a\]

\[\vdots\]

\[G({\theta}^{n-1}) = \frac{1}{2} (s_n - s_{n-2}) G_b\]

\[G({\theta}^{n}) = \frac{1}{2} (s_n - s_{n-1}) G_b\]

Ce qui conduit aux relations :

\[G({\theta}^{1}) = G({\theta}^{2}) \frac{(s_2 - s_1)}{(s_3 - s_1)}\]

\[G({\theta}^{n}) = G({\theta}^{n-1}) \frac{(s_n - s_{n-1})}{(s_n - s_{n-2})}\]

Cas 2 : Discrétisation linéaire et maillage quadratique On fait l’hypothèse que :

\[G(s_1) = G(s_2) = G(s_3) = G_a, \quad G(s_n) = G(s_{n-1}) = G(s_{n-2}) = G_b.\]

Alors, en multipliant la matrice \(A\) par le vecteur \(G(s)\), on obtient :

\[G({\theta}^{1}) = \frac{1}{6} (s_3 - s_1) G_a\]

\[G({\theta}^{2}) = \frac{4}{6} (s_3 - s_1) G_a\]

\[\vdots\]

\[G({\theta}^{n-1}) = \frac{4}{6} (s_n - s_{n-2}) G_b\]

\[G({\theta}^{n}) = \frac{4}{6} (s_n - s_{n-2}) G_b\]

Ce qui conduit aux relations :

\[G({\theta}^{1}) = \frac{1}{4} G({\theta}^{2})\]

\[G({\theta}^{n}) = \frac{1}{4} G({\theta}^{n-1})\]

Implantation de G en thermo-élasticité linéaire

Types d’éléments et de chargements

Il est possible de calculer en thermo-élasticité linéaire le taux de restitution d’énergie :

• En 2D : \(G(\theta )\)

• En 3D (valeur locale) : \(G(s)\)

Ces calculs sont réalisés à l’aide de la commande CALC_G [U4.82.03].

Le champ de propagation virtuel de la fissure \(\theta\) est calculé automatiquement dans la commande CALC_G, à partir des informations fournies par l’utilisateur.

Ces calculs sont valables pour les modélisations suivantes :

• D_PLAN (fissure maillée ou non maillée)

• C_PLAN (fissure maillée ou non maillée)

• AXIS (fissure maillée ou non maillée)

• 3D (fissure maillée ou non maillée)

Les calculs sont valables pour les chargements thermo-mécaniques suivants :

• \(f\) : champ de forces volumiques appliquées sur \(\Omega\) (charges mécaniques du type PESANTEUR, ROTATION, FORCE_INTERNE).

• \(g\) : champ de forces surfaciques appliquées sur une partie \(S\) de \(\delta \Omega\) (y compris sur les lèvres de la fissure : PRES_REP, FORCE_FACE).

• \(T\) : champ de température (transmis via le champ matériau : AFFE_MATERIAU / AFFE_VARC).

• \({\varepsilon}^{0}\) : champ de pré-déformation (PRE_EPSI).

Ces chargements peuvent dépendre du temps et de l’espace. On ne tient pas compte du terme dû aux déplacements imposés sur \({S}_{d}\). Il ne faut donc pas imposer de conditions de Dirichlet sur les lèvres de la fissure.

Les propriétés du matériau \((E, \nu, \alpha )\) peuvent :

• Dépendre de la température \(T\).

• Varier dans l’espace, tout en restant constantes par éléments.

Un champ de contrainte initiale \({\sigma}^{0}\) peut également être pris en compte.

Environnement nécessaire

Fissure maillée:

Pour le calcul du taux de restitution (option G), il faut caractériser préalablement la fissure via le concept de type fond_fiss (en 3D: liste ordonnée des noeuds du fond, mailles de la lèvre supérieure et la lèvre inférieure; en 2D: nœud de fond de fissure et normale à la fissure), avec la commande DEFI_FOND_FISS [U4.82.01].

Lorsque que la fissure est disposée le long d’un axe de symétrie, on peut se contenter de ne représenter que la moitié du modèle, et préciser la symétrie du chargement par le mot-clé SYME dans DEFI_FOND_FISS. Par défaut, on suppose qu’il n’y a pas de symétrie. Si le modèle présente une symétrie, CALC_G la détecte automatiquement via le concept fond_fiss.

Fissure non maillée (méthode X-FEM):

En 2D comme en 3D, la fissure doit être définie, pour le calcul mécanique et pour le post-traitement, à l’aide de la commande DEFI_FISS_XFEM. Le mot-clé FISSURE doit être renseigné dans CALC_G_XFEM.

Si la fissure n’est pas maillée, il n’est pas possible de prendre en compte les éventuelles symétries du modèle par rapport aux lèvres de la fissure.

Le calcul de \(G\) se fait en post-traitement, uniquement à partir du champ de déplacement solution du calcul sur le modèle considéré. En particulier, la densité d’énergie libre et les contraintes sont calculées à partir du champ de déplacement et des caractéristiques du matériau.

Calculs des différents termes du taux de restitution d’énergie

L’expression complète de \(G(\theta )\) est donnée au [§1.3]. Nous allons détailler chaque terme. Le champ \(\theta\) est nul en dehors d’un disque de rayon \({R}_{\sup}(s)\) défini dans le chapitre 3 [Figure 3.3-a]. Remarquons que comme tous les termes font intervenir \(\theta\) ou son gradient, les termes élémentaires sont nuls en dehors de ce disque de rayon \({R}_{\sup}(s)\) . Dans la commande CALC_G, il n’est ainsi pas nécessaire de préciser les chargements qui ne s’appliquent pas dans cette zone.

Terme Classique Élémentaire

Le terme classique élémentaire (TCLA) est défini par :

\[\text{TCLA} = {\sigma}_{ij} {u}_{i,p} {\theta}_{p,j} - \psi (\varepsilon (u),T) {\theta}_{k,k}\]

La densité d’énergie élastique \(\Psi (\varepsilon (u),T)\) s’exprime en thermo-élasticité linéaire comme suit :

1. En 3D et en AXIS

\[\Psi (\varepsilon (u),T) = \frac{1}{2} \lambda ({\varepsilon}_{ii})^2 + \mu {\varepsilon}_{ij} {\varepsilon}_{ij} - {\psi}_{\mathrm{th}}\]

2. En déformations planes

\[\Psi (\varepsilon (u),T) = \frac{(1-\nu)E}{2(1+\nu)(1-2\nu)} ({\varepsilon}_{xx}^{2} + {\varepsilon}_{yy}^{2}) + \frac{\nu E}{(1+\nu)(1-2\nu)} {\varepsilon}_{xx} {\varepsilon}_{yy} + \frac{E}{(1+\nu)} {\varepsilon}_{xy}^{2} - {\psi}_{\mathrm{th}}\]

3. En contraintes planes

\[\Psi (\varepsilon (u),T) = \frac{E}{2(1-\nu^2)} ({\varepsilon}_{xx}^{2} + {\varepsilon}_{yy}^{2}) + \frac{\nu E}{(1-\nu^2)} {\varepsilon}_{xx} {\varepsilon}_{yy} + \frac{E}{(1+\nu)} {\varepsilon}_{xy}^{2} - {\psi}_{\mathrm{th}}\]

où la composante thermique est donnée par :

\[{\Psi}_{\mathrm{th}} = 3K\alpha (T - T_{\mathrm{réf}}) {\varepsilon}_{ii} - \frac{9}{2} K {\alpha}^{2} (T - T_{\mathrm{réf}})^2\]

Les paramètres sont définis comme suit :

\[3K = \frac{E}{1-2\nu}, \quad \lambda = \frac{\nu E}{(1+\nu)(1-2\nu)}, \quad 2\mu = \frac{E}{1+\nu}\]

• \(E\) : Module de Young

• \(\nu\) : Coefficient de Poisson

• \(\lambda, \mu\) : Coefficients de Lamé

• \(\alpha\) : Dilatation thermique

La densité d’énergie élastique \(\psi (\varepsilon (u),T)\) peut s’écrire sous une forme générale :

\[\Psi (\varepsilon (u),T) = \frac{1}{2} K (\varepsilon_{kk} - 3\alpha (T - T_{\mathrm{réf}}))^2 + \frac{2\mu}{3} \varepsilon_{\mathrm{eq}}^2\]

avec :

\[\varepsilon_{\mathrm{eq}}^2 = \frac{3}{2} \varepsilon_{ij}^{D} \varepsilon_{ij}^{D}, \quad \varepsilon_{ij}^{D} = \varepsilon_{ij} - \frac{1}{3} \varepsilon_{kk} \delta_{ij}\]

ou encore :

\[\varepsilon_{\mathrm{eq}}^2 = \frac{1}{2} (3\varepsilon_{ij} \varepsilon_{ij} - \varepsilon_{kk}^2)\]

Ainsi, la densité d’énergie élastique peut aussi s’écrire sous la forme suivante :

\[\Psi (\varepsilon (u),T) = \frac{1}{2} K \varepsilon_{kk}^2 - 3K\alpha (T - T_{\mathrm{réf}}) \varepsilon_{kk} + \frac{9}{2} K \alpha^2 (T - T_{\mathrm{réf}})^2 + \mu \varepsilon_{ij} \varepsilon_{ij} - \frac{\mu}{3} \varepsilon_{kk}^2\]

ou encore :

\[\Psi (\varepsilon (u),T) = \frac{\lambda}{2} \varepsilon_{kk}^2 + \mu \varepsilon_{ij} \varepsilon_{ij} - 3K\alpha (T - T_{\mathrm{réf}}) \varepsilon_{kk} + \frac{9}{2} K \alpha^2 (T - T_{\mathrm{réf}})^2\]

Terme de Force Volumique

Le terme de force volumique (TFOR) est défini par :

\[\mathrm{TFOR} = f_{i} u_{i} \theta_{k,k} + f_{i,k} \theta_{k} u_{i}\]

Terme de Force Surfacique

Le terme de force surfacique (TSUR) s’écrit :

\[\mathrm{TSUR} = g_{i,k} \theta_{k} u_{i} + g_{i} u_{i} \left(\theta_{k,k} - \frac{\partial \theta}{\partial n_{k}} n_{k} \right)\]

Remarque : Dans ce terme surfacique, des dérivations normales à la surface apparaissent. Elles n’ont pas de sens pour les éléments de peau utilisés. On utilise donc la géométrie différentielle et les dérivées contravariantes pour mieux appréhender cet intégrande sur la surface de calcul.

Terme Thermique

Le terme thermique (THER) s’écrit :

\[\mathrm{THER} = -\frac{\partial \psi}{\partial T} T_{,k} \theta_{k}\]

avec :

\[\frac{\partial \Psi}{\partial T} (\varepsilon (u),T) = \left[ \frac{1}{2} \frac{\mathrm{dK}(T)}{\mathrm{dT}} (\varepsilon_{kk} - 3\alpha (T - T_{\mathrm{réf}})) - 3K \left( \alpha + \frac{d\alpha (T)}{\mathrm{dT}} (T - T_{\mathrm{réf}}) \right) \right] (\varepsilon_{kk} - 3\alpha (T - T_{\mathrm{réf}}))\]

Terme de Pré-déformations et Contraintes Initiales

Le terme de pré-déformations et contraintes initiales (TINI) s’écrit :

\[\mathrm{TINI} = \left[ \sigma_{ij} \varepsilon_{ij,k}^{0} - \Lambda_{\mathrm{ijpq}}^{-1} \sigma_{\mathrm{pq}} \sigma_{ij,k}^{0} \right] \theta_{k}\]

On peut remarquer que si \(\sigma = 0\), alors \(\mathrm{TINI} = 0\).

Remarque 1 : Comme le reste de cette documentation, la présence de contraintes initiales n’est valable que dans le cas de l’élasticité linéaire isotrope. Pour être activée, il faut être dans une **relation de comportement incrémentale* (COMPORTEMENT) avec la relation ELAS. Il n’est pas possible de fournir un champ de contraintes initiales dans un autre cas.*

Remarque 2 : Dans tous les cas, le champ de contraintes initiales doit être **auto-équilibré* en l’absence de fissure, avec les seules conditions limites. L’utilisateur pourra (et devrait) vérifier que son champ de contraintes initiales est licite en l’appliquant dans ETAT_INIT de l’opérateur STAT_NON_LINE, avec un comportement élastique incrémental (COMPORTEMENT, RELATION = “ELAS”). Le résultat mécanique doit alors donner le même champ de contrainte sans déformations supplémentaires.*

Remarque 3 : Le champ de contraintes initiales devra avoir été utilisé préalablement dans le calcul thermo-mécanique. La démarche de calcul du taux de restitution d’énergie en présence de contraintes initiales est détaillée ci-dessous. On pourra également se référer au cas test dédié.

Introduction du champ \(\theta\)

Conditions à remplir

Le champ \(\theta\) est un champ de vecteurs, défini sur le solide fissuré, qui représente la transformation du domaine lors d’une propagation de fissure au sens du [§1]. La transformation ne doit modifier que la position du fond de fissure et pas le bord du domaine \(\delta \Omega\) , c’est-à-dire: \(\theta .n=0\) sur (\(n\) normale à \(\delta \Omega\) ). De plus, le champ \(\theta\) doit être régulier sur \(\Omega\) [bib4].

En raison de la singularité du champ de déplacement, il est intéressant du point de vue numérique d’utiliser des champs \(\theta\) constants dans un voisinage de \({\Gamma}_{0}\) , annulant ainsi dans ce voisinage les termes singuliers \(\psi {\theta}_{k,k}-{\sigma}_{ij}{u}_{i,p}{\theta}_{p,k}\) dans \(G(\theta )\) .

Choix du Champ \(\theta\) en Dimension 3

Méthode de Construction

Le champ \(\theta\) doit être construit de manière à vérifier les conditions suivantes :

\[\begin{split}\begin{cases} \theta_n = \theta \cdot n = 0, & \text{sur le bord du domaine } \delta \Omega \text{ (où } n \text{ est la normale à } \delta \Omega \text{)} \\ \bar{\theta} = \bar{\theta}_0, & \text{donné sur le fond de fissure } \Gamma_0 \end{cases}\end{split}\]

où \(\bar{\theta}\) représente la trace de \(\theta\) sur \(\Gamma_0\).

On définit deux surfaces toriques \(T\) et \(S\) (cylindres déformés) entourant le fond de fissure \(\Gamma_0\). Pour vérifier la première condition ci-dessus, il est crucial que la couronne \(S\) soit strictement incluse dans \(\Omega\).

Fig. 306 Figure 3.2.1-a : Construction du champ \(\theta\) en 3D (vue d’ensemble)

On note :

\(R_{\inf}(s)\) : rayon variable de \(T\) et \(R_{\sup}(s)\) : rayon variable de \(S\)

Fig. 307 Figure 3.2.1-b : Construction du champ \(\theta\) en 3D (plan de coupe)

En tout point de \(\Gamma_0\), repéré par son abscisse curviligne \(s\), on définit un plan normal \(P\), dans lequel le champ \(\theta\) est introduit selon les règles suivantes :

• À l’intérieur du tore \(T\) :

  \[\theta_n(r(s)) = \bar{\theta}_0(s), \quad \text{pour } 0 \leq r(s) \leq R_{\inf}(s)\]

• À l’extérieur du cylindre \(S\) :

  \[\theta_n(r(s)) = 0, \quad \text{pour } r(s) \geq R_{\sup}(s)\]

• Dans la couronne intermédiaire \(S \setminus T\) :

  \[\theta_n \text{ varie linéairement en fonction de } r(s)\]

• Condition de continuité :

  \[\theta_n \text{ est continu dans } S(R_{\sup}(s))\]

Cette approche est géométrique : Elle repose sur la définition de deux rayons \(R_{\inf}(s)\) et \(R_{\sup}(s)\), et l’utilisation de calculs de distance entre un point courant et le fond de fissure pour déterminer la valeur de \(\theta\) en ce point.

Algorithmes de calcul

Le module \(|\theta|\) du champ \(\theta\) est donné sur tout le fond de fissure \(\Gamma_0\). Il peut être défini par l’utilisateur ou calculé automatiquement via la méthode \(\theta\).

Le problème consiste à déterminer le champ \(\theta\) (norme et direction) en tout point du domaine \(\Omega\).

On note :

• \(|\theta|_i\) : module du champ \(\theta\) au nœud \(i\) de \(\Gamma_0\)

• \(R_{\inf}^i\) : rayon intérieur pour le nœud \(i\)

• \(R_{\sup}^i\) : rayon extérieur pour le nœud \(i\)

La procédure de calcul du champ \(\theta\) en tout point de \(\Omega\) est la suivante:

Calcul du champ \(\theta\) en chaque point de \(\Gamma_0\)

Le module \(|\theta|_i\) étant donné, il faut déterminer la direction de \(\theta\).Les conditions géométriques à respecter sont : \(\theta\) doit être localement dans le plan tangent aux lèvres de la fissure, \(\theta\) doit être normal à l’arête à laquelle il appartient.

En cas général (fond de fissure non plan), \(\theta\) étant calculé aux nœuds, sa direction est moyennée sur les deux arêtes de \(\Gamma_0\) ayant ce nœud en commun.

Fig. 308 Figure 3.2.2-a : Construction du champ \(\theta\) en 3D (normales)

Soient \(F_1\) et \(F_2\) deux faces appartenant aux lèvres de la fissure et comprenant les arêtes successives \(T_1\) et \(T_2\) du fond de fissure \(\Gamma_0\).

On commence par calculer :

• La normale \(n_1\) à l’arête \(T_1\) dans le plan de la face \(F_1\)

• La normale \(n_2\) à l’arête \(T_2\) dans le plan de la face \(F_2\)

Les normales étant unitaires, on en déduit :

\[n = \frac{n_1 + n_2}{2}\]

Le champ \(\theta\) en un point \(M \in \Gamma_0\) est alors donné par :

\[\theta (M) = |\theta (M)| \cdot n\]

On considère que les faces \(F_i\) sont planes :

• Si \(F_i\) est un triangle, le plan de la face \(F_i\) est directement défini.

• Si \(F_i\) est un quadrangle, on le divise en deux triangles \(F_{i1}\) et \(F_{i2}\). Il faut alors calculer les équations des deux plans contenant ces faces et effectuer deux calculs de normales par arête \(T_i\).

Ce calcul nécessite de connaître les faces appartenant aux lèvres de la fissure et comprenant une arête de \(\Gamma_0\). L’utilisateur doit fournir tous les éléments surfaciques appartenant aux lèvres de la fissure. Ces faces figurent dans un ou plusieurs groupes de mailles et sont décrites dans les connectivités des éléments de surfaces. L’algorithme trie ces faces pour ne conserver que celles ayant deux sommets sur \(\Gamma_0\).

Les étapes de l’algorithme sont :

1. Pour chaque nœud de \(\Gamma_0\), extraction des mailles appartenant aux lèvres de la fissure.

2. Filtrage des mailles pour ne garder que celles ayant deux nœuds sur \(\Gamma_0\).

3. Récupération du type de la face (triangle ou quadrangle) et calcul de l’équation du ou des plans tangents.

4. Pour chaque arête de \(\Gamma_0\) de sommets \(N_i\) et \(N_{i+1}\), calcul des normales \(n_{i,1}\), \(n_{i+1,1}\), \(n_{i,2}\) et \(n_{i+1,2}\).

Le champ \(\theta\) est ensuite calculé selon l’algorithme suivant :

\[\begin{split}\forall N_i \in \Gamma_0 : \begin{cases} n_i = \frac{1}{2} (n_{i,1} + n_{i,2}) \\ \theta (N_i) = |\theta|_i \cdot n_i \end{cases}\end{split}\]

Fin de la boucle sur les sommets \(N_i\) de \(\Gamma_0\).

Fig. 309 Figure 3.2.2-b : Notations des normales au fond de fissure

Calcul de la projection de chaque point de \(\Omega\) sur \(\Gamma_0\)

Pour chaque nœud \(M\) :

• Récupération des coordonnées de \(M\)

• Boucle sur les nœuds \(N_k\) de \(\Gamma_0\), avec \(k=1,...,\mathrm{NNO}-1\) :

  • Récupération des coordonnées de \(N_k\) et \(N_{k+1}\)

  • Calcul de :

    \[s_k = \frac{N_k N_{k+1} \cdot N_k M}{\| N_k N_{k+1} \|}\]

  • Ajustement de \(s_k\) :

    \[\begin{split}\begin{cases} s_k < 0 \Rightarrow s_k = 0 \\ s_k > 1 \Rightarrow s_k = 1 \end{cases}\end{split}\]

  • Calcul des coordonnées de la projection :

    \[\bar{M}_k = N_k + s_k (N_k N_{k+1})\]

  • Calcul de la distance :

    \[d_k = d(M, \bar{M}_k)\]

• Fin de la boucle

• Détermination de \(i\) tel que :

  \[d_i = \min_k (d_k)\]

• Connaissant \(i\), on récupère \(N_i\), \(N_{i+1}\), \(s_i\) et la projection \(\bar{M}\) de \(M\) sur \(\Gamma_0\), telle que :

  \[N_i \bar{M} = s_i N_i N_{i+1}\]

Calcul du champ \(\theta\) en chaque point de \(\Omega\)

Boucle sur les nœuds \(M\) :

1. Calcul de la projection \(\bar{M}\) de \(M\) sur \(\Gamma_0\) (cf. ci-dessus, donne \(N_i, N_{i+1}, s_i\)).

2. Calcul de la distance :

   \[d = d(M, \bar{M})\]

3. Calcul de \(\theta(\bar{M})\) par interpolation linéaire :

   \[\theta (\bar{M}) = (1 - s_i) \theta (N_i) + s_i \theta (N_{i+1})\]

4. Calcul de \(R_{\inf}(s_i)\) et \(R_{\sup}(s_i)\) par interpolation linéaire :

   \[ \begin{align}\begin{aligned}R_{\inf}(s_i) = (1 - s_i) R_{\inf}^i + s_i R_{\inf}^{i+1}\\R_{\sup}(s_i) = (1 - s_i) R_{\sup}^i + s_i R_{\sup}^{i+1}\end{aligned}\end{align} \]

5. Attribution de \(\theta(M)\) selon la valeur de \(d\) :

   \[\begin{split}\begin{cases} d > R_{\sup}(s_i) \Rightarrow \theta(M) = 0 \\ d < R_{\inf}(s_i) \Rightarrow \theta(M) = \theta(\bar{M}) \\ R_{\inf}(s_i) \leq d \leq R_{\sup}(s_i) \Rightarrow \alpha = \frac{d - R_{\inf}(s_i)}{R_{\sup}(s_i) - R_{\inf}(s_i)}, \quad \theta(M) = (1 - \alpha) \theta(\bar{M}) \end{cases}\end{split}\]

• Fin de la boucle

Fig. 310 Figure 3.2.2-d : Calcul du champ \(\theta\) en 3D

Choix du champ \(\theta\) en dimension 2

Il s’agit d’un cas particulier de la dimension 3. Le fond de fissure \(\Gamma_0\) se réduit à un point. L’utilisateur choisit :

• les rayons \(R_{\inf}\) et \(R_{\sup}\)

• le module \(|\theta_0|\) en fond de fissure

Le champ \(\theta\) est construit de manière à respecter les conditions suivantes :

\[\begin{split}\begin{cases} \theta (r) = 0, \quad \text{si } r \geq R_{\sup} \\ \theta (r) = |\theta_0| n, \quad \text{si } r \leq R_{\inf} \\ \theta (r) = \frac{R_{\sup} - r}{R_{\sup} - R_{\inf}} |\theta_0| n, \quad \text{si } R_{\inf} \leq r \leq R_{\sup} \end{cases}\end{split}\]

Pour vérifier la condition \(\theta \cdot n = 0\) sur \(\delta \Omega\), il faut s’assurer que la couronne de rayon \(R_{\sup}\) soit strictement incluse dans \(\Omega\).

Fig. 311 Figure 3.3-a : Calcul du champ \(\theta\) en 2D

Autre méthode : champ \(\theta\) défini en fonction du nombre de couches

Alternativement, l’utilisation des opérandes NB_COUCHE_INF et NB_COUCHE_SUP permet de définir les couronnes d’intégration du champ \(\theta\) non pas en terme de rayon, mais en terme de nombre de couches d’éléments. La notion de couche est ici définie par la connectivité du maillage: la première couche d’éléments est constituée de tous les éléments connectés au front de fissure \({\Gamma}_{0}\) , la seconde couche d’éléments est constituée de tous les éléments connectés à la première couche, et ainsi de suite.

Fig. 312 Figure 3.4: Couches d’éléments 1 à 3 autour du front de fissure (2D)

Le champ \(\theta\) étant un champ aux nœuds, on dérive une notion de couche de nœuds à partir de la notion de couche d’éléments:

• Les nœuds du front de fissure se voient attribuer le numéro de couche 0

• Les nœuds sommets des éléments de la première couche d’éléments (hors nœuds du front) se voient attribuer le numéro de couche 1.

• Les nœuds sommets des éléments de la seconde couche d’éléments (hors nœuds de la première couche) se voient attribuer le numéro de couche 2.

• Et ainsi de suite, jusqu’à atteindre le numéro de couche NB_COUCHE_SUP. Tous les nœuds à l’extérieur de la couronne ainsi définie se voient attribuer le numéro de couche NB_COUCHE_SUP.

• Les éventuels nœuds milieux d’arêtes se voient attribuer un numéro de couche égal à la moyenne des numéros de couches de nœuds sommets de l’arête.

Une fois que chaque nœud du maillage a un numéro de couche, le champ \(\theta\) est calculé grâce aux méthodes présentées en section 3.2 (3D) et 3.3 (2D). R_SUP et R_INF sont alors simplement remplacés par NB_COUCHE_SUP et NB_COUCHE_INF, et la distance du nœud au front de fissure est remplacée par le numéro de couche du nœud.

Autre méthode

L’utilisateur peut spécifier directement le champ \(\theta\) en utilisant la commande CREA_CHAMP [U4.72.04].

Cette commande permet d’affecter le champ \(\theta\) :

• Nœud par nœud

• Par groupe de nœuds

Cela offre une flexibilité maximale pour définir un champ personnalisé en fonction des besoins spécifiques du modèle.

Bibliographie

• BUI H.D., Mécanique de la rupture fragile, Masson, 1977.

• DESTUYNDER Ph, DJAOUA M., Sur une interprétation de l’intégrale de Rice en théorie de la rupture fragile, Mathematics Methods in the Applied Sciences, Vol. 3, pp. 70-87, 1981.

• GRISVARD P., « Problèmes aux limites dans les polygones », Mode d’emploi - EDF- Bulletin de la Direction des Etudes et Recherches, Série C, 1, 1986 pp. 21-59.

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

• MIALON P., Etude du taux de restitution de l’énergie dans une direction marquant un angle avec une fissure, note interne EDF, HI/4740-07-1984.

• GURTIN M.E. An introduction to continuum mechanics. Mathematics in science and engineering. Academic Press, 1981.

• VENTURINI V. et al. Etude PRObabiliste de la cuve par un couplage Mecano-fiabilisTE. Bilan du projet P1-97-04, HP-26/99/012/A, nov. 1999.

• DHATT.G et TOUZOT.G. Une présentation de la méthodes des éléments finis. Ed. Maloine, 1984.

• MURAT.F et SIMON.J. Sur le contrôle par un domaine géométrique. Université de Paris VI, 1976.