r7.02.14 Éléments à discontinuité interne, comportement CZM_EXPet pilotage du chargement

Résumé :

En s’appuyant sur une formulation énergétique exposée brièvement dans une première partie nous présentons un modèle numérique d’amorçage et de propagation de fissures bidimensionnels. Celui-ci s’appuie sur un élément fini avec une discontinuité de déplacement interne, proche des modèles connus sous l’appellation «embedded discontinuityfinite element» dans la littérature.

Pour ce modèle, nous détaillons la loi de comportement utilisée (de type Cohésive Zone Model exponentielle: mot clé CZM_EXP ), les propriétés de l’élément fini ainsi que la résolution du problème de minimisation de l’énergie totale. Par ailleurs nous décrivons une technique de pilotage du chargement permettant de suivre d’éventuelles instabilités dans la réponse globale de la structure.

Ces éléments peuvent être utilisés en 2D plan ou axisymétrique uniquement avec les modélisations PLAN_ELDI et AXIS_ELDI. Ils sont validés sur les cas tests ssnp128a (documentation [V6.03.128]) et ssna115a (documentation [V6.01.115]).

Remarques :

• Cette documentation est largement inspirée des travaux de thèse [6] . Le lecteur intéressé pourra s’y reporter pour avoir une vision plus complète du lien de ce modèle avec l’approche énergétique «Francfort Marigo» ainsi qu’un certain nombre de résultats numériques.

• Le choix de regrouper dans une même documentation la présentation d’un élément fini et d’une loi de comportement est lié à la spécificité du modèle et a pour but de faciliter sa compréhension.

Table des matières

Modèle à discontinuité interne

Ce modèle permet de prendre en compte l’amorçage et la propagation de fissures dans une structure pour une direction donnée. Celui-ci est basé sur des élément finis particuliers que l’on appelle éléments à discontinuité interne (de l’anglais embedded discontinuity finite elements). L’idée principale repose sur l’introduction d’une discontinuité incluse dans l’élément et d’une variable interne de seuil gérant le processus dissipatif ainsi que le caractère irréversible de la fissuration. De plus, le saut de déplacement sera considéré constant par élément ce qui facilitera la résolution numérique. On pourra adopter une technique de condensation statique où le calcul du saut lors de la minimisation se fera au niveau élémentaire.

L’élément à discontinuité a été développé dans le but de s’affranchir de la régularisation de l’énergie en zéro (adhérence pénalisée) de la loi CZM_EXP_REG développée pour l’élément de joint (voir documentation [R3.06.09] pour l’élément et [R7.02.11] pour la loi). En effet, celle-ci est préjudiciable pour les modèles basés sur un principe de minimum d’énergie puisqu’elle conduit à annuler la dérivée de la densité d’énergie de surface en zéro c’est-à-dire la contrainte critique du modèle. Dans une telle situation, l’amorçage d’une fissure se produit dès la mise en charge, aussi faible soit elle (notons toutefois que cela n’empêche pas d’obtenir des résultats intéressant avec une telle loi, dès lors qu’on se fixe le trajet de fissuration).

Pour prendre en compte la discontinuité on effectue un enrichissement des déplacements et des déformations qui permet d’assurer leur compatibilité.

Dans cette partie on présentera dans un premier temps la loi de comportement qui gouverne l’élément à discontinuité. On sera amené à définir un vecteur contrainte dans l’élément qui servira lors de la minimisation de l’énergie. On détaillera ensuite les propriétés de l’élément fini puis la résolution numérique du problème de minimisation d’énergie.

Loi de comportement de l’élément à discontinuité : CZM_EXP

Les éléments à discontinuité ont un comportement mixte: élastique et dissipatif. L’énergie totale de l’élément s’écrit comme la somme d’une énergie élastique et d’une énergie de surface de type CZM. Les parties suivantes seront consacrées à la présentation de cette dernière ainsi qu’à la contrainte qui en dérive.

Énergie totale dans l’élément à discontinuité

On choisit l’énergie totale \({E}_{T}\) comme la somme d’une énergie élastique définie sur l’élément noté \({\Omega}_{e}\) et d’une énergie de surface définie sur la discontinuité de l’élément notée \({\Gamma}_{e}\) :

(4390)\[{E}_{T}(U,\delta ;\kappa )=\Phi (U,\delta )+\Psi (\delta ,\kappa )\]

où \(\Phi\) correspond à l’énergie élastique. Notons que celle-ci dépend du déplacement nodal \(U\) et du saut de déplacement

constant (dans l’élément) dont on verra qu’il induit une déformation supplémentaire. L’énergie de surface \(\Psi\) dépend du saut de déplacement et de \(\kappa\) variable interne permettant de traiter l’irréversibilité de la fissuration.

Remarque : Le trajet de fissuration est connu a priori , la normale à la discontinuité sera donc fixée par l’orientation des éléments dans le maillage (voir figure ).

Énergie élastique

Pour définir l’énergie élastique nous aurons besoin d’introduire dès à présent l’interpolation du champ de déplacement dans l’élément à discontinuité, les propriétés de ce dernier seront toutefois détaillées dans une partie à suivre (voir § 2.2 ). On représente sur la figure l’élément à discontinuité noté \({\Omega}_{e}\) , c’est un quadrangle à quatre nœuds avec une discontinuité interne que l’on note \({\Gamma}_{e}\) .

Figure 1: Élément à discontinuité

L’interpolation du champ de déplacement \(u(x)\) dans l’élément à discontinuité s’écrit:

(4391)\[\begin{split}u(x)=N(x).\left[\begin{array}{c}{U}_{A}\\ {U}_{B}\\ {U}_{C}-\delta \\ {U}_{D}-\delta \end{array}\right]+{H}_{{\Gamma}_{e}}(x).\delta\end{split}\]

avec \(N\) matrice des fonctions de forme bilinéaire classique de type \(\mathrm{Q1}\) , \({H}_{{\Gamma}_{e}}\) la fonction de Heaviside sur la discontinuité de l’élément \({\Gamma}_{e}\) valant 1 si \(x\in {\Omega}_{e}^{+}\) et 0 si \(x\in {\Omega}_{e}^{-}\) . Le vecteur \(\delta\) correspond au saut de déplacement sur la discontinuité. On a bien \(〚u〛=\delta\) puisque \(N\) est une fonction continue et \({H}_{{\Gamma}_{e}}=1\) sur \({\Gamma}_{e}\) . De plus, les \(U\) (indicés par les sommets) correspondent aux déplacements nodaux, en effet on peut écrire:

(4392)\[\begin{split}\left\lbrace \begin{array}{c}{U}_{A}\\ {U}_{B}\\ {U}_{C}\\ {U}_{D}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{U}_{A}\\ {U}_{B}\\ {U}_{C}+\delta -\delta \\ {U}_{D}+\delta -\delta \end{array}\right\rbrace\end{split}\]

A partir du champ approché () on peut définir la déformation associée sur l’élément hors de la discontinuité :

(4393)\[\begin{split}\epsilon ={\nabla}^{\text{s}}N(x).\left\lbrace \begin{array}{c}{U}_{A}\\ {U}_{B}\\ {U}_{C}-\delta \\ {U}_{D}-\delta \end{array}\right\rbrace +\underset{=0}{\underset{\underbrace{}}{{\nabla}^{\text{s}}{H}_{{\Gamma}_{e}}(x)}}\cdot \delta\end{split}\]

\({\nabla}^{\text{s}}\) désigne le gradient symétrisé. On peut récrire () sous la forme synthétique:

(4394)\[\epsilon =B.U-D.\delta\]

où \(B\) est la matrice des gradients symétrisés des fonctions de forme, \(U={({U}_{A},{U}_{B},{U}_{C},{U}_{D})}^{T}\) et \(D.\delta =B{(0,0,\delta ,\delta )}^{T}\) . Dans l’expression () on distingue une partie de la déformation liée aux déplacements \(B.U\) et une autre partie \(-D.\delta\) liée au saut. La contrainte pour une loi de comportement élastique s’écrit:

(4395)\[\sigma =E.\epsilon\]

Avec \(E\) le tenseur d’élasticité. L’énergie élastique dans l’élément est donc définie par:

(4396)\[\Phi (U,\delta )=\frac{1}{2}.{\int}_{{\Omega}_{e}}{(B.U-D.\delta )}^{T}.E.(B.U-D.\delta ).d\Omega\]

Énergie de surface

Comme nous l’avons présenté dans le modèle théorique (voir § 1 ), on choisit une énergie de surface dépendant de la norme du saut de déplacement \(\parallel \delta \parallel\) , de \(\kappa\) variable interne positive (sa valeur initiale \({\kappa}_{0}\) est nulle) et de l’orientation de la fissure (uniquement pour la condition de contact). On prend:

(4397)\[\Psi (\delta ,\kappa )=H(\parallel \delta \parallel -\kappa ).{\Psi}_{\mathrm{dis}}(\parallel \delta \parallel )+\left[1-H(\parallel \delta \parallel -\kappa )\right].{\Psi}_{\mathrm{lin}}(\parallel \delta \parallel ,\kappa )+{I}_{{ℝ}^{+}}(\delta \cdot n)\]

La fonction indicatrice \({I}_{{ℝ}^{+}}\) permet de prendre en compte la condition de non interpénétration des lèvres de la fissure:

(4398)\[\begin{split}{I}_{{ℝ}^{+}}({\delta}_{n})=\lbrace \begin{array}{cc}+\infty & \text{si}\delta \cdot n<0\\ 0& \text{si}\delta \cdot n\ge 0\end{array}\end{split}\]

Suivant la valeur du seuil, l’énergie de surface vaudra \({\Psi}_{\mathrm{dis}}\) ou \({\Psi}_{\mathrm{lin}}\) (plus la fonction indicatrice). Dans le premier cas on parlera de régime dissipatif, dans le second cas de régime linéaire. On peut écrire les deux énergies sous la forme:

(4399)\[{\Psi}_{\mathrm{dis}}(\parallel \delta \parallel )={\int}_{{\Gamma}_{e}}{\psi}_{\mathrm{dis}}(\parallel \delta \parallel ).d\Gamma\]

avec \({\psi}_{\mathrm{dis}}\) et \({\psi}_{\mathrm{lin}}\) densités d’énergie de surface. Présentons en détail les valeurs de ces densités.

Densité d’énergie de surface en régime linéaire

Dans le cas où une fissure existante évolue (ouverture ou refermeture) sans dissiper d’énergie i.e. \(\parallel \delta \parallel <\kappa\) , l’élément sera dans une phase linéaire (charge ou décharge). On choisit donc une densité d’énergie fonction quadratique de la norme du saut:

(4400)\[{\Psi}_{\mathrm{lin}}(\parallel \delta \parallel ,\kappa )=\frac{1}{2}.P(\kappa ).{\parallel \delta \parallel }^{2}+{C}_{0}\]

où \(P(\kappa )=\frac{{\sigma}_{c}}{\kappa}.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\kappa )\) est choisi afin d’assurer la continuité de la dérivée de \(\Psi\) en \(\kappa\) (i.e. la continuité de la contrainte dans l’élément d’un régime à l’autre) et \({C}_{0}\) constante permettant d’assurer la continuité de \(\Psi\) en \(\kappa\) .

Densité d’énergie en régime dissipatif

En s’inspirant de l’idée de Barenblatt [1] de tenir compte du processus de rupture des liaisons inter-atomiques, on suppose que l’énergie de surface est nulle en zéro et croît vers l’énergie de Griffith lorsque la valeur du saut devient importante devant la longueur caractéristique de l’échelle atomique. De plus en s’appuyant sur la forme des potentiels inter-atomiques on choisit une énergie croissante concave et à dérivée convexe. On sait en particulier que la concavité joue un rôle essentiel en dimension 1 pour limiter le nombre de fissures (voir Charlotte et al. [2] ). Ainsi, dans le cas où \(\parallel \delta \parallel <\kappa\) , on choisit une densité d’énergie de surface de la forme suivante, illustrée sur la figure :

(4401)\[{\Psi}_{\mathrm{dis}}(\parallel \delta \parallel )={G}_{\mathrm{c.}}\left[1-\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\parallel \delta \parallel )\right]\]

\({G}_{c}\) représente le taux de restitution d’énergie critique (ou ténacité) du matériau et \({\sigma}_{c}={\Psi}_{\mathrm{rev}}^{'}(0)\) la contrainte critique.

Figure 2: Densité d’énergie de surface en fonction de la norme du saut

Vecteur contrainte dans l’élément à discontinuité

Soit \(\overrightarrow{\sigma}=({\sigma}_{n},{\sigma}_{t})\) le vecteur contrainte dans l’élément. Lorsque l’élément est sain \((\kappa =0)\) , l’énergie de surface est nulle. Le vecteur contrainte est défini à partir du tenseur des contraintes de la loi de comportement élastique et de la normale à la ligne de discontinuité:

(4402)\[\overrightarrow{\sigma}=\sigma \cdot n\]

On verra par la suite (voir § 2.3.1.1 sur le critère d’amorçage ) que le saut reste nul tant que la norme du vecteur contrainte n’atteint pas la contrainte critique \({\sigma}_{c}\) . En revanche, si le seuil \(\kappa\) n’est pas nul, le vecteur contrainte dans l’élément est donné par la dérivée de la densité d’énergie de surface par rapport au saut(le vecteur contrainte reste égal à \(\sigma \cdot n\) puisque l’élément assure la continuité de la contrainte normale, voir § 2.2.2 ).

Régime linéaire

En régime linéaire (\(\parallel \delta \parallel <\kappa\) ), le vecteur contrainte vaut:

(4403)\[\overrightarrow{\sigma}=\frac{\partial {\Psi}_{\mathrm{lin}}}{\partial \delta }=P(\kappa ).\delta\]

Et sa dérivée par rapport au saut:

(4404)\[\frac{\partial \overrightarrow{\sigma}}{\partial \delta }=\mathrm{Id}.P(\kappa )\]

Régime dissipatif

En régime dissipatif (\(\parallel \delta \parallel >\kappa\) ), le vecteur contrainte vaut:

(4405)\[\overrightarrow{\sigma}=\frac{\partial {\Psi}_{\mathit{dis}}}{\partial \delta }={\sigma}_{c}.\frac{\delta}{\mathrm{\parallel }\delta \mathrm{\parallel }}.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\mathrm{\parallel }\delta \mathrm{\parallel })\]

Et sa dérivée par rapport au saut:

(4406)\[\frac{\partial \overrightarrow{\sigma}}{\partial \delta }={\sigma}_{c}.\left[\frac{\mathrm{Id}}{\parallel \delta \parallel }-\frac{\delta}{\parallel \delta \parallel }\otimes \frac{\delta}{\parallel \delta \parallel }.(\frac{{\sigma}_{c}}{{G}_{c}}+\frac{1}{\parallel \delta \parallel })\right].\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\parallel \delta \parallel )\]

Illustration graphique

On représente sur la figure l’évolution de la norme du vecteur contrainte dans l’élément en fonction de la norme du saut. Les flèches représentent les évolutions possibles du vecteur contrainte suivant le cas de figure (élément sain, régime linéaire ou régime dissipatif).

Figure 3: Norme du vecteur contrainte en fonction de la norme du saut

Au seuil en saut \(\kappa\) correspond un seuil en contrainte que l’on note \(R(\kappa )\) . Ce dernier va déterminer à partir de quel niveau de contrainte la fissure dissipera de l’énergie. Ce seuil va évoluer avec l’ouverture de la fissure, il dépend de

et est défini par:

(4407)\[R(\kappa )={\sigma}_{c}.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\kappa )\]

Notons que lorsque \(\kappa =0\) ce seuil correspond au critère d’amorçage que nous présenterons au § 2.3.1.1 .

Propriétés de l’élément à discontinuité

L’élément à discontinuité est un quadrangle à quatre nœuds avec un saut de déplacement interne. Une première partie sera consacrée à la description géométrique de l’élément ainsi qu’à la présentation d’un paramétrage en s’appuyant sur l’élément de référence. Ensuite, nous verrons que l’élément assure la continuité de la contrainte normale à travers le saut de déplacement. Puis, nous montrerons l’unicité du saut de déplacement pourvu que la taille de l’élément soit suffisamment petite. Pour finir, nous verrons que le choix d’un saut constant introduit une énergie de surface parasite qui tend vers zéro quand on raffine le maillage.

Géométrie de l’élément et paramétrage

Géométrie

Soit \((X,Y)\) la base cartésienne constituant un repère global de \({ℝ}^{2}\) . L’élément à discontinuité, noté \({\Omega}_{e}\) est un élément fini quadrangle à quatre nœuds. Il est constitué d’une région élastique et d’une discontinuité \({\Gamma}_{e}\) passant par le centre de l’élément (segment passant par les milieux des côtés \([\mathrm{AD}]\) et \([\mathrm{BC}]\) ) et de longueur \(l\) (voir figure ).

Figure 4: Élément à discontinuité L’orientation de la discontinuité définit un repère local à l’élément \((n,t)\) . L’élément de référence correspondant est un carré \(\stackrel{ˆ}{{\Omega}_{e}}\) défini par le domaine \([-1,1]\times [-1,1]\) (voir figure ). Chacun des éléments possède quatre points de Gauss. Ceux de l’élément de référence ont pour coordonnées \((\pm \sqrt{(3)}/3,\pm \sqrt{(3)}/3)\) . Notons \({\omega}_{g}\) et \(\stackrel{ˆ}{{\omega}_{g}}\) les poids des points de Gauss respectivement dans la configuration réelle et dans la configuration de référence. Enfin, rappelons que l’approximation du champ de déplacement dans l’élément a été présenté dans la partie§ 2.1.1.1 .

Figure 5: Élément de référence

Paramétrage

Figure 6: Paramétrage de l’élément à discontinuité Soit \(T\) la transformation géométrique qui à un point \((x,y)\) de l’élément à discontinuité associe un point \((\eta ,\xi )\) de élément de référence. Cherchons à paramétrer la position d’un point \(M\) quelconque de l’élément à discontinuité par \((\eta ,\xi )\) ses coordonnées dans le repère de l’élément de référence.

Soient \(P\) , \(Q\) et \(M\) les points (voir figure ) tels que:

(4408)\[\begin{split}\lbrace \begin{array}{c}\overrightarrow{\mathrm{AP}}=\frac{1}{2}.(1+\xi ).\overrightarrow{\mathrm{AB}}\\ \overrightarrow{\mathrm{DQ}}=\frac{1}{2}.(1+\xi ).\overrightarrow{\mathrm{DC}}\\ \overrightarrow{\mathrm{PM}}=\frac{1}{2}.(1-\eta ).\overrightarrow{\mathrm{PQ}}\end{array}\end{split}\]

En utilisant la somme des vecteurs, on a:

(4409)\[\overrightarrow{\mathrm{AM}}=\overrightarrow{\mathrm{AP}}+\overrightarrow{\mathrm{PM}}=\frac{1}{2}.(1+\xi ).\overrightarrow{\mathrm{AB}}+\frac{1}{2}.(1-\eta ).\overrightarrow{\mathrm{PQ}}\]

On ré-écrit le vecteur \(\overrightarrow{\mathrm{PQ}}\) :

(4410)\[\overrightarrow{\mathrm{PQ}}=\overrightarrow{\mathrm{PA}}+\overrightarrow{\mathrm{AD}}+\overrightarrow{\mathrm{DQ}}\]

En développant avec les expressions de ():

(4411)\[\overrightarrow{\mathrm{PQ}}=\frac{1}{2}.(1+\xi ).\overrightarrow{\mathrm{BA}}+\overrightarrow{\mathrm{AD}}+\frac{1}{2}.(1+\xi ).\overrightarrow{\mathrm{DC}}\]

En regroupant les termes, on obtient finalement:

(4412)\[\overrightarrow{\mathrm{AM}}=\frac{1+\eta }{2}.\frac{1+\xi }{2}.\overrightarrow{\mathrm{AB}}+\frac{1-\eta }{2}.\overrightarrow{\mathrm{AD}}+\frac{1-\eta }{2}.\frac{1+\xi }{2}.\overrightarrow{\mathrm{DC}}\]

Et donc:

(4413)\[\begin{split}dM=\left\lbrace \begin{array}{c}\mathrm{dx}\\ \mathrm{dy}\end{array}\right\rbrace =\frac{N(\xi )}{2}.d\eta +\frac{T(\eta )}{2}.d\xi\end{split}\]

Avec les deux vecteurs du repère sur la discontinuité:

(4414)\[\begin{split}\lbrace \begin{array}{c}N(\xi )=\left\lbrace \begin{array}{c}{N}_{x}\\ {N}_{y}\end{array}\right\rbrace =\frac{(1+\xi )}{2}.\overrightarrow{\mathrm{AB}}-\frac{(1-\xi )}{2}.\overrightarrow{\mathrm{DC}}-\overrightarrow{\mathrm{AD}}\\ T(\eta )=\left\lbrace \begin{array}{c}{T}_{x}\\ {T}_{y}\end{array}\right\rbrace =\frac{(1+\eta )}{2}.\overrightarrow{\mathrm{AB}}+\frac{(1-\eta )}{2}.\overrightarrow{\mathrm{DC}}\end{array}\end{split}\]

On a alors le paramétrage:

(4415)\[\begin{split}\left\lbrace \begin{array}{c}\mathrm{dx}\\ \mathrm{dy}\end{array}\right\rbrace =\left[J\right].\left\lbrace \begin{array}{c}d\eta \\ d\xi \end{array}\right\rbrace\end{split}\]

Où \(\left[J\right]\) est la matrice jacobienne de la transformation:

(4416)\[\begin{split}\left[J\right]=\frac{1}{2}.\left[\begin{array}{cc}{N}_{x}& {T}_{x}\\ {N}_{y}& {T}_{y}\end{array}\right]\end{split}\]

Remarque : La discontinuité se situant au centre de l’élément, on peut définir les vecteurs du repère local à celle-ci ainsi que sa longueur à partir du vecteur \(T\) en \(\eta =0\) :

(4417)\[\begin{split}\lbrace \begin{array}{c}n={R}_{-\pi /2}T(0)/\parallel T(0)\parallel \\ t=T(0)/\parallel T(0)\parallel \\ l=\parallel T(0)\parallel \end{array}\end{split}\]

Continuité de la contrainte normale

En s’appuyant sur le paramétrage précédent, on montre dans [6] comment l’élément à discontinuité assure la continuité de la contrainte normale à travers le saut de déplacement quand la contrainte est homogène dans l’élément. Dit autrement, on montre que le vecteur contrainte \(\overrightarrow{\sigma}={\psi}^{'}(\delta )\) de la loi de comportement CZM_EXP (voir § 2.1.2 ) est égal à la contrainte normale dans l’élément \(\sigma .n\) .

Condition d’existence et d’unicité du saut dans l’élément

D’un point de vue numérique il est important de se placer dans un cas de figure où la recherche du saut dans l’élément conduit à une solution unique. En d’autres termes il est nécessaire que la solution du problème de minimisation () possède une seule solution à \(U\) et \(\kappa\) fixés.

(4418)\[\underset{\delta \in {ℝ}^{2}}{\min}{E}_{T}(U,\delta ,\kappa )\]

Dans [6] on montre que l’existence et l’unicité du saut sont assurés dès que la condition suivante, portant sur la géométrie ainsi que sur les paramètres matériau de l’élément, est assurée:

(4419)\[\frac{\mu}{16}.\sum_{g}(\frac{1}{{\omega}_{g}}.\underset{{\eta}_{g}}{\min}{\parallel T({\eta}_{g})\parallel }^{2})>l.\frac{{\sigma}_{c}^{2}}{{G}_{c}}\]

Remarque : Dans le cas particulier où l’élément à discontinuité est rectangulaire, les poids des points de Gauss sont égaux à un quart de la surface de l’élément. De plus, la longueur de l’élément est égale à la longueur de la discontinuité \(l\) . Si on note \(e\) la largeur de l’élément, on a pour tout \(g\) :

(4420)\[{\omega}_{g}=\frac{l.e}{4}\]

De plus \(\parallel T\parallel =l\) sur chaque point de Gauss, donc la condition d’unicité devient une condition sur la largeur de l’élément:

(4421)\[e<\mu .\frac{{G}_{c}}{{\sigma}_{c}^{2}}\]

C’est cette dernière condition qui permet d’assurer l’unicité du saut pour les éléments à discontinuité dans Code_Aster . Il faut donc manipuler avec précaution ces derniers quand ils ne sont pas de forme rectangulaire et s’assurer que la condition donnée par () est bien vérifiée.

Énergie parasite

Le saut constant dans l’élément à discontinuité conduit à introduire une énergie parasite à l’interface entre deux éléments adjacents. On montre dans [6] que cette énergie tend vers zéro quand on raffine le maillage.

Minimisation de l’énergie totale

En adoptant le principe de minimisation d’énergie, le but de cette partie est de présenter le calcul du saut de déplacement dans les éléments à discontinuité ainsi que celui du champ de déplacement. L’énergie totale (voir § 2.1.1 ) n’est pas convexe vis-à-vis du couple \((U,\delta )\) . La recherche d’un minimum global pour une telle fonctionnelle n’est pas possible avec la méthode numérique que nous utiliserons. De plus, à force imposée, l’énergie totale n’étant pas bornée inférieurement, le minimum global n’existe pas. Ces deux arguments nous amènent a faire le choix de rechercher un minimum local. Le problème de minimisation de l’énergie totale s’écrit:

(4422)\[\text{Chercher}({U}^{\text{*}},{\delta}^{\text{*}})\text{minimum local de l'énergie totale}{E}_{T}(U,\delta ,\kappa )\]

Dans une première section (§ 2.3.1 ) on cherchera à calculer \({\delta}^{\text{*}}\) minimum local de l’énergie totale à \(U\) fixé:

(4423)\[{\delta}^{\text{*}}(U)=\underset{\delta}{\text{argmin}}{E}_{T}(U,\delta ,\kappa )\]

Dans la seconde section (§ 2.3.2 ), on se contentera de chercher une condition nécessaire pour que \({U}^{\text{*}}\) soit minimum local de l’énergie totale avec \(\delta ={\delta}^{\text{*}}(U)\) :

(4424)\[{U}^{\text{*}}=\underset{U}{\text{argmin}}{E}_{T}(U,{\delta}^{\text{*}}(U),\kappa )\]

Cette seconde étape revient à résoudre les équations d’équilibre en tenant compte des éventuels travaux des efforts extérieurs et des déplacements imposés à la structure.

Remarque : Le traitement numérique de l’indicatrice \({I}_{{ℝ}^{+}}(\delta \cdot n)\) traduisant la condition de non-interpénétration des lèvres de la fissure, s’effectuera en ajoutant la contrainte \(\delta \cdot n\ge 0\) au problème de minimisation.

Minimisation de l’énergie totale par rapport au saut

L’objet de cette section est de calculer les sauts de déplacement, sur chacun des éléments à discontinuité d’un maillage, qui minimisent l’énergie totale. Si on indice par \(i\) les éléments à discontinuité et que l’on note \({\delta}_{i}\) les sauts de déplacement sur chacun de ces éléments, le problème de minimisation s’écrit:

(4425)\[\mathrm{Chercher}\mathrm{les}{\delta}_{i}\text{minima locaux de l'énergie totale, solutions de}\underset{{\delta}_{i},{\delta}_{i}\cdot n\ge 0}{\min}\sum_{i}{E}_{T}(U,{\delta}_{i},\kappa )\]

Deux points importants vont permettre de simplifier sensiblement le problème. Tout d’abord, le choix d’un saut de déplacement indépendant d’un élément à l’autre permet de minimiser l’énergie totale par rapport au saut à un niveau élémentaire (technique de condensation statique). On a l’égalité suivante:

(4426)\[\underset{{\delta}_{i},{\delta}_{i}\cdot n\ge 0}{\min}\sum_{i}{E}_{T}(U,{\delta}_{i},\kappa )=\sum_{i}\underset{{\delta}_{i},{\delta}_{i}\cdot n\ge 0}{\min}{E}_{T}(U,{\delta}_{i},\kappa )\]

Par ailleurs, on a vu dans le § 2.2.3 que, pour des éléments suffisamment petits, le saut de déplacement est unique. Ainsi, le problème se ramène à une recherche de minimum global sur chaque élément:

(4427)\[\text{Sur chaque élément chercher}\delta \in S\text{minimum global de}{E}_{T}(U,{\delta}_{i},\kappa )\]

avec \(S\) l’ensemble des saut admissibles:

(4428)\[S=\lbrace s/s\text{constant par élément et}s\cdot n\ge 0\rbrace\]

La condition d’équilibre dite condition d’optimalité d’ordre un devient nécessaire et suffisante pour déterminer la solution du problème. Celle-ci s’écrit sous la forme:

(4429)\[{E}_{T}(U,\delta ,\kappa )\le {E}_{T}(U,\delta +h\phi ,\kappa )\]

Pour tout \(\phi\) constant par élément et pour \(h\) suffisamment petit de sorte que \(\delta +h\phi \in S\) . En passant à la limite quand \(h\) tend vers zéro, devient une condition sur la dérivée directionnelle:

(4430)\[\underset{h\to 0}{\lim}\frac{1}{h}.({E}_{T}(U,\delta +h\phi ,\kappa )-{E}_{T}(U,\delta ,\kappa ))\ge 0\]

Pour tout \(\phi\) constant par élément tel que \(\delta +h\phi \in S\) . D’une façon générale, la dérivée directionnelle est une fonction positivement homogène de degré un en \(\phi\) .

Par la suite nous présenterons le calcul du saut déplacement sur un élément, en s’appuyant sur la condition d’optimalité. On commencera par mettre en évidence une condition nécessaire et suffisante, que l’on appellera critère d’amorçage, pour que le saut nul soit solution du problème. Puis nous détaillerons le calcul du saut après amorçage, pour les deux types de comportement de l’élément: linéaire et dissipatif.

Critère d’amorçage

Dans cette partie nous chercherons à déterminer à quelle condition le saut de déplacement nul est solution du problème (). Notons que le choix de faire dépendre l’énergie de surface de la norme du saut implique que l’énergie totale n’admette pas de dérivée en zéro. De ce fait la dérivée de l’énergie dans la direction \(\phi\) ne sera pas une fonction linéaire mais positivement homogène de degré un en \(\phi\) . En utilisant la définition de l’énergie totale (voir § 2.1.1 ) lorsque l’élément est sain (i.e. quand \(\kappa =0\) ), on calcule sa dérivée dans la direction \(\phi\) en zéro:

(4431)\[\underset{h\to 0}{\lim}\frac{1}{h}.({E}_{T}(U,h\phi ,\kappa )-{E}_{T}(U,0,\kappa ))=-\sum_{g}{\omega}_{g}.{\phi }^{t}.{D}_{g}^{t}.E.{B}_{g}.U+l.{\sigma}_{c}.\parallel \phi \parallel\]

La condition d’optimalité () conduit donc à l’inégalité:

(4432)\[-\sum_{g}{\omega}_{g}.{\phi }^{t}.{D}_{g}^{t}.E.{B}_{g}.U+l.{\sigma}_{c}.\parallel \phi \parallel \ge 0\]

Pour tout \(\phi \in S\) . Or, pour montrer la continuité de la contrainte normale (voir [6] ), on a vu que:

(4433)\[\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.{B}_{g}.U=\mathrm{l.}\overrightarrow{\sigma}\]

De plus \(\phi\) est constant sur l’élément, donc () devient:

(4434)\[-\overrightarrow{\sigma}\cdot \phi +{\sigma}_{c}.\parallel \phi \parallel \ge 0\]

En posant \(\lbrace \begin{array}{c}{\phi }_{n}=\rho .\cos\theta \\ {\phi }_{t}=\rho .\sin\theta \end{array}\) avec \(\rho >0\) et \(\frac{-\pi }{2}\le \theta \le \frac{\pi}{2}\) afin que \({\phi }_{n}\ge 0\) , et en utilisant la définition \(\overrightarrow{\sigma}=({\sigma}_{n},{\sigma}_{t})\) , on obtient la condition:

(4435)\[{\sigma}_{n}.\cos\theta +{\sigma}_{t}.\sin\theta \le {\sigma}_{c}\forall \theta \in \left[-\frac{\pi}{2},\frac{\pi}{2}\right]\]

Une étude de la fonction \(f(\theta )={\sigma}_{n}.\cos\theta +{\sigma}_{t}.\sin\theta\) donne:

(4436)\[\underset{-\frac{\pi}{2}\le \theta \le \frac{\pi}{2}}{\sup}f(\theta )=\sqrt{\text{<}{\sigma}_{n}{\text{>}}_{+}^{2}+{\sigma}_{t}^{2}}\]

Avec la définition de la partie positive d’une quantité \(\text{<}.{\text{>}}_{+}=\max(.,0)\) . L’inégalité () conduit donc au critère d’amorçage en contrainte suivant:

(4437)\[\sqrt{\text{<}{\sigma}_{n}{\text{>}}_{+}^{2}+{\sigma}_{t}^{2}}\le {\sigma}_{c}\]

On représente ce critère dans le plan \(({\sigma}_{n},{\sigma}_{t})\) sur la figure .

Figure 7: Critère d’amorçage en contrainte de l’élément à discontinuité

Remarques : On se place dans le cas de figure où l’orientation de la discontinuité \(n\) est fixée (par le maillage). Cela explique que l’on obtient une condition sur le vecteur contrainte et non sur les composantes principales de la contrainte.

Calcul du saut

Cherchons à résoudre le problème de minimisation dans le cas où le saut n’est pas nul. Dans ce cas, la dérivée directionnelle est une fonction linéaire de \(\phi\) , la condition d’optimalité () s’écrit:

(4438)\[\underset{h\to 0}{\lim}\frac{1}{h}.({E}_{T}(U,\delta +h\phi ,\kappa )-{E}_{T}(U,\delta ,\kappa ))=\frac{\partial {E}_{T}}{\partial \delta }\cdot \phi \ge 0\]

Pour tout \(\phi\) constant par élément tel que \(\delta +h\phi \in S\) . En jouant sur les directions admissibles on en déduit les conditions:

(4439)\[\begin{split}\lbrace \begin{array}{ccc}\frac{\partial {E}_{T}}{\partial {\delta}_{n}}.{\phi }_{n}\ge 0& \forall {\phi }_{n}\ge 0& \text{si}{\delta}_{n}=0\\ \frac{\partial {E}_{T}}{\partial {\delta}_{n}}.{\phi }_{n}=0& \forall {\phi }_{n}& \text{si}{\delta}_{n}>0\\ \frac{\partial {E}_{T}}{\partial {\delta}_{t}}.{\phi }_{t}=0& \forall {\phi }_{t}& \end{array}\end{split}\]

Par ailleurs, en approchant l’intégrale par une somme discrète sur les points de Gauss, la dérivée de l’énergie totale dans la direction \(\phi\) s’écrit:

(4440)\[\frac{\partial {E}_{T}}{\partial \delta }\cdot \phi =-\sum_{g}{\omega}_{g}.{\phi }^{t}.{D}_{g}^{t}.E.({B}_{g}.U-{D}_{g}\delta ).\phi +l.\frac{\partial \psi }{\partial \delta }\cdot \phi\]

Avec \(\psi\) la densité d’énergie de surface, égale à \({\psi}_{\mathrm{lin}}\) en régime linéaire et à \({\psi}_{\mathrm{dis}}\) en régime dissipatif. Pour alléger les notations on définit \(S\) de la façon suivante:

(4441)\[\begin{split}S=\left\lbrace \begin{array}{c}{S}_{n}\\ {S}_{t}\end{array}\right\rbrace =-\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.{B}_{g}.U\end{split}\]

et \(Q\) :

(4442)\[\begin{split}Q=\left[\begin{array}{cc}{Q}_{\mathrm{nn}}& {Q}_{\mathrm{nt}}\\ {Q}_{\mathrm{nt}}& {Q}_{\mathrm{tt}}\end{array}\right]=-\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.{D}_{g}\end{split}\]

On a donc l’écriture simplifiée de ():

(4443)\[\frac{\partial {E}_{T}}{\partial \delta }\cdot \phi =(S+Q.\delta )\cdot \phi +l.\frac{\partial \psi }{\partial \delta }\cdot \phi\]

Cherchons maintenant à calculer le saut déplacement à partir des conditions en utilisant la dérivée de l’énergie écrite sous la forme (). On distinguera le calcul dans le régime linéaire de celui dans le régime dissipatif et pour chacun d’entre eux on distinguera le cas ou le saut normal est nul de celui où il ne l’est pas.

Calcul du saut en régime linéaire

En régime linéaire la densité d’énergie de surface est donnée par:

(4444)\[\psi ={\psi}_{\mathrm{lin}}=\frac{1}{2}.P(\kappa ).\delta \cdot \delta\]

La dérivée de l’énergie totale par rapport au saut () devient:

(4445)\[\frac{\partial {E}_{T}}{\partial \delta }\cdot \phi =(S+Q.\delta )\cdot \phi +l.P(\kappa ).\delta \cdot \phi\]

Dans le cas d’un saut normal nul \((0,{\delta}_{t})\) , les conditions d’optimalité () conduisent à:

(4446)\[\begin{split}\lbrace \begin{array}{c}{S}_{n}+{Q}_{\mathrm{nt}}.{\delta}_{t}\ge 0\\ {\delta}_{t}=-\frac{{S}_{t}}{{Q}_{t}+l.P(\kappa )}\end{array}\end{split}\]

Le saut tangent est donné explicitement par la seconde condition. La première impose que la contrainte normale dans l’élément soit négative ou nulle. Si tel n’était pas le cas, il y aurait nécessairement un saut normal. Dit autrement, le saut dans l’élément sera nul dès que l’élément sera mis en compression. Ceci traduit la prise en compte de la non-interpénétration des lèvres de la fissure.

Dans le cas d’un saut normal strictement positif \(({\delta}_{n}>0,{\delta}_{t})\) , les conditions d’optimalité () donnent directement le saut de déplacement:

(4447)\[\delta =-{(Q+\mathrm{Id}.l.P(\kappa ))}^{-1}.S\]

En régime linéaire la variable interne \(\kappa\) n’évolue pas.

Calcul du saut en régime dissipatif

En régime dissipatif la densité d’énergie de surface est donnée par:

(4448)\[\psi ={\psi}_{\mathrm{dis}}={G}_{c}.\left[1-\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\parallel \delta \parallel )\right]\]

La dérivée de l’énergie totale par rapport au saut () devient:

(4449)\[\frac{\partial {E}_{T}}{\partial \delta }\cdot \phi =(S+Q.\delta )\cdot \phi +l.{\sigma}_{c}.\frac{\delta \cdot \phi }{\parallel \delta \parallel }.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\parallel \delta \parallel )\]

Dans le cas d’un saut normal nul \((0,{\delta}_{t})\) , les conditions d’optimalité () conduisent à:

(4450)\[\begin{split}\lbrace \begin{array}{c}{S}_{n}+{Q}_{\mathrm{nt}}.{\delta}_{t}\ge 0\\ {S}_{t}+{Q}_{\mathrm{tt}}.{\delta}_{t}+l.{\sigma}_{c}.\text{sign}({\delta}_{t}).\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.∣{\delta}_{t}∣)=0\end{array}\end{split}\]

La première condition implique que le saut normal soit nul quand la contrainte normale dans l’élément est négative. Ce qui traduit la prise en compte de la non-interpénétration des lèvres de la fissure. De la seconde condition, on déduit que \(\text{sign}({\delta}_{t})=-\text{sign}({S}_{t})\) , donc \({\delta}_{t}=-\text{sign}({S}_{t}).\beta\) avec \(\beta >0\) est solution de l’équation scalaire non linéaire suivante:

(4451)\[∣{S}_{t}∣+{Q}_{n}.\beta +l.{\sigma}_{c}.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\beta )=0\]

que l’on résout par un algorithme de Newton.

Dans le cas d’un saut normal strictement positif \(({\delta}_{n}>0,{\delta}_{t})\) , les conditions d’optimalité () conduisent à l’équation non linéaire suivante:

(4452)\[S+Q.\delta +l.{\sigma}_{c}.\frac{{\delta}_{t}}{\parallel \delta \parallel }.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\parallel \delta \parallel )=0\]

On note \(\delta =r.\tilde{\delta}\) avec \(r>0\) et \(\tilde{\delta}=(\cos\theta ,\sin\theta )\) où \(\frac{-\pi }{2}\le \theta \le \frac{\pi}{2}\) de façon à ce que \(\tilde{\delta}\cdot n\ge 0\) . D’après () on a:

(4453)\[\tilde{\delta}=-{(r.Q+l.{\sigma}_{c}.\mathrm{Id}.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.r))}^{-1}.S\]

Il reste à calculer \(r\) . Pour cela il suffit de résoudre l’équation scalaire non linéaire \(\parallel \tilde{\delta}\parallel =1\) par un algorithme de Newton. En régime dissipatif la variable interne \(\kappa\) évolue, elle correspond à la norme du saut calculée \(\kappa =\parallel \tilde{\delta}\parallel\) .

Algorithme de calcul

Dans cette partie, on présente l’algorithme de calcul du saut \({\delta}^{\text{*}}\) sur un élément donné ainsi que l’évolution de la variable interne \(\kappa\) .

Début du calcul Test d’existence et d’unicité de la solution (voir §:ref:2.2.3 <Ref78379485> ) 3. Si (\(\kappa =0\) et \(\sqrt{\text{<}{\sigma}_{n}{\text{>}}_{+}^{2}+{\sigma}_{t}^{2}}\le {\sigma}_{c}\) ) alors * \({\delta}^{\text{*}}=0\) , le seuil \(\kappa\) reste nul → Fin du calcul 1. Si \(\kappa >0\) alors * Calcul du seuil en régime linéaire (voir 2.3.1.2 )→ \({\delta}^{c}\) * Si \(\parallel {\delta}^{c}\parallel <\kappa\) , alors \({\delta}^{\text{*}}={\delta}^{\text{c}}\) , le seuil \(\kappa\) n’évolue pas → Fin du calcul 1. Calcul du saut en régime dissipatif (voir §:ref:2.3.1.2 <Ref78379665> ) → \({\delta}^{c}\) * \({\delta}^{\text{*}}={\delta}^{\text{c}}\) , actualisation du seuil \(\kappa =\parallel {\delta}^{\text{*}}\parallel\) → Fin du calcul 1. Fin du calcul

Minimisation de l’énergie totale par rapport aux déplacements

Après avoir calculé les sauts de déplacement \({\delta}^{\text{*}}(U)\) sur chacun des éléments à discontinuité, l’objectif est de calculer le champ de déplacements qui minimise l’énergie totale. Le problème se formule de la façon suivante:

(4454)\[\text{Chercher}{U}^{\text{*}}\text{minimum local de l'énergie totale}{E}_{T}(U,{\delta}^{\text{*}}(U),\kappa )\]

La fonctionnelle \({E}_{T}(U,{\delta}^{\text{*}}(U),\kappa )\) n’est pas convexe en \(U\) . Pour être minimum local, la solution doit vérifier les conditions d’équilibre et de stabilité. La seconde condition porte sur la dérivée seconde de la fonctionnelle et nécessite l’utilisation d’une méthode spécifiquement dédiée à la minimisation. L’algorithme de Newton utilisé ne permet pas de vérifier cette seconde condition. On se contentera donc de vérifier la condition d’équilibre, condition nécessaire pour que le déplacement soit solution du problème ().

Considérons une structure \(\Omega\) avec le chargement suivant:

• \(f\) densité de force volumique sur \(\Omega\) ;

• \(F\) densité de force surfacique sur \({\Gamma}_{N}\) ;

• \({U}^{d}\) déplacements imposés sur \({\Gamma}_{D}\)

\({\Gamma}_{D}\) et \({\Gamma}_{N}\) sont des parties disjointes de la frontière de \(\Omega\) . Le travail des efforts extérieurs s’écrit:

(4455)\[{W}^{\mathrm{ext}}(U)={\int}_{\Omega}f.U.d\Omega +{\int}_{{\Gamma}_{N}}F.U.d{\Gamma}_{N}\]

L’énergie totale de la structure s’écrit alors:

(4456)\[{E}_{T}(U,{\delta}^{\text{*}}(U),\kappa )=\Phi (U,{\delta}^{\text{*}}(U))+\Psi ({\delta}^{\text{*}}(U),\kappa )-{W}^{\mathrm{ext}}(U)\]

Avec \(U\) appartenant à l’espace des champs de déplacement cinématiquement admissibles. La condition d’optimalité d’ordre un, condition nécessaire pour que \(U\) soit minimum local, s’écrit:

(4457)\[\frac{\partial {E}_{T}}{\partial U}(U,{\delta}^{\text{*}}(U),\kappa ).V\ge 0\]

Pour tout champ test \(V\) admissible. On sait que :

(4458)\[\frac{\partial {E}_{T}}{\partial \delta }(U,{\delta}^{\text{*}}(U),\kappa )=0\]

Ce qui donne:

(4459)\[{\int}_{\Omega}{B}^{t}.E.(B.U-D.{\delta}^{\text{*}}(U)).V.d\Omega -{\int}_{\Omega}f.V.d\Omega -{\int}_{{\Gamma}_{N}}F.V.d{\Gamma}_{N}\ge 0\forall V\]

Cette inégalité devient une égalité en prenant des \(V\) bien choisis, et, comme cette égalité est vraie pour tout \(V\) admissible, on obtient l’équilibre entre les efforts intérieurs et extérieurs:

(4460)\[{F}^{int}(U)={F}^{\text{ext}}\]

Avec:

(4461)\[\begin{split}\lbrace \begin{array}{c}{F}^{int}(U)={\int}_{\Omega}{B}^{t}.E.(B.U-D.{\delta}^{\text{*}}).d\Omega \\ {F}^{\text{ext}}={\int}_{\Omega}f.d\Omega +{\int}_{{\Gamma}_{N}}F.d{\Gamma}_{N}\end{array}\end{split}\]

Notons \(C\) l’opérateur linéaire traduisant les conditions de déplacement imposés. La dualisation de ces conditions conduit au système suivant:

(4462)\[\begin{split}\lbrace \begin{array}{c}{F}^{int}(U)+{C}^{t}.\lambda ={F}^{\text{ext}}\\ C.U={U}^{d}\end{array}\end{split}\]

Les inconnues sont à présent à tout instants le couple \((U,\lambda )\) où \(\lambda\) représente les multiplicateurs de Lagrange associés aux conditions de Dirichlet. La méthode utilisée pour résoudre est une méthode de Newton (voir documentation de STAT_NON_LINE [R5.03.01}). Le calcul de la matrice tangente est expliqué en annexe.

Méthode de pilotage du chargement

La méthode de pilotage du chargement a pour but de prendre en compte les instabilités de la structure, et suivre ainsi d’éventuels «retour arrière» de la réponse globale force déplacement. Un cas particulier de cette méthode, initialement développé par Lorentz et Badel [7] , a été adapté au modèle de fissuration précédent.

Le principe de la méthode de pilotage ainsi que la résolution du nouveau système global qui en découle est expliqué en détail dans [R5.03.80]. Nous détaillerons la partie du pilotage spécifique à la loi CZM_EXP en expliquant le choix de l’équation de contrôle du pilotage ainsi que la méthode pour la résoudre.

Équation de contrôle du pilotage

Le but de cette partie est de présenter l’équation de contrôle du pilotage ainsi que sa résolution pour le modèle à discontinuité interne CZM_EXP. L’inconnue de cette équation est l’intensité du chargement \({\eta}_{i}^{n}\) à l’instant de calcul \(i\) et à l’itération de Newton \(n\) que l’on notera dorénavant

pour simplifier la notation. Comme nous l’avons vu au § 2.1 la loi de comportement de l’élément à discontinuité est gouvernée par un seuil. Notons

la fonction seuil en contrainte suivante :

(4463)\[{F}_{\mathit{el}}({\sigma}_{n},{\sigma}_{t})=\sqrt{\text{<}{\sigma}_{n}{\text{>}}_{+}^{2}+{\sigma}_{t}^{2}}-R(\kappa )\]

Avec le vecteur saut de contrainte qui vaut:

(4464)\[\begin{split}\overrightarrow{\sigma}=\left\lbrace \begin{array}{c}{\sigma}_{n}\\ {\sigma}_{t}\end{array}\right\rbrace =\frac{1}{l}.\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.({B}_{g}.U-{D}_{g}\delta )\end{split}\]

Et \(R(\kappa )\) la variable de seuil en contrainte définie par ():

(4465)\[R(\kappa )={\sigma}_{c}.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\kappa )\]

Notons \({G}_{\mathrm{el}}\) la fonction seuil en saut:

(4466)\[{G}_{\mathit{el}}({\delta}_{n},{\delta}_{t})=\sqrt{{\delta}_{n}^{2}+{\delta}_{t}^{2}}-\kappa\]

Avec \({\delta}_{n}\ge 0\) et \(\kappa\) le seuil en saut. Choisissons une équation de contrôle du pilotage de telle sorte que l’intensité du chargement \(\eta\) fasse sortir du critère au moins un élément à discontinuité de la structure d’une quantité proportionnelle à \(\Delta \tau\) . Les fonctions seuils dépendent alors de \(\eta\) et nous les noterons \(\tilde{{F}_{\mathrm{el}}}(\eta )\) et \(\tilde{{G}_{\mathrm{el}}}(\eta )\) . L’équation de contrôle du pilotage s’écrit:

• Pour le critère en contrainte:

(4467)\[\tilde{P}(\eta )\underset{\text{def}}{=}\underset{j}{max}\tilde{{F}_{\mathit{el}}^{j}}(\eta )=R({\kappa}^{-}).\Delta \tau\]

• Pour le critère en saut :

(4468)\[\tilde{P}(\eta )\underset{\text{def}}{=}\underset{j}{max}\tilde{{G}_{\mathit{el}}^{j}}(\eta )={\kappa}^{-}.\Delta \tau\]

Où \(j\) désigne l’indice des éléments à discontinuité du maillage et \({\kappa}^{-}\) la variable seuil à l’instant \(i-1\) .

Remarques :

• Dans le cas où \({\kappa}^{-}=0\) , le saut étant nul, nous utiliserons l’équation de contrôle du pilotage exprimé avec le critère en contrainte ().

• Dans le cas où , \({\kappa}^{-}>0\) nous utiliserons l’équation de contrôle du pilotage exprimé avec le critère en saut (). En effet, le saut n’étant pas connu, on ne peut utiliser les contraintes \(({\sigma}_{n},{\sigma}_{t})\) qui dépendent explicitement de celui-ci.

Résolution de l’équation avec le critère en contrainte

Nous nous plaçons donc dans le cas où \({\kappa}^{-}=0\) .Le but ici est d’expliquer la résolution de () pour un élément donné, sans tenir compte du \(max\) sur tous les éléments (la prise en compte du \(max\) est expliquée dans [R5.03.80]). On cherche donc à résoudre:

(4469)\[\tilde{{F}_{\mathit{el}}}(\eta )=R(0).\Delta \tau\]

Faisons d’abord apparaître explicitement \(\eta\) dans l’expression du critère avant \(\tilde{{F}_{\mathrm{el}}}(\eta )\) d’expliquer la résolution de l’équation. A l’instant \(i\) on cherche le déplacement \({U}_{i}\) qui s’exprime à l’itération de Newton \(n\) :

(4470)\[{\mathrm{U}}_{i}={\mathrm{U}}_{i-1}+\Delta {\mathrm{U}}_{i}^{n}+\delta {\mathrm{U}}_{\text{impo},i}^{n}+{\eta}_{i}.\delta {\mathrm{U}}_{\text{pilo},i}^{n}\]

On sépare la partie connue (au sens de non-pilotéee) du déplacement:

(4471)\[{\mathrm{U}}_{\text{impo}}={\mathrm{U}}_{i-1}+\Delta {\mathrm{U}}_{i}^{n}+\delta {\mathrm{U}}_{\text{impo},i}^{n}\]

Et la partie pilotée:

(4472)\[{\mathrm{U}}_{\text{pilo}}=\delta {\mathrm{U}}_{\text{pilo},i}^{n}\]

Directement:

(4473)\[{\mathrm{U}}_{i}={\mathrm{U}}_{\text{impo}}+\eta .{\mathrm{U}}_{\text{pilo}}\]

Le saut de déplacement étant nul on a:

(4474)\[\begin{split}\left\lbrace \begin{array}{c}{\sigma}_{n}\\ {\sigma}_{t}\end{array}\right\rbrace =\frac{1}{l}.\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.{B}_{g}.{U}_{i}\end{split}\]

En utilisant la décomposition ():

(4475)\[\begin{split}\left\lbrace \begin{array}{c}{\sigma}_{n}\\ {\sigma}_{t}\end{array}\right\rbrace ={\mathrm{S}}_{\text{impo}}+\eta .{\mathrm{S}}_{\text{pilo}}\end{split}\]

Avec:

(4476)\[{\mathrm{S}}_{\text{impo}}=\frac{1}{l}.\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.{B}_{g}.{U}_{\text{impo}} {\mathrm{S}}_{\text{pilo}}=\frac{1}{l}.\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.{B}_{g}.{U}_{\text{pilo}}\]

L’équation de contrôle du pilotage s’écrit donc, en tenant compte du fait que \({\kappa}^{-}=0\) :

(4477)\[\tilde{{F}_{\mathit{el}}}(\eta )\underset{\text{def}}{=}\sqrt{{\langle ({\mathrm{S}}_{\text{impo}}+\eta .{\mathrm{S}}_{\text{pilo}})\cdot n\rangle }_{+}^{2}+{(({\mathrm{S}}_{\text{impo}}+\eta .{\mathrm{S}}_{\text{pilo}})\cdot t)}^{2}}-{\sigma}_{c}={\sigma}_{c}.\Delta \tau\]

Cette équation correspond à un polynôme de degré deux, que l’on note \({p}_{1}\) , si \(({\mathrm{S}}_{\text{impo}}+\eta .{\mathrm{S}}_{\text{pilo}})\cdot n>0\) et à un polynôme de degré deux que l’on note \({p}_{2}\) dans le cas contraire.

Figure 8: Sortie du critère en contrainte

Pour trouver \(\eta\) , on résout ces deux polynômes. Cela revient, dans l’espace \(({\sigma}_{n},{\sigma}_{t})\) , à trouver l’intersection de deux demi-droites avec le critère en contrainte. La figure () représente la «sortie de critère» contrôlée par le paramètre \(\Delta \tau\) .

Chacun de ces deux polynômes possède 0, 1 ou 2 solutions, on admettra les solutions de \({p}_{1}\) qui vérifient \(({\mathrm{S}}_{\text{impo}}+\eta .{\mathrm{S}}_{\text{pilo}})\cdot n>0\) et les solutions de \({p}_{2}\) qui vérifient \(({\mathrm{S}}_{\text{impo}}+\eta .{\mathrm{S}}_{\text{pilo}})\cdot n\le 0\) . Les solutions admises seront au plus deux, on les note \({\eta}_{k=0,1,2}\) .

Résolution de l’équation avec le critère en saut

Le but ici est de résoudre l’équation de contrôle du pilotage pour un élément donné dans le cas où \({\kappa}^{-}>0\) , toujours sans tenir compte du \(max\) sur tous les éléments. On cherche à résoudre:

(4478)\[\tilde{{G}_{\mathit{el}}}(\eta )={\kappa}^{-}.\Delta \tau\]

Faisons d’abord apparaître explicitement \(\eta\) dans l’expression du critère \(\tilde{{G}_{\mathrm{el}}}(\eta )\) puis détaillons le calcul de (). On procède comme dans le cas du critère en contrainte (§ 3.1.1 ). La différence concerne l’écriture du saut de contrainte dans le cas général (quand le saut n’est pas nul):

(4479)\[\begin{split}\left\lbrace \begin{array}{c}{\sigma}_{n}\\ {\sigma}_{t}\end{array}\right\rbrace =\frac{1}{l}.\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.({B}_{g}.{U}_{i}-{D}_{g}.\delta )\end{split}\]

Que l’on décompose en trois parties:

(4480)\[\begin{split}\left\lbrace \begin{array}{c}{\sigma}_{n}\\ {\sigma}_{t}\end{array}\right\rbrace ={\mathrm{S}}_{\text{impo}}+\eta .{\mathrm{S}}_{\text{pilo}}+Q.\delta\end{split}\]

Avec:

(4481)\[{\mathrm{S}}_{\text{impo}}=\frac{1}{l}.\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.{B}_{g}.{U}_{\text{impo}} {\mathrm{S}}_{\text{pilo}}=\frac{1}{l}.\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.{B}_{g}.{U}_{\text{pilo}} \mathrm{Q}=\frac{1}{l}.\sum_{g}{\omega}_{g}.{D}_{g}^{t}.E.{D}_{g}\]

Par ailleurs pour un seuil en saut \({\kappa}^{-}>0\) fixé on a, d’après la loi de comportement CZM_EXP:

(4482)\[\begin{split}\left\lbrace \begin{array}{c}{\sigma}_{n}\\ {\sigma}_{t}\end{array}\right\rbrace =P({\kappa}^{-}).\delta\end{split}\]

Avec \(P({\kappa}^{-})=\frac{{\sigma}_{c}}{{\kappa}^{-}}.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.{\kappa}^{-})\) . De () et () on déduit:

(4483)\[{\mathrm{S}}_{\text{impo}}+\eta .{\mathrm{S}}_{\text{pilo}}+Q.\delta =P({\kappa}^{-}).\delta\]

Ce qui permet de connaître \(\delta\) en fonction du paramètre de pilotage \(\eta\) :

(4484)\[\begin{split}\delta =\left\lbrace \begin{array}{c}{\delta}_{n}\\ {\delta}_{t}\end{array}\right\rbrace ={\delta}_{\text{impo}}+\eta .{\delta}_{\text{pilo}}\end{split}\]

Avec:

(4485)\[{\delta}_{\text{impo}}={(P({\kappa}^{-}).\mathit{Id}-Q)}^{\text{-1}}.{S}_{\text{impo}} {\delta}_{\text{pilo}}={(P({\kappa}^{-}).\mathit{Id}-Q)}^{\text{-1}}.{S}_{\text{pilo}}\]

L’équation de contrôle du pilotage s’écrit donc:

(4486)\[\tilde{{G}_{\mathit{el}}}(\eta )\underset{\text{def}}{=}\sqrt{{(({\delta}_{\text{impo}}+\eta .{\delta}_{\text{pilo}})\cdot n)}^{2}+{(({\delta}_{\text{impo}}+\eta .{\delta}_{\text{pilo}})\cdot t)}^{2}}-{\kappa}^{-}={\kappa}^{-}.\Delta \tau\]

Avec la condition de saut normal positif:

(4487)\[({\delta}_{\text{impo}}+\eta .{\delta}_{\text{pilo}})\cdot n\ge 0\]

Le critère se réduit à un demi cercle dans l’espace \(({\delta}_{n},{\delta}_{t})\) . Pour pouvoir utiliser le pilotage lorsque le saut normal tend à être négatif (compression de l’élément) nous allons autoriser une petite sortie de critère pour de tels sauts. La «sortie de critère», contrôlée par le paramètre \(\Delta \tau\) , peut être représentée dans l’espace \(({\delta}_{n},{\delta}_{t})\) par un arc de cercle et un segment (voir figure ).

Figure 9: Sortie du critère en saut

Remarque : Il est entendu que seule la prédiction de l’algorithme de Newton sortira du critère, donc en particulier pourrait violer la condition de contact. La solution du problème, quant à elle, la respectera.

La résolution de l’équation de contrôle du pilotage revient donc à trouver l’intersection entre une droite et un arc de cercle ou l’intersection entre une droite et un segment.

En résumé, on résout:

• Le polynôme (), de degré deux en \(\eta\) (intersection droite/cercle): on admet les \(\eta\) qui vérifient \(({\delta}_{\text{impo}}+\eta .{\delta}_{\text{pilo}})\cdot n\ge -{\kappa}^{-}.\Delta \tau\) ;

• Le polynôme \(({\delta}_{\text{impo}}+\eta .{\delta}_{\text{pilo}})\cdot n=-{\kappa}^{-}.\Delta \tau\) , de degré un en \(\eta\) : \(({\delta}_{p}+\eta {\delta}_{d})\cdot n=-{\kappa}^{-}\Delta \tau\) (intersection droite/segment), on admet les \(\eta\) qui vérifient \(∣({\delta}_{\text{impo}}+\eta .{\delta}_{\text{pilo}})∣\cdot t\le d\) avec \(d={\kappa}^{-}.\sqrt{1+2.\Delta \tau }\) .

On notera \({\eta}_{k=0,1,2}\) les solutions admises. S’agissant du choix parmi ces solutions on se reportera à [R5.03.80].

ANNEXE : Calcul de la matrice tangente pour l’élément à discontinuité

Dans cette partie nous allons détailler le calcul de la dérivée des efforts intérieurs par rapport au déplacement: \(\frac{\partial {F}^{int}}{\partial U}\) qui apparaît dans \({K}_{T}\) matrice tangente du système intervenant dans l’algorithme de Newton. On a:

(4488)\[{F}^{int}(U)={\int}_{\Omega}{B}^{t}.E.(B.U-D.{\delta}^{\text{*}}).d\Omega\]

Donc:

(4489)\[\frac{\partial {F}^{int}(U)}{\partial U}={\int}_{\Omega}{B}^{t}.E.(B-D.\frac{\partial {\delta}^{\text{*}}}{\partial U}).d\Omega\]

Il faut calculer le terme \(\frac{\partial {\delta}^{\text{*}}}{\partial U}\) sur chaque élément à discontinuité \({\Omega}_{e}\) (condensation statique). La minimisation par rapport à \(\delta\) de l’énergie totale nous a conduit à calculer \({\delta}^{\text{*}}\) solution de l’équation:

(4490)\[{\int}_{{\Omega}_{e}}{D}^{t}.E.(B.U-D.{\delta}^{\text{*}}).d\Omega =l.{\psi}^{'}({\delta}^{\text{*}})\]

Dérivons cette équation par rapport à \(U\) :

(4491)\[{\int}_{{\Omega}_{e}}{D}^{t}.E.(B.U-D.\frac{d{\delta}^{\text{*}}}{dU}).d\Omega =l.{\psi}^{\text{''}}({\delta}^{\text{*}}).\frac{d{\delta}^{\text{*}}}{dU}\]

Donc:

(4492)\[\frac{d{\delta}^{\text{*}}}{\mathrm{dU}}={\left[l.{\Psi}^{''}(\delta \text{*})+{\int}_{{\Omega}_{e}}{D}^{t}.E.D.d\Omega \right]}^{-1}.{\int}_{{\Omega}_{e}}{D}^{t}.E.B.d\Omega\]

Il suffit alors de calculer \({\Psi}^{''}\) , on distingue trois cas de figure:

1. Si \(\delta =0\) , alors \({\Psi}^{''}=0\) ;

2. Si l’élément est en régime linéaire, \({\Psi}^{''}={\Psi}_{\mathrm{lin}}^{''}=\frac{{\sigma}_{c}}{\kappa}.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\kappa ).{\mathrm{Id}}_{2\times 2}\) ;

3. Si l’élément est en régime dissipatif, \({\Psi}^{''}={\Psi}_{\mathrm{dis}}^{''}=C.{\mathrm{Id}}_{2\times 2}+\tilde{C}.\delta \otimes \delta\) avec

(4493)\[\begin{split}\lbrace \begin{array}{c}C=\frac{{\sigma}_{c}}{\parallel \delta \parallel }.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\parallel \delta \parallel )\\ \tilde{C}=-(\frac{1}{\parallel \delta \parallel }+\frac{{\sigma}_{c}}{{G}_{c}}).\frac{{\sigma}_{c}}{{\parallel \delta \parallel }^{2}}.\exp(-\frac{{\sigma}_{c}}{{G}_{c}}.\parallel \delta \parallel )\end{array}\end{split}\]

Références Bibliographiques

1. BarenblattG. I., The mathematical theory if equilibrium cracks in brittle fracture. Adv. Appl. Mech ., 7, pp. 55-129 (1962).

2. CharlotteM., FrancfortG.A. , MarigoJ.-J. and TruskinovskyL., Revisiting brittle fracture as an energy minimization problem : comparison of Griffith and Barenblatt surface energy models. Proceedings of the Symposium on « Continuous Damage and Fracture » The data science library, Elsevier, edited by A. Benallal , Paris , pp. 7-18, (2000).

3. Francfort G.A. and Marigo J.-J., Revisiting brittle fracture as an energy minimization problem. J. Mech. Phys. Solids , 46 (8), pp. 1319-1342 (1998).

4. JirasekM., Comparative study on elements with embedded discontinuities. Comp. Meth. Appl. Mech. Engng. , 188 pp. 307-330 (2000).

5. JirasekM., Embedded crack models for concrete fracture. In R. de Borst, H. N. Mang Bicanic, etG. Meschke, editors, Comp. Mod. of concrete Struct. (EURO-C) , pp. 291-300 (1998).

6. Laverne J. Formulation énergétique de la rupture par des modèles de forces cohésives, considérations théoriques et implantations numériques, Thèse de doctorat de l’Université Paris XIII , (Novembre 2004).

7. LorentzE. and BadelP., A new path-following constraint for strain-softening finite element simulations. Int. J. Num. Meth. Eng. 60 pp. 499-526 (2004).

Description des versions du document

Version Aster	Auteur(s) ou contributeur(s), organisme	Description des modifications
9.4	J.Laverne	Texte initial
10.2	M.Abbas	Mise en forme de la partie pilotage