1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2%% Fichier : ViscoCompressible 3%% Auteur : th202608@pleiades077.intra.cea.fr 4%% Date : 29 avril 2014 5%% Répertoire : /home/th202608/Documents/notes/2014/TutorielMFront/ 6%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7 8\section{Implantation d'une loi viscoplastique isotrope compressible} 9\label{sec:impl-dune-loi} 10 11Cette section décrit les étapes de l'implantation d'une loi de 12comportement élasto-viscoplastique typique du combustible nucléaire~: 13\begin{itemize} 14\item nous commençons par décrire la loi de comportement en la 15 décomposant en parties élémentaires (élasticité, viscoplasticité)~; 16\item nous détaillons ensuite comment implanter les propriétés 17 élastiques de la loi (module d'\nom{Young}, coefficient de 18 \nom{Poisson}, coefficient de dilatation thermique)~; 19\item une première implantation est alors proposée. Elle se base sur 20 une méthode de résolution par perturbation qui ne nécessite pas de 21 fournir le jacobien du système. L'implantation obtenue n'est pas 22 optimale mais permet de réaliser rapidement des tests~; 23\item nous montrons ensuite comment réaliser des premiers tests de 24 vérifications~; 25\item nous reprenons ensuite le travail d'implantation de la loi pour 26 obtenir une implantation optimisée de la loi~: 27 \begin{itemize} 28 \item par le calcul des termes du jacobien du système implicite~; 29 \item par le calcul de la matrice tangente cohérente. 30 \end{itemize} 31% \item nous montrons ensuite comment supporter les hypothèses de 32% contraintes planes. 33% \item nous concluons en montrant quelques résultats sur un calcul de 34% structure en grandes transformations~; 35\end{itemize} 36 37\subsection{Description de la loi de comportement} 38 39Nous considérons ici l'implantation d'une loi de comportement 40viscoplastique isotrope et compressible. Cette loi couple l'écoulement 41viscoplastique à l'évolution de la porosité du matériau. 42 43\subsubsection{Partition des déformations} 44 45On fait l'hypothèse classique en petites transformations que la 46déformation totale \(\tepsilonto\) se répartit additivement en une 47déformation élastique \(\tepsilonel\), une déformation viscoplastique 48\(\tepsilonvis\) et une déformation thermique \(\tepsilonth\)~: 49\begin{equation} 50 \label{eq:1} 51 \tepsilonto=\tepsilonel+\tepsilonvis+\tepsilonth 52\end{equation} 53 54Soulignons dès à présent que cette relation, d'apparence anodine, sera 55en réalité centrale dans l'implantation que nous allons proposer~: 56\begin{itemize} 57\item elle permet une grande modularité. Nous décrivons ailleurs qu'il 58 \og~suffit~\fg{} d'ajouter un nouveau terme de déformation pour 59 représenter un nouveau phénomène (écoulement viscoplastique sous 60 irradiation, plasticité dilatante représentant la décohésion des 61 joints de grains, fissuration macroscopique) sans que cela impacte 62 l'implantation de la partie viscoplastique 63 (voir~\cite{helfer14:_prise_mfron}). 64\item elle permet un calcul quasi-automatique de la matrice tangente 65 cohérente. Si la matrice tangente cohérente n'est pas utilisée par 66 le code éléments finis \castem{}, c'est un élément important pour 67 les performances de la loi dans l'application combustible 68 Cyrano~\cite{thouvenin_edf_2010}~; 69\item elle permet un support simple et non intrusif des hypothèses en 70 contraintes planes, notamment l'hypothèse d'axisymétrie en 71 contrainte plane généralisée utilisée par l'application combustible 72 Cyrano~\cite{thouvenin_edf_2010}. 73\end{itemize} 74 75Ces différents points justifient largement de conserver cette équation 76dans le système implicite que nous allons étudier dans la 77suite. L'implantation proposée pourra apparaître comme moins optimisée 78que ce que nous avons proposé par 79ailleurs~\cite{helfer_implantation_2005}. 80 81\subsubsection{Élasticité} 82 83Le comportement élastique relie les déformations élastiques 84\(\tepsilonel\) aux contraintes \(\tsigma\) par la loi de \nom{Hooke} 85isotrope~: 86\[ 87\tsigma=\tDq\,\tepsilonel{}=\lambda\,\trace{\tepsilonel}\,\tenseur{I}+2\,\mu\,\tepsilonel 88\] 89où nous avons utilisé les notations suivantes~: 90\begin{itemize} 91\item 92 \(\tDq\paren{T,p}=\lambda\,\tenseur{I}\,\otimes\,\tenseur{I}+2\,\mu\,\tenseurq{I}\) 93 est la matrice de raideur, supposée dépendre de la température \(T\) 94 et de la porosité~; 95\item \(\lambda\paren{T,p}\) et \(\mu\paren{T,p}\) sont respectivement les 96premier et second coefficients de \nom{Lamé} du matériau~; 97\item \(\trace{\tepsilonel}\) désigne la trace du tenseur \(\tepsilonel\), c'est 98 à dire la somme de ses termes diagonaux~; 99\item \(\tenseur{I}\) est le tenseur identité d'ordre \(2\)~; 100\item \(\tenseurq{I}\) est le tenseur identité d'ordre \(4\)~; 101\item \(\tenseur{I}\otimes\tenseur{I}\) est le tenseur d'ordre \(4\) 102 tel que pour tout tenseur \(\tenseur{s}\)~: 103\[ 104\tenseur{I}\otimes\tenseur{I}\,\tenseur{s} = \trace{s}\,\tenseur{I} 105\] 106\end{itemize} 107 108\paragraph{Propriétés élastiques} En pratique, les coefficients de 109\nom{Lamé} se déduisent du module d'\nom{Young} \(E\paren{T,p}\), 110fonction de la température et de la porosité \(p\) du matériau, et du 111coefficient de \nom{Poisson} \(\nu\), supposé constant~: 112\[ 113\lambda\paren{T,p}=\Frac{\nu\,E\paren{T,p}}{\paren{1+\nu}\paren{1-2\,\nu}} \quad\quad \mu\paren{T,p}=\Frac{E\paren{T,p}}{2\,\paren{1+\nu}} 114\] 115 116Pour illustrer notre propos, nous utiliserons une loi qui n'est pas 117celle des applications combustible \cea{}, mais une loi librement 118accessible~\cite{fink_thermodynamic_1981}~: 119\[ 120E\paren{T,p}=2,26\,10^{11}\,\paren{1-2,62\,p}\,\paren{1-1,31\,10^{-4}\,\paren{T-273,15}} 121\] 122Cette même référence préconise un coefficient de \nom{Poisson} 123constant d'une valeur de \(0,316\). 124 125\subsubsection{Dilatation thermique} 126 127Dans le cas des petites transformations, la dilatation thermique d'un 128corps est généralement traitée par le code appelant, en amont de la 129loi de comportement qui reçoit en entrée une déformation totale 130parfois appelée \og~mécanique~\fg. Nous avons choisi de traiter 131explicitement la dilatation thermique pour pouvoir~: 132\begin{itemize} 133\item traiter de manière exacte la dilatation thermique, la plupart 134 des codes réalisant une approximation que nous décrirons plus 135 loin~; 136\item utiliser notre loi en grandes déformations. 137\end{itemize} 138 139L'annexe~\ref{sec:rappels-sur-la} donne des précisions sur le 140traitement de la dilatation thermique. Dans la suite, le coefficient 141de dilatation thermique \(\alpha\) sera supposé constant et égal à 142\(10^{-5}\,K^{-1}\). 143 144\subsubsection{Écoulement viscoplastique} 145 146Nous supposons qu'il existe un potentiel de dissipation \(\Phi\) 147fonction des contraintes dont dérive l'écoulement 148viscoplastique~\cite{moreau_sur_1970,besson_mecanique_2001,chaboche_mecanique_2009}. 149 150La vitesse de déformation viscoplastique \(\tdepsilonvis\) s'exprime 151alors~: 152\[ 153\tepsilonvis = \deriv{\Phi}{\tsigma} 154\] 155 156\paragraph{Contrainte équivalente} Dans la suite, nous faisons 157l'hypothèse que le potentiel de dissipation \(\Phi\) est en fait une 158fonction d'une contrainte équivalente (scalaire) \(s_{eq}\) telle 159que~: 160\[ 161s_{eq}=\sqrt{A\paren{f}\,p^{2}+B\paren{f}\,\sigmaeq^{2}} 162\] 163où \(p\) est la pression hydrostatique et \(\sigmaeq\) est la 164contrainte équivalente au sens de \textsc{Von Mises}~: 165\[ 166\begin{aligned} 167 p = \Frac{1}{3}\trace{\sigma}\\ 168 \sigmaeq = \sqrt{\Frac{3}{2}\mathop{dev}\paren{\tenseur{s}}\,\colon\,\mathop{dev}\paren{\tenseur{s}}} 169\end{aligned} 170\] 171où \(\mathop{dev}\paren{\tenseur{s}}\) est la partie déviatorique des contraintes~: 172\[ 173\mathop{dev}\paren{\tenseur{s}}=\sigma-\Frac{1}{3}\,\trace{\sigma}\tenseur{I} 174\] 175Par définition, \(\mathop{dev}\paren{\tenseur{s}}\) est de trace nulle. 176 177Les fonctions \(A\paren{f}\) et \(B\paren{f}\) dépendent de la porosité. Un 178choix possible pour ces fonctions est~: 179\[ 180 \begin{aligned} 181 A\paren{f} &= \Frac{9}{4}\paren{E\,\paren{f^{-1/E}-1}}^{-\frac{2\,E}{E+1}}\\ 182 B\paren{f} &= \paren{1+\frac{2}{3}\,f}\,\paren{1-f}^{-\frac{2\,E}{E+1}} 183 \end{aligned} 184\] 185où \(E\) est une constante. Nous verrons dans la suite qu'il est 186intéressant d'introduire un coefficient \(C_{A}\) tel que~: 187\[ 188A\paren{f} = \Frac{9\,C_{A}}{4}\paren{E\,\paren{f^{1/E}-1}}^{-\frac{2\,E}{E+1}} 189\] 190Si \(C_{A}\) est nul, on retrouve une loi viscoplastique isotrope 191incompressible (sans effet de la porosité), ce qui est utile lors de 192l'étape de vérification de la loi. 193 194Cette forme de potentiel de dissipation provient de l'homogénisation 195du problème d'une matrice poreuse incompressible isotrope supposée 196obéir à une loi viscoplastique de 197\nom{Norton}~\cite{Michel1992783}. Elle a été utilisée à plusieurs 198reprises pour décrire le combustible 199nucléaire~\cite{monerie_overall_2006}. La constante \(E\) apparaissant 200dans la définition des fonctions \(A\paren{f}\) et \(B\paren{f}\) est 201l'exposant de la loi de de \nom{Norton} de la matrice. 202 203\paragraph{Potentiel elliptique} Les mêmes résultats d'homogénéisation 204conduisent à proposer comme potentiel de dissipation un potentiel dit 205\og~elliptique~\fg{} qui généralise l'expression d'une loi de 206\nom{Norton}~: 207\[ 208\Phi\paren{s_{eq}}=\Frac{A_{\Phi}}{E+1}\,s_{eq}^{E+1} 209\] 210 211\paragraph{Expression usuelle de l'écoulement plastique} Développons 212maintenant l'expression de la vitesse d'écoulement viscoplastique~: 213\[ 214\tepsilonvis = \deriv{\phi}{\tsigma} = \deriv{\phi}{s_{eq}}\deriv{s_{eq}}{\tsigma} 215= \deriv{\phi}{s_{eq}}\left[\deriv{s_{eq}}{p}\deriv{p}{\tsigma}+\deriv{s_{eq}}{\sigmaeq}\deriv{\sigmaeq}{\tsigma}\right] 216= \deriv{\phi}{s_{eq}}\left[\Frac{1}{3}\deriv{s_{eq}}{p}\tenseur{I}+\deriv{s_{eq}}{\sigmaeq}\tenseur{n}\right] 217\] 218où nous avons introduit le classique tenseur normal \(\tenseur{n}\) 219égal à la dérivée de la contrainte équivalente de \nom{Von Mises} par 220rapport au tenseur des contraintes~: 221\[ 222\tenseur{n}=\deriv{\sigmaeq}{\tsigma}=\Frac{3}{2\,\sigmaeq}\mathop{dev}\paren{\tenseur{s}} 223\] 224où \(s\) est le déviateur des contraintes introduit plus haut. Le 225tenseur \(n\) est donc également déviatorique (de trace nulle). 226 227Nous pouvons alors donner les expressions des différentes dérivées 228apparaissant dans cette expression~: 229\[ 230\deriv{\phi}{s_{eq}} = s_{eq}^{E}\quad\quad 231\deriv{s_{eq}}{p} = \Frac{1}{s}\,A\paren{f}\,p \quad\quad 232\deriv{s_{eq}}{\sigmaeq} = \Frac{1}{s}\,B\paren{f}\,\sigmaeq 233\] 234 235\subsubsection{Évolution de la porosité} 236 237L'évolution de la porosité est donnée par~: 238\[ 239\dot{f}=\paren{1-f}\,\trace{\tdepsilonvis} 240\] 241 242\subsection{Propriétés matériau} 243 244Telle que nous l'avons présentée, notre loi de comportement s'appuie 245sur différentes propriétés du matériau~: 246\begin{itemize} 247\item le module d'\nom{Young}~; 248\item le coefficient de \nom{Poisson}~; 249\item le coefficient de dilatation thermique. 250\end{itemize} 251 252Afin de traiter la question de leur implantation en \mfront{}, prenons 253le temps d'examiner deux questions~: 254\begin{itemize} 255\item ces propriétés doivent-elles être évaluées par la loi de 256 comportement ou par le code appelant ? 257\item vaut-il mieux utiliser des coefficients matériaux ou des 258 paramètres ? 259\end{itemize} 260 261\subsubsection{Évaluation par la loi ou par le code appelant ?} 262 263Ces propriétés doivent-elles être implantées dans la loi de 264comportement ou implantées à part, par exemple pour être accessibles 265dans d'autres lois de comportement ou, dans le cas du code aux 266éléments finis \castem{}, pour permettre la construction d'une matrice 267de raideur ? 268 269La première solution nous apparaît comme la plus satisfaisante sur le 270principe~: l'identification de la loi se fait pour des propriétés 271élastiques données qui font donc partie intégrante de la loi. En 272pratique, la seconde solution a pour l'instant toujours été retenue 273car elle minimisait le travail d'implantation d'une propriété. 274 275Dans le cas du code \castem{}, cette propriété était souvent écrite en 276langage \gibiane{} dans le jeu de données de la simulation, ce qui 277n'était pas sans poser de nombreux problèmes~: 278\begin{itemize} 279\item l'écriture du code était fastidieuse pour les propriétés 280 complexes et d'une mise en \oe{}uvre délicate dès lors qu'une 281 propriété dépendait de plusieurs paramètres~; 282\item le code était \og~copier/coller~\fg{} d'un jeu de données à un 283 autre. Des tentatives de mutualisation ont bien été menées en se 284 basant sur un fichier \og~UTILPROC~\fg{} commun, mais force est de 285 constater que malgré les efforts déployés, ces tentatives ne se sont 286 pas révélées pérennes~: 287\item l'implantation ne peut être réutilisée dans d'autres 288 solveurs. Or, nous verrons dans la suite qu'il est intéressant et 289 pratique de réaliser des tests unitaires sur la base du logiciel 290 \mtest{} puis de réaliser les calculs de structure à l'aide de 291 \castem{}. 292\end{itemize} 293 294\mfront{} apporte une solution quelque peu différente qui concilie les 295deux solutions évoquées plus haut~: un même fichier source peut servir 296à générer une bibliothèque appelée depuis \castem{} (ou \aster{}) ou 297directement être intégrée à l'implantation de la loi. 298 299\begin{figure}[htbp] 300 \centering 301 \lstinputlisting{@top_srcdir@/docs/tutorial/mfront/UO2_PoissonRatio_Fink1981.mfront} 302 \caption{Implantation du coefficient de \nom{Poisson} du dioxide d'uranium.} 303 \label{fig:tutorial:uo2:Poisson} 304\end{figure} 305 306\subsubsection{Implantation du coefficient de \nom{Poisson}} La 307figure~\ref{fig:tutorial:uo2:Poisson} présente l'implantation du 308coefficient de Poisson de la loi. Nous renvoyons à la documentation 309générale de \mfront{} pour les détails liés à l'implantation d'une 310propriété matériau~\cite{helfer_generateur_2013-1}. 311 312Il peut sembler lourd d'écrire un fichier source pour une propriété 313constante, à la fois en terme d'implantation, mais aussi en terme 314d'efficacité numérique (les codes de calcul savent gérer des valeurs 315constantes !). Cela présente cependant plusieurs avantages~: 316\begin{enumerate} 317\item il est possible de stocker ce fichier dans une base de données, 318 la base de données \sirius{} de la plate-forme \pleiades{} 319 notamment~; 320\item cette valeur n'est pas dans les jeux de données des codes 321 appelant, mais conservée à un endroit unique~: on ne duplique pas 322 d'information physique à chaque nouveau calcul~; 323\item certains codes permettent d'optimiser l'évaluation d'une 324 fonction sans argument (le code \licos{} 325 notamment~\cite{helfer_licos_2011}) et l'appel ne se fait qu'une 326 fois. 327\end{enumerate} 328 329\paragraph{Choix du nom du fichier \mfront{}} Le nom de la fonction 330générée sera \texttt{UO2\_\-Poisson\-Ratio\_\-Fink\-1981} (combinaison 331du nom du matériau donné en ligne \(2\) et du nom de la loi donnée en 332ligne \(3\)). Les noms des sources \cxx{} générées commenceront par ce 333nom de fonction. Pour simplifier l'intégration de ce fichier \mfront{} 334dans un environnement de compilation, il est recommandé de nommer le 335fichier source \mfront{} 336\og~\texttt{UO2\_\-Poisson\-Ratio\_\-Fink\-1981\-.mfront}~\fg{} car 337cela permet d'utiliser des règles de compilation génériques. 338 339\begin{figure}[htbp] 340 \centering 341 \lstinputlisting{@top_srcdir@/docs/tutorial/mfront/UO2_YoungModulus_Fink1981.mfront} 342 \caption{Implantation du module d'\nom{Young} du dioxide d'uranium.} 343 \label{fig:tutorial:uo2:YoungModulus} 344\end{figure} 345 346\subsubsection{Implantation du module d'\nom{Young}} La 347figure~\ref{fig:tutorial:uo2:YoungModulus} montre comment est implanté 348le module d'\nom{Young} du dioxide d'uranium. Notons quelques 349nouveautés~: 350\begin{itemize} 351\item la définition de bornes physiques de la loi en lignes \(21\) et 352 \(22\) (mot clé {\tt @Physical\-Bounds}). En dehors de ces bornes, 353 les arguments ont des valeurs non physiques~; 354\item la définition de bornes de la loi en ligne \(24\). Cette borne 355 indique la limite de validité de la corrélation (mot clé {\tt 356 @Bounds})~; 357\item la définition de noms de glossaire pour les variables d'entrée 358 en lignes \(16\) et \(17\). Ces noms de glossaire permettent à des 359 codes comme \mtest{} et \licos{} d'évaluer la propriété. 360\end{itemize} 361 362\begin{figure}[htbp] 363 \centering 364 \lstinputlisting{@top_srcdir@/docs/tutorial/mfront/UO2_ThermalExpansion_MATPRO.mfront} 365 \caption{Implantation du coefficient de dilatation thermique du 366 dioxide d'uranium.} 367 \label{fig:tutorial:uo2:thermalexpansion} 368\end{figure} 369 370\paragraph{Implantation du coefficient de dilatation thermique} 371L'implantation du coefficient de dilatation thermique est similaire à 372celle du coefficient de Poisson. Elle est indiquée en 373figure~\ref{fig:tutorial:uo2:thermalexpansion}. Notons que nous avons 374explicitement indiqué la température de référence \(T_{\alpha}\) dans 375une constante nommée {\tt Reference\-Temperature}. Cette constante 376sera utilisée par la loi de comportement pour le calcul de la 377dilatation thermique. 378 379\paragraph{Compilation} 380Pour pouvoir être utilisé avec \mtest{}, l'outil de test unitaire 381livré avec \mfront{}, les propriétés matériau précédentes doivent être 382compilées avec l'interface \castem{}~: 383 384\begin{flushleft} 385 {\tt 386 \$ \small mfront -{}-obuild -{}-interface=castem UO2\_PoissonRatio\_Fink1981.mfront UO2\_YoungModulus\_Fink1981.mfront UO2\_ThermalExpansion\_MATPRO.mfront\\ 387 Treating target : all\\ 388 The following library has been build :\\ 389 - libCastemUO2.so 390 } 391\end{flushleft} 392 393Cette commande crée une librairie nommée {\tt libUO2Castem.so} dans un 394sous-répertoire nommé {\tt src}. 395 396Dans une phase de vérification, d'autres interfaces peuvent être 397intéressantes. Par l'exemple, l'interface \python{} crée un module qui 398peut être utilisé pour la vérification. 399 400\begin{figure}[htbp] 401 \centering 402 \lstinputlisting[language=gnuplot]{@top_srcdir@/docs/tutorial/mfront/UO2.tplot} 403 \caption{Script \tplot{} permettant de tracer l'évolution du module 404 d'\nom{Young} de l'uranium au travers de son implantation 405 \mfront{} et de la comparer à une expression analytique.} 406 \label{fig:tutorial:uo2:youngmodulus:tplot} 407\end{figure} 408 409\begin{figure}[htbp] 410 \centering 411 \includegraphics[width=10cm]{@top_srcdir@/docs/tutorial/images/UO2_YoungModulus_Fink1981.png} 412 \caption{Évolution du module d'\nom{Young} de l'$UO_{2}$. Données 413 libres~\cite{fink_thermodynamic_1981}} 414 \label{fig:tutorial:uo2:youngmodulus:check} 415\end{figure} 416 417\paragraph{Vérification} Il existe différentes façons de vérifier 418l'implantation d'une loi \mfront{}. \licos{} propose un utilitaire 419graphique nommée \tplot{} qui propose un langage de script proche de 420celui du classique \gnuplot{} et qui permet de comparer les valeurs 421retournées par la fonction générée avec \mfront{} à l'expression 422analytique. La figure~\ref{fig:tutorial:uo2:youngmodulus:tplot} 423présente le script utilisé pour tracer la 424figure~\ref{fig:tutorial:uo2:youngmodulus:check}~: la seule 425spécificité par rapport à \gnuplot{} est la ligne \(4\) (mot clé {\tt 426 import}) qui permet d'utiliser la fonction générée par \mfront{}. 427 428Pour une vérification systématique, les applications de la plate-forme 429\pleiades{} peuvent utiliser l'utilitaire \mpplot{} qui permet de 430sortir des fichiers colonnes qui peuvent être comparés à des fichiers 431de référence à l'aide de l'utilitaire \pcheck{}. 432 433\subsubsection{Les coefficients de la loi~: propriété matériau ou paramètre ?} 434 435En plus des propriétés thermoélastiques, la loi de comportement qui 436nous intéresse dans cette section présente de nombreux coefficients~: 437\begin{itemize} 438\item exposant \(E\) de la loi de \nom{Norton}~; 439\item intensité \(A_{\phi}\) de l'écoulement~; 440\item le coefficient \(C_{A}\) de la loi de comportement. 441\end{itemize} 442Par rapport aux propriétés thermoélastiques traitées au paragraphes 443précédents, ces coefficients sont propres à la loi et ont peu de 444chances d'être réutilisés dans d'autres lois. 445 446Ces coefficients peuvent être déclarés soit comme des propriétés 447matériau {\em que le code appelant doit fournir} soit comme des {\em 448 paramètres}. Chacune de ces modalités présente des avantages et des 449inconvénients que nous allons maintenant détailler~: le fait de 450déclarer un coefficient comme étant un paramètre ou une propriété 451matériau est une affaire de compromis et d'habitude de travail. 452 453\paragraph{Propriétés matériau} 454Afin de pouvoir modifier les valeurs des coefficients d'une loi depuis 455un jeu de données \castem{} ou \aster{}, l'interface {\tt castem} 456\og~standard~\fg{} impose de déclarer ces coefficients comme des 457propriétés matériau (mot clé {\tt @Material\-Property}) que doit 458fournir le code appelant. 459 460Le principal avantage de cette méthode est de permettre l'écriture de 461lois de comportement génériques. En effet, l'écriture d'une loi de 462comportement avec un langage de bas niveau ({\tt fortran} ou {\tt C} 463par exemple) est une tâche longue, difficile et fastidieuse. L'usage 464de \mfront{} limite considérablement la portée de cet argument. 465 466Le second avantage de cette méthode est que l'utilisateur a une grande 467liberté et peut faire des tests en dehors de la loi de comportement, 468par exemple faire varier \(A_{\phi}\) avec le taux de combustion. 469 470Or, cet avantage est à notre avis le principal défaut de cette méthode 471puisqu'une partie de l'information physique est dissociée de 472l'implantation de la loi~: nous considérons préférable de rassembler 473toutes les informations physiques dans un unique fichier \mfront{} (ou 474éventuellement dans des fichiers dédiés comme pour le module 475d'\nom{Young}) et cela d'autant plus que l'on souhaite capitaliser les 476lois dans des bases de données (la base \sirius{} de la plate-forme 477\pleiades{} par exemple). 478 479L'autre défaut de cette méthode est que la gestion des propriétés 480matériau par le code appelant peut être coûteuse~: une implantation 481d'une loi utilisant de nombreuses propriétés matériau est 482significativement plus lente à l'exécution que l'implantation de la 483même loi utilisant des propriétés matériau \og~en dur~\fg{}. 484 485\paragraph{La notion de paramètres} Le besoin de modifier à 486l'exécution des coefficients de la loi demeure cependant, par exemple 487dans le cadre d'études de sensibilité. Dans ce cas, \mfront{} propose 488une alternative efficace aux propriétés matériau~: la notion de 489paramètre (mot clé {\tt @Parameter}). Un paramètre se comporte comme 490une variable globale dont la valeur par défaut est spécifiée dans le 491fichier \mfront{} (ce qui est impossible avec une propriété 492matériau). La valeur d'un paramètre peut être modifiée à l'exécution 493sans passer par l'appel standard à la loi de comportement et donc sans 494induire de coût supplémentaire, si le code appelant le permet. Seul le 495code \castem{} ne le permet pas, mais les applications de la 496plate-forme \pleiades{}, le code \aster{} ou \mtest{}, peuvent 497modifier les paramètres de la loi depuis le langage \python{} ou le 498\cxx{}. 499 500Insistons sur le fait qu'un paramètre est alors une caractéristique 501propre à une loi et non à un matériau~: si deux matériaux distincts 502utilisent la même loi, la modification d'un paramètre affecte le 503comportement des deux matériaux. 504 505Dans la suite, nous déclarerons les coefficients \(E\), \(A_{\phi}\) 506et \(C_{A}\) comme des paramètres. 507 508\subsection{Discrétisation et implantation par une $\theta$ méthode} 509 510Nous avons choisi d'intégrer notre loi de comportement par une 511$\theta$-méthode pour différentes raisons~: 512\begin{itemize} 513\item c'est la méthode la plus performante~; 514\item on a accès à la matrice tangente cohérente~; 515\item les autres lois du combustible (modèle de fissuration DDIF2 516 notamment) sont déjà implantées dans ce formalisme. 517\end{itemize} 518 519\subsubsection{Gestion de la dilatation thermique} 520 521Nous utilisons le mot clé {\tt @Compute\-Thermal\-Expansion} pour que 522\mfront{} gère la dilatation thermique. Celle-ci sera gérée 523différemment selon que la loi est utilisée en petites déformations ou 524en grandes déformations. Les variables {\tt eto} et {\tt deto} qui 525apparaîtront dans l'implantation sont respectivement la déformation 526totale mécanique (déformation totale moins la dilatation 527thermique) en début de pas et son incrément. 528 529\subsubsection{Réduction du nombre de variables internes} 530 531Telle que présentée au paragraphe précédent, notre loi possède \(3\) 532variables internes~: 533\begin{itemize} 534\item la déformation élastique \(\tepsilonel\)~; 535\item la déformation viscoplastique \(\tepsilonvis\)~; 536\item la porosité~; 537\end{itemize} 538 539\paragraph{Cas de l'écoulement viscoplastique} 540Commençons par traiter le cas de la déformation viscoplastique 541\(\tepsilonvis\). Il serait assez simple de retirer complètement cette 542variable du système et économiser ainsi \(6\) variables internes dans 543le système. Pour simplifier l'écriture du système implicite et pour 544simplifier les post-traitements, nous avons choisi dans ce tutoriel un 545choix moins optimal qui introduit deux variables internes 546supplémentaires. L'origine de ces variables provient du fait qu'il est 547possible de décomposer la vitesse d'écoulement viscoplastique suivant 548deux directions orthogonales dans l'espace des contraintes: notons 549respectivement \(\dot{p}_{v}\) et \(\dot{p}_{d}\) la trace de 550\(\tepsilonvis\) et la projection de \(\tepsilonvis\) suivant 551\(\tenseur{n}\). Nous avons~: 552\[ 553\begin{aligned} 554 \dot{p}_{v} &= \deriv{\phi}{s_{eq}}\deriv{s_{eq}}{p}\\ 555 \dot{p}_{d} &= \deriv{\phi}{s_{eq}}\deriv{s_{eq}}{\sigmaeq}\\ 556\end{aligned} 557\] 558 559Ces équations introduisent les deux variables internes recherchées~: 560\begin{itemize} 561\item la variable \(p_{v}\) qui représente le changement de volume dû 562 à l'écoulement viscoplastique (liés à la croissance ou la diminution 563 de la porosité)~; 564\item la variables \(p_{d}\) qui représente la déformation 565 viscoplastique déviatorique cumulée. 566\end{itemize} 567 568\paragraph{Traitement de la porosité} En utilisant la variable 569\(p_{v}\), l'équation d'évolution de la porosité est donnée par~: 570\[ 571\dot{f}=\paren{1-f}\dot{p}_{v} 572\] 573Cette relation explicite entre \(p_{v}\) et \(f\) va nous permettre 574d'exclure la porosité du système implicite. Dans \mfront{}, ce type de 575variable interne est qualifiée d'auxiliaire. 576 577\subsubsection{Écriture du système implicite} 578 579Connaissant les valeurs des variables internes à l'instant \(t\), nous 580cherchons à connaître leurs valeurs en fin de pas de temps, à 581l'instant \(t+\Delta\,t\). 582 583\paragraph{Quelques notations} 584Dans la suite, pour chacune des variables internes \(x\), nous 585noterons~: 586\begin{enumerate} 587\item \(\Delta x\) son incrément sur le pas de temps~; 588\item \(\bts{x}\) sa valeur en début de pas de temps~; 589\item \(\mts{x}=\bts{x}+\theta\,\Delta\,x\) l'estimation de sa valeur 590 au temps intermédiaire \(t+\theta\,\Delta\,t\). Nous parlerons 591 souvent, de manière abusive, de la valeur en milieu de pas de temps 592 de la variable \(x\) (en pratique, \(\theta\) est souvent pris égal 593 à \(0,5\), sauf pour le traitement des lois plastiques)~; 594\item \(\ets{x}=\bts{x}+\Delta\,x\) sa valeur en fin de pas de temps~; 595\end{enumerate} 596 597\paragraph{Équation relative à la déformation élastique} La partition 598des déformations s'écrit, sous forme incrémentale~: 599\[ 600\Delta\tepsilonel+\Delta\tepsilonvis-\Delta\tepsilonto=0 601\] 602En introduisant les incréments des variables \(p_{v}\) et \(p_{d}\), 603nous obtenons l'équation \(f_{\tepsilonel}\) associée à la déformation 604élastique~: 605\[ 606f_{\tepsilonel}=\Delta\tepsilonel+\Frac{1}{3}\,\Delta\,p_{v}\,\tenseur{I}+\Delta\,p_{d}\mts{\tenseur{n}}-\Delta\tepsilonto 607\] 608où \(\mts{\tenseur{n}}\) est l'évaluation du tenseur normal à 609l'instant \(t+\theta\Delta\,t\). 610 611\paragraph{Équation relative au changement de volume viscoplastique} 612L'équation \(f_{p_{v}}\) associée à la variable \(p_{v}\) est~: 613\[ 614f_{p_{v}}=\Delta\,p_{V}-\Delta\,t\,\mts{\deriv{\phi}{s}}\,\mts{\deriv{s}{p}} 615\] 616 617\paragraph{Équation relative à la déformation viscoplastique cumulée} 618L'équation \(f_{p_{d}}\) associée à la variable \(p_{d}\) est~: 619\[ 620f_{p_{d}}=\Delta\,p_{d}-\Delta\,t\,\mts{\deriv{\phi}{s}}\,\mts{\deriv{s}{\sigmaeq}} 621\] 622 623\paragraph{Mise à jour de la porosité} La porosité est mise en jour au 624cours des itérations de \nom{Newton} par la relation~: 625\[ 626\mts{f}=\bts{f}+\Frac{\paren{1-f}\,\theta\,\Delta\,p_{v}}{1+\theta\,\Delta\,p_{v}} 627\] 628 629\subsubsection{Implantation} 630 631Les équations définissant \(f_{\tepsilonel}\) \(f_{p_{v}}\) et 632\(f_{p_{d}}\) forment un système de trois équations pour trois 633inconnues~: elles sont donc suffisantes pour déterminer les incréments 634des variables internes sur un pas de temps. Reste donc à réaliser 635l'implantation de la loi. Nous commençons par utiliser un algorithme 636de \nom{Newton-Raphson} dont la jacobienne est calculée numériquement. 637 638\begin{figure}[htbp] 639 \centering 640 \lstinputlisting[firstline=1,lastline=41]{@top_srcdir@/docs/tutorial/mfront/EllipticCreep-1.mfront} 641 \caption{Première implantation de la loi d'écoulement viscoplastique 642 isotrope incompressible, première partie~: déclaration des 643 variables.} 644 \label{fig:tutorial:uo2:EllipticCreep1-1} 645\end{figure} 646 647\paragraph{La partie déclarative} Un fichier \mfront{} commence 648généralement par une partie déclarative reprise en 649figure~\ref{fig:tutorial:uo2:EllipticCreep1-1}~: 650\begin{itemize} 651\item la ligne \(1\) précise que nous allons utiliser une méthode 652 d'intégration implicite~; 653\item la ligne \(2\) donne le nom de la loi~; 654\item les lignes \(3\) et \(4\) renseignent les noms des auteurs et 655 la date~; 656\item les lignes \(6\) à \(7\) demandent l'inclusion d'un fichier 657 d'entête de la bibliothèque \TFEL{} permettant de calculer les 658 coefficients de \nom{Lamé} et la matrice d'élasticité à partir de 659 ces coefficients~; 660\item la ligne \(10\) précise que l'on utilise un algorithme de 661 \nom{Newton} avec jacobienne numérique~; 662\item les lignes \(12\) à \(20\) déclarent les variables internes et 663 leur associent des noms pour l'extérieur. Ces noms de glossaire sont 664 importants pour une adhérence aux applications de la plate-forme 665 \pleiades{} et pour \mtest{}~; 666\item les lignes \(22\) à \(27\) définissent les paramètres de la loi 667 et les noms à utiliser pour les modifier depuis \mtest{}, \aster{} 668 ou les applications de la plate-forme \pleiades{}~; 669\item les lignes \(26\) à \(36\) définissent quelques variables 670 locales qui seront utilisables dans l'ensemble des blocs de 671 code. Notons que pour traiter l'évolution de la porosité, nous 672 introduisons une variable interne représentant sa valeur au cours du 673 pas~: \mfront{} ne définit pas d'incrément pour les variables 674 auxiliaires~; 675\item les lignes \(38\) et \(39\) indiquent dans quels fichiers sont 676 implantés le module d'\nom{Young} et le coefficient de 677 \nom{Poisson}~; 678\item la ligne \(41\) demande le calcul de la dilatation thermique. 679\end{itemize} 680 681\paragraph{Remarque sur le type des variables} Certaines variables 682sont déclarées de type {\tt stress} et d'autres de type {\tt real}. Il 683s'agit dans les deux cas d'un alias vers un type de nombre réel 684(double précision dans toutes les interfaces considérées jusqu'à 685présent). Le choix d'utiliser le type {\tt stress} est uniquement 686informatif et vise à rendre le code plus explicite. 687 688\begin{figure}[htbp] 689 \centering 690 \lstinputlisting[firstline=43,lastline=48,firstnumber=43]{@top_srcdir@/docs/tutorial/mfront/EllipticCreep-1.mfront} 691 \caption{Première implantation de la loi d'écoulement viscoplastique 692 isotrope incompressible, seconde partie~: initialisation des 693 variables locales.} 694 \label{fig:tutorial:uo2:EllipticCreep1-2} 695\end{figure} 696 697\paragraph{Initialisation des variables locales} La 698figure~\ref{fig:tutorial:uo2:EllipticCreep1-2} montre comment 699initialiser les variables locales de la loi. Le coefficient de 700\nom{Poisson} étant constant, nous l'évaluons un fois pour l'ensemble 701de l'intégration. Nous y sauvegardons également la valeur de la 702porosité en début de pas. 703 704\begin{figure}[htbp] 705 \centering 706 \lstinputlisting[firstline=50,lastline=71,firstnumber=50]{@top_srcdir@/docs/tutorial/mfront/EllipticCreep-1.mfront} 707 \caption{Première implantation de la loi d'écoulement viscoplastique 708 isotrope incompressible, troisième partie~: calcul des contraintes 709 en milieu de pas et en fin de pas.} 710 \label{fig:tutorial:uo2:EllipticCreep1-3} 711\end{figure} 712 713\paragraph{Évaluation des contraintes en milieu et fin de pas de 714 temps} La figure~\ref{fig:tutorial:uo2:EllipticCreep1-3} montre 715comment sont évaluées les contraintes en milieu de pas et en fin de 716pas. Nous tombons ici sur une spécificité de l'implantation proposée. 717 718En effet, la plupart des lois peuvent se contenter de ne définir que 719le bloc {\tt @Compute\-Stress}~: \mfront{} évalue dans ce cas les 720différentes variables à l'instant considéré (milieu de pas ou fin de 721pas). 722 723\begin{figure}[htbp] 724 \centering 725 \includegraphics{@top_srcdir@/docs/tutorial/images/ImplicitDSL.pdf} 726 \caption{Ordre d'évaluation des différents blocs de code fournis par 727 l'utilisateur dans le cas d'une intégration implicite.} 728 \label{fig:ImplicitOrder} 729\end{figure} 730 731Ce n'est pas possible ici du fait de l'utilisation de la porosité 732comme variable auxiliaire~: celle-ci doit être mise à jour 733explicitement avant l'évaluation des propriétés élastiques. Pour bien 734comprendre ce point, on peut se référer à la 735figure~\ref{fig:ImplicitOrder} qui montre l'ordre d'appel des 736différents blocs de code définis par l'utilisateur. 737 738 739Le bloc de code {\tt @Compute\-Stress} sera utilisé pour évaluer les 740contraintes en milieu de pas et le bloc de code {\tt 741 @Compute\-Final\-Stress} en fin de pas de temps. 742 743\begin{figure}[!htbp] 744 \centering 745 \lstinputlisting[firstline=304,lastline=352,firstnumber=304]{@top_srcdir@/docs/tutorial/mfront/EllipticCreep.hxx} 746 \caption{Première implantation de la loi d'écoulement viscoplastique 747 isotrope incompressible, code généré par \mfront{} pour 748 l'évaluation des contraintes en milieu et en fin de pas.} 749 \label{fig:tutorial:uo2:EllipticCreep1-3b} 750\end{figure} 751 752Pour se convaincre que les contraintes sont évaluées correctement, il 753peut être intéressant de regarder quel est le code généré par 754\mfront{} pour ces deux blocs 755(figure~\ref{fig:tutorial:uo2:EllipticCreep1-3b}). 756 757\begin{figure}[!htbp] 758 \centering 759 \lstinputlisting[firstline=74,firstnumber=74]{@top_srcdir@/docs/tutorial/mfront/EllipticCreep-1.mfront} 760 \caption{Première implantation de la loi d'écoulement viscoplastique 761 isotrope incompressible, quatrième partie~: système implicite.} 762 \label{fig:tutorial:uo2:EllipticCreep1-4} 763\end{figure} 764 765La figure~\ref{fig:tutorial:uo2:EllipticCreep1-4} constitue le 766c\oe{}ur de la résolution. 767 768Le code contenu dans le bloc \texttt{@Integrator} est appelé après 769l'évaluation des contraintes en milieu de pas. La variable \(f\) 770contient l'estimation de la porosité en milieu de pas de temps et la 771variable \(sig\) l'estimation du tenseur des contraintes. 772 773Nous commençons donc par évaluer la pression \(p\), la contrainte 774équivalente au sens de \nom{Von Mises} \(\sigmaeq\) et la contrainte 775équivalente \(s\) en milieu de pas à l'aide de la contrainte en milieu 776de pas. 777 778Dans le cas où les contraintes ne sont pas trop faibles, nous 779définissons les équations associées aux variables viscoplastiques 780\(p_{v}\) et \(p_{d}\). Nous utilisons ici les conventions implicites 781de \mfront{}~: les variables \(f_{x}\) sont initialisées avant chaque 782appel à \(\Delta\,x\). Ainsi, si les contraintes sont trop faibles, 783les équations associées aux variables \(p_{v}\) et \(p_{d}\) seront~: 784\[ 785\begin{aligned} 786 f_{p_{v}}&=\Delta\,p_{v}&=0 \\ 787 f_{p_{d}}&=\Delta\,p_{d}&=0 \\ 788\end{aligned} 789\] 790qui se résolvent trivialement en~: 791\[ 792\Delta\,p_{v}=\Delta\,p_{d}=0 793\] 794 795\subsection{Premiers tests} 796 797\subsubsection{Compilation} 798 799La première étape des tests est la compilation de la loi. Nous 800utilisons l'interface {\tt castem} utilisée par le code \castem{}~: 801\begin{flushleft} 802 {\tt 803 \$ \small mfront -{}-obuild -{}-interface=castem EllipticCreep.mfront \\ 804 Treating target : all \\ 805 UO2\_ThermalExpansion-mfront.cxx: In function ‘double mfront::UO2\_ThermalExpansion()’: \\ 806 UO2\_ThermalExpansion-mfront.cxx:21:22: warning: unused variable ‘ReferenceTemperature’ [-Wunused-variable] \\ 807 The following libraries have been build : \\ 808 - libCastemUO2.so \\ 809 - libMFrontMaterialLaw.so \\ 810 - libCastemBehaviour.so 811 } 812\end{flushleft} 813Le message d'avertissement relatif à la variable {\tt 814 Reference\-Temperature} est sans importance. 815 816\paragraph{Analyse des messages d'erreurs} La phase de compilation 817peut montrer des problèmes de syntaxe. Par exemple, remplaçons la 818ligne \(59\) de la figure~\ref{fig:tutorial:uo2:EllipticCreep1-3} par 819celle-ci, contenant une faute de frappe ({\tt lamdba} au lieu de {\tt 820 lambda})~: 821\begin{center} 822 \lstinputlisting[firstline=59,lastline=59,firstnumber=59]{@top_srcdir@/docs/tutorial/mfront/EllipticCreep-1.mfront} 823\end{center} 824Dans ce cas, le compilateur GNU 825g++~\cite{fsf_gnu_2014} affichera le message 826d'erreur suivant~: 827\begin{flushleft} 828 {\tt 829 In file included from EllipticCreep.cxx:13:0: 830 EllipticCreep-1.mfront: In member function ‘void tfel::material::EllipticCreep<hypothesis, Type, false>::computeStress()’: 831 EllipticCreep-1.mfront:59:13: error: ‘lamdba’ was not declared in this scope 832 compilation terminated due to -Wfatal-errors. 833 make: *** [EllipticCreep.o] Erreur 1 834 terminate called after throwing an instance of 'std::runtime\_error' 835 what(): MFront::buildLibraries : libraries building went wrong 836 Abandon 837 } 838\end{flushleft} 839 840La ligne importante ici est la troisième~: 841\begin{flushleft} 842 {\tt 843 EllipticCreep-1.mfront:59:13: error: ‘lamdba’ was not declared in this scope 844 } 845\end{flushleft} 846qui indique qu'en ligne \(59\) à la colonne \(13\) la variable {\tt 847 lamdba} n'est pas connue. Un effort particulier a été fait pour que 848le compilateur affiche l'erreur dans le fichier source \mfront{} et 849non dans le fichier généré. Le code étant transformé par \mfront{}, le 850numéro de colonne n'est pas fiable. 851 852\subsubsection{Essai de compression hydrostatique} 853 854\paragraph{Description} Une des spécificités de la loi est d'être 855compressible. Pour vérifier cet aspect de la loi, nous considérons un 856essai de compression hydrostatique à pression imposée constante de 857\(70\,MPa\) durant \(1\) heure. La porosité initiale est de 858\(5\,\%\). La température est de \(293,15\,K\). 859 860\begin{figure}[htbp] 861 \centering 862 \lstinputlisting{@top_srcdir@/docs/tutorial/mtest/castemellipticcreep.mtest} 863 \caption{Modélisation d'un essai de compression hydrostatique avec 864 \mtest{}.} 865 \label{fig:tutorial:uo2:castemellipticcreep} 866\end{figure} 867 868\paragraph{Mise en données} La 869figure~\ref{fig:tutorial:uo2:castemellipticcreep} présente la mise en 870donnée utilisée. La principale difficulté est ici la définition de 871l'état initial (en contraintes, déformations totales et variables 872internes). 873 874\paragraph{Solution approchée} 875Afin de vérifier la solution obtenue, il est intéressant de disposer 876d'une solution analytique. Il s'avère que même dans ce cas très 877simple, les calculs deviennent rapidement rédhibitoires. Nous 878considérons donc une solution approchée considérant le coefficient 879\(A\paren{f}\) constant. La contrainte équivalente \(s\) est alors 880également constante et la variable \(p_{v}\) est linéaire en temps~: 881\[ 882p_{v}\paren{t}=\left.\deriv{\phi}{s}\right|_{t=0}\,\left.\deriv{s}{p}\right|_{t=0}\,\,t 883\] 884 885\begin{figure} 886 \centering 887 \includegraphics[width=12cm]{@top_srcdir@/docs/tutorial/mtest/castemellipticcreep-porosity.png} 888 \caption{Évolution de la porosité au cours de l'essai de compression hydrostatique.} 889 \label{fig:tutorial:uo2:mtest1:porosity} 890\end{figure} 891 892\paragraph{Comparaison} La 893figure~\ref{fig:tutorial:uo2:mtest1:porosity} compare la solution 894\mtest{} à la solution approchée. Les deux solutions coïncident bien 895au début de l'essai. Les deux solutions s'écartent en fin d'essai car 896l'approximation faite dans la solution approchée cesse d'être 897satisfaisante. 898 899\begin{figure} 900 \centering 901 \includegraphics[width=12cm]{@top_srcdir@/docs/tutorial/mtest/castemellipticcreep-convergence.png} 902 \caption{Convergence vers l'équilibre mécanique.} 903 \label{fig:tutorial:uo2:mtest1:convergence} 904\end{figure} 905 906\paragraph{Convergence} Il est intéressant d'étudier la convergence de 907l'algorithme d'équilibre de \mtest{}. Celui-ci est représenté en 908figure~\ref{fig:tutorial:uo2:mtest1:convergence} pour le dernier pas 909de temps. Il faut environ \(50\) itérations pour que la norme du 910résidu soit inférieure à la valeur du critère utilisé par \mtest{}. La 911vitesse de convergence est visiblement linéaire. Cela est dû au fait 912que \mtest{} utilise un opérateur tangent élastique (c'est le choix 913par défaut pour les lois compilées avec l'interface {\tt castem}). 914 915\begin{figure} 916 \centering 917 \includegraphics[width=12cm]{@top_srcdir@/docs/tutorial/mtest/castemellipticcreep-convergence-2.png} 918 \caption{Convergence vers l'équilibre mécanique avec l'algorithme 919 d'accélération de \castem{}.} 920 \label{fig:tutorial:uo2:mtest1:convergence2} 921\end{figure} 922 923\paragraph{Accélération de convergence} Il est possible d'activer dans 924\mtest{} l'algorithme d'accélération utilisée par les procédures de 925résolution de \castem{}~\cite{verpeaux_algorithmes_2014}. Il s'agit 926d'une généralisation de la méthode des sécantes, la plus simple à 927mettre en \oe{}uvre et la moins coûteuse (à notre connaissance). Pour 928cela, il suffit de rajouter la ligne suivante dans le fichier de 929donnée précédent~: 930\begin{center} 931 \lstinputlisting[firstline=10,lastline=10,numbers=none]{@top_srcdir@/docs/tutorial/mtest/castemellipticcreep-2.mtest} 932\end{center} 933 934La figure~\ref{fig:tutorial:uo2:mtest1:convergence2} compare les 935vitesses de convergence avec et sans l'utilisation de l'accélération 936de convergence. L'effet est net~: avec l'accélération de convergence, 937il ne faut que \(11\) itérations pour atteindre la précision requise 938(en fait, dans cet exemple, le résidu est indiqué nul à la précision 939machine absolue près (\(\approx{} 2,3\,10^{-308}\)). 940 941\paragraph{Le mode {\tt debug}} Il est 942également intéressant de regarder comment converge localement 943l'intégration de la loi. Pour cela, \mfront{} dispose d'une option 944{\tt -{}-{}debug}. Recompilons la loi~: 945\begin{flushleft} 946 {\tt 947 \$ \small mfront -{}-obuild -{}-debug -{}-interface=castem EllipticCreep.mfront 948 } 949\end{flushleft} 950 951Lors de l'exécution du cas test, la loi affiche certains détails sur 952sa convergence locale~: 953\begin{flushleft} 954 {\tt 955 EllipticCreep::integrate() : beginning of resolution\\ 956 EllipticCreep::integrate() : iteration 1 : 2.14848e-05\\ 957 EllipticCreep::integrate() : iteration 2 : 1.48821e-10\\ 958 EllipticCreep::integrate() : convergence after 2 iterations 959 } 960\end{flushleft} 961 962Cet exemple montre que l'algorithme de \nom{Newton} local converge en 963\(2\) itérations (la précision par défaut est de \(10^{-8}\)). 964 965\paragraph{Choix des paramètres numériques} Nous allons maintenant 966nous intéresser à deux paramètres numériques~: 967\begin{itemize} 968\item la valeur du critère de convergence \(\varepsilon\)~; 969\item la valeur de la perturbation \(\varepsilon_{p}\) pour le calcul 970 de la jacobienne numérique. 971\end{itemize} 972Pour comprendre ce que représente ce dernier paramètre, il faut 973rappeler comment la valeur de la jacobienne est approchée~: 974\[ 975J\paren{i,j}=\deriv{f_{i}}{\Delta\,x_{j}}\approx\Frac{f_{i}\paren{x_{j}+\varepsilon_{p}}-f_{i}\paren{x_{j}-\varepsilon_{p}}}{2\,\varepsilon_{p}} 976\] 977 978Par défaut, la valeur du critère de convergence \(\varepsilon\) est de 979\(10^{-8}\). Il est intéressant de remarquer qu'elle est plus lâche 980que la précision utilisée pour l'équilibre global. 981 982Dans la suite, nous verrons qu'il est nécessaire de baisser ce seuil 983pour obtenir une matrice tangente cohérente de bonne qualité, 984typiquement vers \(10^{-11}\). 985 986La valeur par défaut du second seuil, \(\varepsilon_{p}\), est dix 987fois inférieure à celle du critère de convergence. Dans notre exemple, 988nous avons donc une valeur de \(10^{-9}\), ce qui est une valeur 989acceptable dans la plupart des cas. Précisons que cette valeur est 990fixée à la compilation et n'est pas affectée par une modification de 991la valeur de \(\varepsilon\) à l'exécution. 992 993Ces deux valeurs peuvent être modifiées par les paramètres {\tt 994 epsilon} et {\tt numerical\_\-jaco\-bian\_\-epsilon}. Il est possible 995de modifier ces valeurs depuis \mtest{} à l'aide du mot clé {\tt 996 @Para\-meter}~: 997\begin{center} 998 \lstinputlisting[firstline=12,lastline=13,numbers=none]{@top_srcdir@/docs/tutorial/mtest/castemellipticcreep-2.mtest} 999\end{center} 1000 1001Avec un seuil à \(10^{-13}\) et une perturbation de \(10^{-9}\), nous 1002obtenons la sortie suivante~: 1003\begin{flushleft} 1004 {\tt 1005 EllipticCreep::integrate() : beginning of resolution\\ 1006 EllipticCreep::integrate() : iteration 1 : 2.14849e-05\\ 1007 EllipticCreep::integrate() : iteration 2 : 1.48417e-10\\ 1008 EllipticCreep::integrate() : iteration 3 : 9.80423e-16\\ 1009 EllipticCreep::integrate() : convergence after 3 iterations 1010 } 1011\end{flushleft} 1012 1013Avec un seuil à \(10^{-13}\) et une perturbation de \(10^{-18}\), nous 1014pouvons voir deux effets. En premier lieu, la qualité de la jacobienne 1015numérique est moindre, ce qui augmente le nombre d'itération dans 1016l'algorithme local. 1017\begin{flushleft} 1018 {\tt 1019 EllipticCreep::integrate() : beginning of resolution\\ 1020 EllipticCreep::integrate() : iteration 1 : 2.14849e-05\\ 1021 EllipticCreep::integrate() : iteration 2 : 4.17309e-08\\ 1022 EllipticCreep::integrate() : iteration 3 : 8.06865e-10\\ 1023 EllipticCreep::integrate() : iteration 4 : 6.27906e-11\\ 1024 EllipticCreep::integrate() : iteration 5 : 1.30492e-12\\ 1025 EllipticCreep::integrate() : iteration 6 : 1.78654e-14\\ 1026 EllipticCreep::integrate() : convergence after 6 iterations 1027 } 1028\end{flushleft} 1029Dans ce cas, on note que la valeur finale du résidu est moins bonne 1030que dans le cas précédent. Cela a un effet très fort sur la 1031convergence globale, qui nécessite dans ce cas \(53\) itérations (avec 1032une précision de \(10^{-12}\) le calcul ne converge même pas). Ceci 1033est dû au critère très serré utilisé par \mtest{} sur la valeur de 1034contraintes. Qualitativement, nous pouvons dire que la réponse de la 1035loi de comportement est moins \og~stable~\fg{} que précédemment, 1036malgré une valeur plus basse du critère de convergence. 1037 1038La diminution de la qualité de l'approximation de la dérivée avec la 1039valeur de la perturbation est un phénomène numérique 1040connu~\cite{fortin_analyse_2001}. 1041 1042Nous avons présenté ici un cas favorable où il faut donner une valeur 1043très basse à \(\varepsilon_{p}\) pour voir un effet notable. Pour 1044d'autres lois, des effets de ce type peuvent se faire sentir dès 1045\(10^{-12}\). Nous recommandons donc de toujours explicitement 1046préciser une valeur de l'ordre de \(10^{-9}\) à l'aide du mot clé {\tt 1047 @Pertubation\-Value\-For\-Numerical\-Jacobian\-Computation}~: 1048\begin{center} 1049 \lstinputlisting[firstline=10,lastline=11,firstnumber=10]{@top_srcdir@/docs/tutorial/mfront/EllipticCreep-2.mfront} 1050\end{center} 1051 1052% \subsubsection{Essai de traction uniaxiale} 1053% 1054% \subsection{Les hypothèses de contraintes planes} 1055% 1056% \subsubsection{L'hypothèse de contraintes planes} 1057% 1058% \subsubsection{L'hypothèse d'axisymétrie et de contraintes planes} 1059% 1060% \subsection{Optimisations numériques} 1061% 1062% \subsubsection{Matrice tangente cohérente} 1063% 1064% \subsubsection{Jacobienne analytique} 1065% 1066% \subsection{Essai de structure avec \castem{}} 1067 1068%%% Local Variables: 1069%%% mode: latex 1070%%% TeX-master: "tutoriel" 1071%%% End: 1072