xref: /dports/science/tfel/tfel-3.4.0/docs/mfront/behaviours/behaviours.tex.in

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% File      : behaviours.tex
% Author    : th202608@pleiades068.intra.cea.fr
% Date      : 15 oct. 2012
% Directory : /home/th202608/codes/tfel/tests/Broyden/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% \documentclass[rectoverso,pleiades,pstricks,leqno,anti]{note_technique_2010}
\documentclass[rectoverso,pleiades,pstricks,leqno,anti,projet]{note_technique_2010}

% \usepackage{draftcopy}
% \draftcopySetGrey{0.8}
% \draftcopyName{Version provisoire}{80}

\usepackage[dvips]{graphicx}
\usepackage[dvips,breaklinks]{hyperref}

\usepackage{mathematiques}
\usepackage{mecanique}
\usepackage{couleurs}
\usepackage{presentation}

\usepackage{pst-plot}
\usepackage{array}
\usepackage{subfigure}
\usepackage{relsize}
\usepackage{multind}

% one column index
\makeatletter
\def\printindex#1#2{\section*{#2}
\addcontentsline{toc}{section}{#2}
\@input{#1.ind}}
\makeatother

\usepackage[frenchb]{babel}

\newcommand{\pleiades}{\texttt{pleiades}}
\newcommand{\TFEL}{\texttt{tfel}}
\newcommand{\mfront}{\texttt{mfront}}
\newcommand{\mtest}{\texttt{mtest}}
\newcommand{\alcyone}{\texttt{alcyone}}
\newcommand{\germinal}{\texttt{germinal}}
\newcommand{\licos}{\texttt{licos}}
\newcommand{\celaeno}{\texttt{celaeno}}
\newcommand{\mfm}{\texttt{MFront\-Materials}}
\newcommand{\cyrano}{\texttt{cyrano}}
\newcommand{\galileo}{\texttt{galileo}}
\newcommand{\castem}{\texttt{Cast3M}}
\newcommand{\gibiane}{\texttt{gibiane}}
\newcommand{\tmfft}{\texttt{TMFFT}}
\newcommand{\aster}{\href{http://www.code-aster.org/}{\texttt{Aster}}}
\newcommand{\pycastem}{\texttt{pyCast3M}}
\newcommand{\umat}{\texttt{umat}}
\newcommand{\sirius}{\texttt{sirius}}
\newcommand{\fortran}{\texttt{fortran}}

\newcommand{\cmake}{\href{http://www.cmake.org/}{\texttt{cmake}}}
\newcommand{\autotools}{\href{http://fr.wikipedia.org/wiki/Autotools}{\texttt{autotools}}}
\newcommand{\python}{\href{http://python.org}{\texttt{python}}}
\newcommand{\gnuplot}{\href{http://www.gnuplot.info}{\texttt{gnuplot}}}
\newcommand{\latex}{\href{http://www.latex-project.org}{\LaTeX2e{}}}
\newcommand{\make}{\href{http://www.gnu.org/software/make/}{\texttt{make}}}
\newcommand{\doxygen}{\href{http://www.stack.nl/~dimitri/doxygen/}{\texttt{doxygen}}}
\newcommand{\valgrind}{\href{http://www.valgrind.org/}{\texttt{valgrind}}}

\newcommand{\mkey}[1]{\index{bmkeys}{#1@\symbol{64}#1}{\texttt{@#1}}}
\newcommand{\mkeyb}[2]{\index{bmkeys}{#1@\symbol{64}#1}{\texttt{@#2}}}

\newcommand{\header}[1]{\index{bheaders}{#1}{\texttt{#1}}}
\newcommand{\headerb}[2]{\index{bheaders}{#1}{\texttt{#2}}}

\newcommand{\tfel}[1]{\index{btfel}{#1}{\texttt{#1}}}
\newcommand{\tfelb}[2]{\index{btfel}{#1}{\texttt{#2}}}

%\newcommand{\env}[1]{\index{benv}{#1}{\texttt{#1}}}
%\newcommand{\envb}[2]{\index{benv}{#1}{\texttt{#2}}}

\newcommand{\moption}[1]{\texttt{-{}-#1}}

\newcommand{\bigO}[1]{\ensuremath{\mathop{}\mathopen{}O\mathopen{}\left(#1\right)}}

%c from texinfo.tex
\def\ifmonospace{\ifdim\fontdimen3\font=0pt }

%c C plus plus
\def\cpp{%
\ifmonospace%
    C++%
\else%
    C\kern-.1667em\raise.30ex\hbox{\smaller{++}}%
\fi%
\spacefactor1000 }

\newcommand{\varcpp}[1]{\texttt{#1}}

\newcommand{\sigmaH}{\ensuremath{\sigma_{H}}}

\newcommand{\nbzrc}{$NbZrC$}
\newcommand{\upuc}{$\paren{U,Pu}C$}
\newcommand{\sic}{$SiC$}

\newcommand{\cea}{CEA}
\newcommand{\windows}{\href{http://www.microsoft.com/france/windows/default.mspx}{\texttt{Windows}}}
\newcommand{\unix}{\href{http://www.kernel.org/}{\texttt{unix}}}
\newcommand{\msys}{\href{http://www.mingw.org/wiki/MSYS}{\texttt{msys}}}
\newcommand{\cygwin}{\href{http://www.cygwin.com/}{\texttt{cygwin}}}
\newcommand{\linux}{\href{http://www.kernel.org/}{\texttt{linux}}}
\newcommand{\debian}{\href{http://www.debian.org/}{\texttt{Debian}}}
\newcommand{\ubuntu}{\href{http://www.ubuntu.com}{\texttt{Ubuntu}}}
\newcommand{\redhat}{\href{http://www.redhat.com}{\texttt{Red Hat}}}
\newcommand{\mandriva}{\href{http://www.mandriva.com}{\texttt{Mandriva}}}
\newcommand{\excel}{\href{http://www.microsoft.com/france/office/2007/programs/excel/overview.mspx}{\texttt{Microsoft Office Excel}}}

\newcommand{\debutpas}[1]{\ensuremath{\left.#1\right|_{t}}}
\newcommand{\milieupas}[1]{\ensuremath{\left.#1\right|_{t+\theta\, \Delta\, t}}}
\newcommand{\finpas}[1]{\ensuremath{\left.#1\right|_{t+\Delta\, t}}}
\newcommand{\demipas}[1]{\ensuremath{\left.#1\right|_{t+\frac{\Delta\, t}{2}}}}

\newcommand{\uod}{\ensuremath{UO_{2}}}

\newcommand{\code}[1]{
  \psframebox[linecolor=ceaorange,shadow=true,blur=true]{
    \begin{minipage}[htbp]{1.0\linewidth}
      \ttfamily\small#1
    \end{minipage}
  }
}

\newcommand{\bash}[1]{
  \begin{center}
    \begin{minipage}{0.8\linewidth}
      \footnotesize{}
      \texttt{\$#1}
    \end{minipage}
  \end{center}
}

\input{LSC}

\auteurs{T.~Helfer, É. Castelier, V.~Blanc, J.~Julien}
\redacteur{T.~Helfer}
\verificateur{J.-M.~Proix}
\approbateur{R.~Masson}
\emetteur{E.~Touron}

\titre{Le générateur de code \mfront{}~: écriture de lois de
  comportement mécanique}

\date{Décembre 2013}
\numero{13-020}
\indice{0}
\dateversion{12/2013}
\numeroaffaire{A-SICOM-A1-01}
\domaine{DEN/DISN/SIMU}
\accords{tripartite}
\clients{AREVA - EDF}
\programmerecherche{SICOM}
\classification{DO}
\motsclefs{
  \mfront{} - \pleiades{}
}

% \codebarre{images/code_barre}
\diffusionexterne{
{EDF/R\&D}              & O. Marchand     & 1 & Diffusion par\\
{EDF/R\&D}              & P. Vasseur      & 1 & courriel     \\
{EDF/R\&D/AMA}          & É. Lorentz      & 1 & \\
                        & C. Durand       & 1 & \\
{EDF/R\&D/AMA/T64}      & T. de Soza      & 1 & \\
                        & J. Delmas       & 1 & \\
                        & J.M. Proix      & 1 & \\
                        & F. Hammon       & 1 & \\
                        & N. Sellenet     & 1 & \\
{EDF/R\&D/MMC}          & P. Ollar        & 1 & \\
{EDF/R\&D/MMC/MAESTRO}  & N. Rupin        & 1 & \\
                        & F. Latourte     & 1 & \\
{EDF/R\&D/MMC/CPM}      & N. Prompt       & 1 & \\
                        & N. Barnel       & 1 & \\
{EDF/R\&D/MMC/CPM}      & G. Thouvenin    & 1 & \\
                        & F. Douchin      & 1 & \\
                        & R. Largenton    & 1 & \\
                        & C. Petry        & 1 & \\
EDF/SEPTEN              & N. Waeckel      & 1 & \\
                        & C. Chauliac     & 1 & \\
                        & H. Billat       & 1 & \\
                        & C. Bernaudat    & 1 & \\
AREVA NP/LA DEFENSE     & L. Catalani     & 1 & \\
                        & L. Brunel       & 1 & \\
AREVA NP/LYON           & P. Melin        & 1 & \\
                        & V. Bessiron     & 1 & \\
                        & C. Garnier      & 1 & \\
                        & V. Garat        & 1 & \\
                        & F. Arnoux       & 1 &
}
\CoupeListeDiffusion{}
\diffusioninterne{
  DEN/DISN/SIMU       & J.P. Deffain       & 1 & Diffusion par\\
                      & D. Caruge          & 1 & courriel     \\
  DEN/DM2S/SEMT       & X. Averty          & 1 & \\
  DEN/DM2S/SEMT/LM2S  & J.L. Fayard        & 1 & \\
                      & P. Verpeaux        & 1 & \\
                      & A. Millard         & 1 & \\
                      & S. Pascal          & 1 & \\
                      & O. Fandeur         & 1 & \\
  DEN/DMN             & P. Yvon            & 1 & \\
                      & J.L. Seran         & 1 & \\
                      & F. Dalle           & 1 & \\
  DEN/DMN/SRMA        & P. Chapelot        & 1 & \\
                      & S. Carassou        & 1 & \\
                      & B. Marini          & 1 & \\
  DEN/DMN/SRMA/LC2M   & L. Nicolas         & 1 & \\
                      & J. Garnier         & 1 & \\
                      & S. Vincent         & 1 & \\
                      & L. Vincent         & 1 & \\
                      & L. Gelebart        & 1 & \\
                      & M. Sauzay          & 1 & \\
                      & L. Dupuy           & 1 & \\
                      & P. Forget          & 1 & \\
                      & A. Hellouin de Menibus  & 1 & \\
                      & M. Le Saux         & 1 & \\
                      & C. Robertson       & 1 & \\
  DEN/DMN/SRMA/LA2M   & J.-L. Bechade      & 1 & \\
  DEN/DMN/SRMA/LTMEX  & L. Chaffron        & 1 & \\
                      & D. Sornin          & 1 & \\
  DEN/DMN/SEMI        & C. Poussard        & 1 & \\
                      & B. Tanguy          & 1 & \\
  DEN/DMN/SEMI/LCMI   & V. Vandenberghe    & 1 & \\
                      & A. Courcelle       & 1 & \\
                      & F. Hure            & 1 & \\
                      & D. Leboulch        & 1 & \\
                      & Q. Auzoux          & 1 & \\
                      & Y. Robert          & 1 & \\
  DEN/DER/SESI/LE2S   & P. Lamagnère       & 1 & \\
                      & D. Gentet          & 1 & \\
                      & Y. Lejeail         & 1 & \\
                      &                    &  & \\
  DEN/DEC             & P. Brossard        &  & Document disponible\\
  DEN/DEC/PROJETS     & P. Obry            &  & sur intradec\\
  DEN/DEC/SESC        & E. Touron          &  & \\
                      & M. Casella         &  & \\
                      & M. Agard           &  & \\
  DEN/DEC/SESC/LIPA   & C. Nonon-Solaro    &  & \\
                      & C. Bassi           &  & \\
                      & O. Bremond         &  & \\
  DEN/DEC/SESC/LLCC   & V. Basini          &  & \\
                      & J.-M. Escleine     &  & \\
  DEN/DEC/SESC/LC2I   & D. Plancq          &  & \\
                      & S. Bejaoui         &  & \\
                      & V. Blanc           &  & \\
                      & T. Beck            &  & \\
                      & F. Biscarrat       &  & \\
                      & D. Lorenzo         &  & \\
                      & I. Guénot-Delahaie &  & \\
                      & P. Masoni          &  & \\
                      & B. Valentin        &  & \\
                      & M. Zabiego         &  & \\
  DEN/DEC/SESC/LSC    & R. Masson          &  & \\
                      & B. Michel          &  & \\
                      & M. Pelletier       &  & \\
                      & M. Lainet          &  & \\
                      & V. Bouineau        &  & \\
                      & V. Marelle         &  & \\
                      & T. Helfer          &  & \\

}

% \signatures{-0.}{-39.2}{0.12}{images/signatures.eps}

\stylebib{@abs_top_srcdir@/docs/tex/texmf/bibtex/fr-insa}
\fichierbib{@abs_top_srcdir@/docs/tex/texmf/bibtex/bibliographie}

\resumecea{
  L'une des missions de la plate-forme \pleiades{} est de capitaliser
  les connaissances matériau utilisées pour la simulation des éléments
  combustibles et absorbants.

  Cette note s'intéresse aux lois de comportement mécanique et plus
  particulièrement à leur écriture à l'aide du générateur de code
  \mfront{}. Elle complète la présentation générale de \mfront{} qui
  traite des propriétés matériau et des modèles
  physico-chimiques~\cite{helfer_generateur_2013-1}. Dans le cadre de
  la plate-forme \pleiades{}, ces lois de comportement peuvent être
  stockées dans la base de données \sirius{}.

  La variété des phénomènes traités fait que \mfront{} propose
  différentes façons d'écrire des lois de comportement mécaniques.

  Certaines sont spécifiques à des lois de comportement qui sont
  d'usage courant et pour lesquelles des algorithmes d'intégration
  performants existent (plasticité et viscoplasticité isotropes
  notamment).

  Pour les lois les plus complexes, \mfront{} permet d'utiliser
  des méthodes de \nom{Runge-Kutta} ou des méthodes implicites.

  % Les différents tests menés avec les codes aux éléments finis
  % \castem{} et \aster{} montrent que les lois de comportement générées
  % sont numériquement très
  % efficaces~\cite{michel_etude_2009,helfer_bilan_2011,proix_integration_2013,olagnon_analysis_2013}.
  % Il s'agit d'un outil mature, qui se veut de qualité industrielle et
  % qui intéresse fortement différents départements du
  % \cea{}~\cite{helfer_presentation_2011,darrigo_notice_2012} ainsi que ses
  % partenaires, Areva~\cite{olagnon_analysis_2013} et
  % ÉDF~\cite{proix_integration_2013}. Cet intérêt confirme la crédibilité et
  % la rigueur de la gestion des connaissances matériau au sein de la
  % plate-forme. Le partage de cet outil ouvre des perspectives variées
  % et enrichissantes de co-développement et de mutualisation.

}

%\makeindex{env}
\makeindex{btfel}
\makeindex{bheaders}
\makeindex{bmkeys}

\begin{document}

\clearpage
\newpage
\section{Introduction}

Ce document décrit comment écrire des lois de comportement mécanique
avec le générateur de code \mfront{}. \mfront{} vise à garantir une
gestion pérenne, robuste, efficace et évolutive des connaissances
matériau dans la plate-forme
\pleiades{}~\cite{michel_etude_2009,helfer_bilan_2011}. Il permet
également à des utilisateurs non développeurs d'écrire leurs propres
connaissances matériau~\cite{helfer_ajout_2010}. Cette note intègre
des éléments d'une note antérieure décrivant un algorithme numérique
particulier~\cite{blanc_integration_2011}.

Elle complète la présentation générale de
\mfront{}~\cite{helfer_generateur_2013-1}. Deux autres notes
traitent\footnote{Citons également qu'un guide de référence de la
  librairie \TFEL{} est également en cours de rédaction. Cette note
  étant assez ambitieuse, sa rédaction s'étalera dans le temps. Elle
  est cependant disponible en version projet dans les sources de la
  librairie. Elle contient en particulier le guide d'installation de
  la librairie.}~:
\begin{itemize}
\item de l'interface {\tt umat} utilisée par le codes aux éléments
  finis \castem{}~\cite{helfer_interface_2013}~;
\item de l'interface utilisée par le codes aux éléments finis
  \aster{}~\cite{helfer_interface_2013-1}~;
\end{itemize}
Ces notes sont intégrées à la gestion de configuration de \TFEL{} et
évoluent continûment avec les développements de \TFEL{}\footnote{ La
  version la plus à jour est obtenue dans le répertoire de compilation
  des sources par la commande~:
\begin{center}
  {\tt make doc-pdf}
\end{center}}. Le présent document a été généré à partir de la
révision @TFEL_SVN_REVISION@.

La lecture de ce document suppose que le lecteur est déjà familier de
\mfront{}. A minima, la lecture de la présentation générale de
\mfront{} semble nécessaire~\cite{helfer_generateur_2013-1}.

\paragraph{De multiples phénomènes} Les matériaux solides réagissent
aux sollicitations mécaniques par différents phénomènes~: élasticité,
viscoplasticité, plasticité, endommagement. Nous renvoyons aux
ouvrages classiques pour la description physique de ces
phénomènes~\cite{francois_comportement_1995,chaboche_mecanique_2009,besson_mecanique_2001}.

\paragraph{Rôle de la loi de comportement} Nous décrivons en
annexe~\ref{sec:mfront:mechanical_equilibrium} un algorithme simplifié
de recherche de l'équilibre mécanique statique non-linéaire qui
précise la place de la loi de comportement mécanique. Le lecteur
intéressé pourra se reporter aux documentations des codes éléments
finis pour une description plus
précise~\cite{pascal_notice_2005,abbas_algorithme_2013,besson_mecanique_2001}.

En résumé, connaissant l'état mécanique du matériau à un instant
\(t\), les lois de comportement doivent, en réponse à un incrément de
sollicitation mécanique représentée par un incrément de déformations
totales \(\Delta\,\tepsilonto\), fournir, en chaque point
d'intégration~:
\begin{itemize}
\item l'évolution microstructurelle du matériau, décrit par un
  ensemble de variables internes \(y_{i}\) sur le pas de temps
  \(\Delta\, t\)~;
\item la valeur de la contrainte en fin de pas de temps (ou de manière
  équivalente l'incrément de contraintes \(\Delta\,\tsigma\)).
\end{itemize}

Si l'algorithme de recherche de l'équilibre mécanique le nécessite, la
loi de comportement peut également fournir~:
\begin{itemize}
\item la matrice tangente cohérente
  \(\deriv{\Delta\tsigma}{\Delta\tepsilonto}\) ou une estimation de
  celle-ci.
\item une matrice de prédiction en début de pas de temps permettant de
  fournir une première estimation de la solution en fin de pas de
  temps, avant de débuter la résolution.
\end{itemize}

Dans certains cas, les lois de comportement peuvent également donner
des indications sur la qualité de la discrétisation temporelle du
problème~\cite{proix_etude_2011}.

\subsection{Analyseurs disponibles}

L'écriture de lois de comportement nécessite de prendre à la fois en
compte les mécanismes physiques décrits (certains algorithmes sont
plus adaptés à certains phénomènes) et des besoins des codes éléments
finis (matrice tangente cohérente). La variété des situations
rencontrées en pratique expliquent la variété des analyseurs proposés
par \mfront{}, qui sont actuellement au nombre de \(8\)\footnote{La
  notion d'analyseur est décrite plus en détails dans la première
  partie de la documentation de
  \mfront~\cite{helfer_generateur_2013-1}.}. La description de ces
analyseurs est l'objet de ce document.

Plusieurs analyseurs dédiés aux lois de comportement sont
actuellement disponibles~:
\begin{itemize}
\item l'analyseur \texttt{DefaultParser} qui permet de traiter
  tous types de lois de comportements~;
\item l'analyseur \texttt{IsotropicMisesCreep} qui gère
  exclusivement les lois de comportement mécanique viscoplastique
  incompressible sans écrouissage des matériaux isotropes~;
\item l'analyseur \texttt{IsotropicStrainHardeningMisesCreep} qui gère
  exclusivement les lois de comportement mécanique viscoplastique
  incompressible avec écrouissage des matériaux isotropes~;
\item l'analyseur \texttt{IsotropicPlasticMisesFlow} qui gère
  exclusivement les lois de comportement mécanique plastique
  incompressible des matériaux isotropes~;
\item l'analyseur \texttt{MultipleIsotropicMisesFlows} qui gère
  une combinaison arbitraire d'écoulements des trois types précédents.
  Les différents écoulements sont supposés non couplés~;
\item les analyseurs \texttt{Implicit} et \texttt{ImplicitII} qui
  simplifient la résolution d'une loi de comportement mécanique
  quelconque à l'aide d'une intégration implicite.
\end{itemize}

Ces différents analyseurs sont décrits dans les sections suivantes.

\subsection{Plan de la note}

La section~\ref{mfront:behaviours:gen} décrit des généralités sur les
lois de comportements mécaniques et introduit certaines notions
nécessaires à la suite.

Les analyseurs \texttt{Isotropic\-Mises\-Creep},
\texttt{IsotropicStrainHardening\-Mises\-Creep},
\texttt{Isotropic\-Plastic\-Mises\-Flow} et
\texttt{Multiple\-Isotropic\-Mises\-Flows} sont très proches et sont
décrits en section~\ref{sec:mfront:isotropic:analyser}.

Les analyseurs \texttt{DefaultParser}, \texttt{RungeKutta} et
\texttt{Implicit} sont dits {\em génériques} car ils permettent de
traiter n'importe quelle loi de comportement. Ils sont respectivement
décrits dans les sections~\ref{sec:defaultparser}, \ref{sec:RK}
et~\ref{sec:Implicite}.

\clearpage
\newpage
\section{Généralités}
\label{mfront:behaviours:gen}

Nous abordons dans cette section des points qui sont utiles pour la
lecture de la suite. Nous commençons par préciser les définitions
utilisées dans la suite. Nous donnons ensuite quelques conseils sur le
choix de l'analyseur à utiliser. Nous traitons quelques points qui
sont indépendants de l'analyseur utilisé. Enfin, nous précisons le
rôle des interfaces aux lois de comportement.

\subsection{Quelques définitions}

Les lois de comportement peuvent être complexes. Pour les décrire, il
nous faut introduire quelques définitions.

\paragraph{Propriétés matériau} Afin de pouvoir adapter des lois de
comportements à différents matériaux, celles-ci peuvent utiliser des
propriétés matériau, qui sont définies, dans \mfront{}, comme des
fonctions des valeurs actuelles des variables d'état du
matériau~\cite{helfer_generateur_2013-1}.

Les lois de comportement mécaniques peuvent~:
\begin{itemize}
\item demander à ce qu'un certain nombre de propriétés matériau leur
  soient fournies par le code appelant~\footnote{Voir le mot clé
    \mkey{MaterialProperty}.}. Dans certains cas, le code \castem{}
  notamment~\cite{helfer_interface_2013}, le code appelant impose que la
  loi de comportement utilise des propriétés matériau prédéfinies.
\item utiliser des lois définies dans des fichiers \mfront{}
  dédiés~\footnote{Voir le mot clé \mkey{MaterialLaw}.}.
\end{itemize}

\paragraph{Les variables internes} Les variables internes décrivent
l'état mécanique local du matériau.

Pour \mfront{}, les variables internes peuvent être soit des {\em
  scalaires} soit des {\em tenseurs} d'ordre \(2\) symétriques.

Certains analyseurs, dédiés à des lois de comportement spécifiques, ne
permettent pas de déclarer de nouvelles variables internes.

\paragraph{Les variables internes auxiliaires} Les variables internes
auxiliaires désignent des variables internes qui ne sont pas nécessaires
pour l'intégration de la loi de comportement. Ces variables ont des
utilités diverses~:
\begin{itemize}
  \item elles peuvent désigner des variables qui peuvent être éliminées
  de l'intégration~;
  \item des variables uniquement destinées aux posttraitements.
\end{itemize}

Les variables internes auxiliaires sont mises à jour après
l'intégration des variables internes. Elles peuvent être soit des {\em
  scalaires} soit des {\em tenseurs}.

\paragraph{Les variables externes} Les variables externes
désignent des variables dont l'évolution est donnée par ailleurs et
connue sur le pas de temps. Ces variables peuvent ou être des variables
d'état du matériau ou des paramètres externes (flux de neutrons,
fluence, densité de fissions).

Parmi les variables externes, nous pouvons citer la température, la
déformation totale du matériau (qui représente la sollicitation locale
et dont la valeur est donnée par la résolution de l'équilibre global du
matériau). Ces deux variables sont traitées de manière particulière par
\mfront{} et déclarées automatiquement.

\paragraph{Les variables locales} Les variables locales permettent
généralement de calculer des variables avant de débuter l'intégration
afin d'éviter des calculs superflus.

Une utilisation typique de variable locale est de calculer avant
l'intégration des termes de type \nom{Arrhenius} (termes de la forme
\(\exp\paren{-\Frac{Q}{R\,T}}\)) afin de ne pas les réévaluer au cours
des itérations. Ces termes sont souvent évalués en milieu de pas de
temps, ce qui est cohérent avec une intégration par une
\(\theta\)-méthode, et une approximation généralement suffisante pour
les autres méthodes d'intégration.

Il n'y a pas de limite sur le type des variables locales.

La directive \mkey{InitLocalVariables}\footnote{La directive
  \mkey{InitLocalVars} est synonyme de la directive
  \mkey{InitLocalVariables}.} permet d'initialiser ces variables
locales.

\paragraph{Tableaux de variables internes, de propriétés matériau et
  de variables externes} Afin de pouvoir regrouper des équations dont
le {\em formalisme} était similaire, nous avons introduit la
possibilité de définir des tableaux de variables internes, de
propriétés matériaux et de variables externes. Il faut noter que la
taille de ces tableaux, c'est à dire le nombre de variables internes,
de propriétés matériau et/ou de variables externes est fixés \og~en
dur~\fg{} dans le fichier d'entrée.

Il est alors possible d'écrire des boucles sur les variables internes
constituant le tableau pour condenser l'écriture des lois.

La possibilité d'utiliser des tableaux de variables internes est
particulièrement utile en homogénéisation. Les lois de comportement
homogénéisées peuvent avoir un grand nombre de variables internes qui
partagent des lois d'évolutions similaires. Un exemple de cela est
donné par les lois issues de l'homogénéisation de poly-cristaux qui
peuvent conduire à plusieurs milliers de variables
internes~\cite{proix_comportements_2013-1}. Grâce aux tableaux de variables
internes, l'implantation de ce type de lois peut être très courte (une
centaine de lignes).

\subsection{Conseils sur le choix de l'analyseur à utiliser}

Pour pouvoir traiter tous les types de lois de comportement,
différents analyseurs sont disponibles. Les analyseurs propres à des
formes de lois de comportement sont dits spécifiques. Par opposition,
les autres sont dits génériques.

Nous pouvons donner quelques conseils généraux sur le choix de
l'analyseur à utiliser~:
\begin{itemize}
\item si un intégrateur spécifique existe, il vaut mieux l'utiliser~:
  il utilise un algorithme optimisé et robuste et le nombre
  d'informations à fournir est réduit.
\item si l'on doit recourir à un intégrateur générique, il vaut mieux
  préférer l'intégration implicite, surtout s'il s'agit de lois
  indépendantes du temps (plasticité, endommagement)~:
  \begin{itemize}
    \item l'équation différentielle pour la plasticité ou
    l'endommagement doit être remplacée par la nullité du critère en fin
    de pas de temps~;
    \item les temps de calculs sont souvent {\em très} avantageux~;
  \end{itemize}
  \item il ne faut utiliser l'analyseur \nom{Runge-Kutta} que~:
  \begin{itemize}
    \item si {\em vraiment rien d'autre} n'est possible (impossibilité
    de calculer la jacobienne)
    \item si le temps de développement est limité~;
    \item si le temps d'intégration de la loi de comportement ne pose
    pas de problème de performance~;
  \item si le nombre de variables internes est très grand (loi de
    comportement issus de l'homogénéisation de
    poly-cristaux~\cite{proix_comportements_2013-1} notamment).
  \end{itemize}
\end{itemize}

Insistons sur le fait que l'utilisation de l'analyseur
\nom{Runge-Kutta} est fortement déconseillée. De nombreux
développements ont été faits pour rendre l'utilisation des méthodes
implicites plus simples (calcul automatique de la jacobienne par
différentiation numérique, algorithmes de \nom{Broyden}) et si
l'expression de la loi reste encore un peu plus complexe, les gains en
performance valent largement l'effort supplémentaire.

L'analyseur \texttt{DefaultParser} n'est utile que dans ces cas très
particuliers, par exemple pour des lois dépendant explicitement des
déformations (loi de \nom{Mazars} par
exemple~\cite{hamon_modeendommagement_2013}).

\subsection{Points particuliers}

\subsubsection{Hypothèses de modélisation}

Les lois de comportement sont implantées par des classes {\tt
  template} paramétrée par l'hypothèse de modélisation. Ce choix
permet de~:
\begin{itemize}
\item proposer une implantation optimisée et fiable de la loi de
  comportement pour chacune des hypothèses de modélisation~;
\item traiter les cas particuliers, les contraintes planes notamment~;
\item de spécifier quelles sont les hypothèses valides. Par exemple,
  les lois issues de l'homogénéisation de poly-cristaux doivent
  nécessairement être intégrées en
  \(3D\)~\cite{proix_comportements_2013-1}
\end{itemize}

Le fichier d'entête
\headerb{TFEL/MaterialLaw/ModellingHypothesis.hxx}{TFEL/\-Material\-Law/\-Modelling\-Hypo\-thesis\-.hxx}
définit une structure
\tfelb{ModellingHypothesis}{Modelling\-Hypo\-thesis}. Dans cette
structure, un objet de type énumération nommé
\tfelb{ModellingHypothesis::Hypothesis}{Hypo\-thesis} définit les
différentes hypothèses de modélisation ajourd'hui supportées~:
\begin{itemize}
 \item \texttt{AXISYMETRICAL\-GENERALISED\-PLANE\-STRAIN} qui désigne
 une modélisation \(1D\) axisymétrique plan généralisée. Dans ce cas,
 le tenseur des contraintes a \(3\) composantes et est représenté par
 le vecteur suivant~:
\[
 \tsigma = \left(
   \begin{array}{c}
     \sigma_{rr} \\
     \sigma_{zz} \\
     \sigma_{\theta\theta} \\
   \end{array}
   \right)
 \]
 \item \texttt{AXI\-SYMETRICAL} qui désigne une modélisation \(2D\)
 axisymétrique. Dans ce cas, le tenseur des contraintes a \(4\)
 composantes et est représenté par le vecteur suivant~:
\[
 \tsigma = \left(
   \begin{array}{c}
     \sigma_{rr} \\
     \sigma_{zz} \\
     \sigma_{\theta\theta} \\
     \sqrt{2}\sigma_{rz} \\
   \end{array}
   \right)
 \]
 \item \texttt{PLANE\-STRESS}, \texttt{PLANE\-STRAIN} et
 \texttt{GENERALISED\-PLANE\-STRAIN} qui désigne différentes
 modélisations \(2D\) qui se distinguent par le traitement de la
 direction axiale. Dans ces cas, le tenseur des contraintes a \(4\)
 composantes et est représenté par le vecteur suivant~:
\[
 \tsigma = \left(
   \begin{array}{c}
     \sigma_{xx} \\
     \sigma_{yy} \\
     \sigma_{zz} \\
     \sqrt{2}\sigma_{xy} \\
   \end{array}
   \right)
 \]
\item \texttt{TRIDIMENSIONAL} qui désigne la modélisation la plus
  générale. Dans ce cas, le tenseur des contraintes a \(6\)
  composantes et est représenté par le vecteur suivant~:
\[
 \tsigma = \left(
   \begin{array}{c}
     \sigma_{xx} \\
     \sigma_{yy} \\
     \sigma_{zz} \\
     \sqrt{2}\sigma_{xy} \\
     \sqrt{2}\sigma_{xz} \\
     \sqrt{2}\sigma_{yz} \\
   \end{array}
   \right)
 \]
\end{itemize}

\paragraph{Les directives \mkey{ModellingHypothesis} et
  \mkey{ModellingHypotheses}} Les directives
\mkeyb{ModellingHypothesis}{Modelling\-Hypothesis} et
\mkeyb{ModellingHypotheses}{Modelling\-Hypotheses} permettent de
spécifier les hypothèses de modélisation valides. La première attend
le nom d'une hypothèse, la seconde un tableau contenant une ou
plusieurs hypothèses.

Les noms d'hypothèses valides sont donc les suivants~:
\begin{itemize}
\item AxisymmetricalGeneralisedPlaneStrain~;
\item Axisymmetrical~;
\item PlaneStress~;
\item PlaneStrain~;
\item GeneralisedPlaneStrain~;
\item Tridimensional.
\end{itemize}

\paragraph{Cas des contraintes planes}
Le cas des contraintes planes appelle certaines remarques. Le plus
souvent, les lois de comportements ne peuvent être utilisées telles
quelles en contraintes planes~: une implantation spécifique doit être
faite. Pour cette raison, les contraintes planes sont exclues des
hypothèses de modélisation supportées par défaut.

Plusieurs stratégies peuvent être mises en place pour éviter d'avoir à
écrire une implantation spécifique~:
\begin{itemize}
\item le code aux éléments finis {\tt Zebulon} introduit des éléments
  finis spécifiques, possédant un degré de liberté supplémentaire (la
  déformation axiale) dont la force associée (liée à la contrainte
  axiale) est nulle~\cite{besson_object-oriented_1998-1}~;
\item le code aux éléments finis \aster{} modifie l'algorithme de
  résolution global pour converger vers une solution en contraintes
  planes~\cite{proix_prise_2012}.
\end{itemize}

Pour le code aux éléments finis \castem{}, si une implantation en
contraintes planes n'est pas disponible, l'interface \umat{} fournie
par \mfront{} gère les contraintes planes en introduisant une variable
interne supplémentaire (la déformation axiale) qu'elle ajuste (par
appels successifs à la loi de comportement \(2D\)) pour trouver une
contrainte axiale nulle~\cite{helfer_interface_2013}.

\subsubsection{Traitement des matériaux orthotropes}

Ce paragraphe décrit des fonctionnalités qui ne sont pas accessibles aux
analyseurs spécifiques (dédiés à des lois isotropes).

Écrire les lois de comportement orthotrope de manière indépendante
de l'hypothèse de modélisation est rendue difficile par certaines
restrictions imposées par les codes aux éléments finis \castem{} et
\aster{}~:
\begin{itemize}
  \item il est nécessaire de gérer la rotation des contraintes
  dans le repère propre au matériau.
\item il n'est pas possible de choisir de manière cohérente une
  convention sur l'ordonnancement des axes valable quelle que soit la
  dimension. Le cas des tubes est traité en détail dans
  l'annexe~\ref{sec:annexe:orthotropie}.
\end{itemize}

\subsubsection{Traitement des dilatations libres}

\mfront{} ne traite en général pas des dilatations libres, et plus
particulièrement ne traite en général pas des dilatations
thermiques~:
\begin{enumerate}[-]
\item en petites transformations, Dans la plupart des codes aux
éléments finis~\cite{cea_site_2013,edf_site_2013}, les dilatations
libres sont généralement traitées en amont de l'intégration de la loi
de comportement. Cela se traduit par le fait que la déformation totale
passée à la loi de comportement est en fait la déformation mécanique
totale, c'est à dire la déformation totale à laquelle ont été
retranchées toutes les dilatations libres.
\item en grandes transformations, le traitement des dilatations libres
  dépend du formalisme choisi et en cadre général ne peut être fourni.
\end{enumerate}

Il peut cependant être intéressant de traiter certains gonflements dans
la loi de comportement. Donnons quatres exemples pratiques~:
\begin{enumerate}[-]
\item certaines lois de comportement viscoplastique de gaine de
  réacteurs refroidis au sodium relient l'intensité de l'écoulement
  viscoplastique à la vitesse de gonflement~;
\item certaines modélisations plastiques des aciers considèrent la
  température comme une variable interne~: l'évolution de la
  température, supposée adiabatique le temps de l'expérience, est
  supposée due à la dissipation mécanique locale. Il peut paraître
  intéressant de coupler finement mécanique et thermique~;
\item il peut être nécessaire de traiter des changements de phases au
  cours du calcul mécanique, en particulier la transition
  martensite-austénite de certains aciers qui dépend de l'état de
  contraintes~\cite{milliard_mechanical_2014}~;
\item certaine stratégies \og~grandes transformations~\fg{} permettent
  d'utiliser un formalisme de lois en petites déformations, voir
  d'utiliser les lois identifiées en petites déformations
  (voir~\cite{helfer_ecriture_2014}). Dans ce cas, la dissymétrie de
  traitement par les codes appelant (notée plus haut) entre les
  petites déformations (traitement des dilatations libres) et les
  grandes déformations (non traitement des dilatations libres)
  apparaît d'autant plus gênante que l'expression usuelle des
  dilatations libres (voir équation~\eqref{eq:behaviour:epsilonth}
  pour la dilatation thermique), convenablement post-traitée, peut
  être réutilisée.
\end{enumerate}

\paragraph{Cas des lois en petites transformations, cas général}

En petites transformations, une solution simple pour prendre en compte
une dilatation libre, est de calculer un incrément de dilatation et de
le soustraire à l'incrément de dilatation totale dans la partie
\mkey{InitLocalVariables}. Cette stratégie n'est valide que dans le
cas où la dilatation libre ne dépend pas du résultat de la mécanique

Il est possible de garder une trace de cette dilatation en définissant
une variable auxiliaire associée. Il est cependant nécessaire de prendre
garde à la manière dont sont mises à jour les variables auxiliaires par
l'analyseur utilisé.

\paragraph{Cas des lois en petites transformations, dilatation
  thermique} La dilatation thermique est un cas particulier par son
coté systématique~: il suffit de savoir calculer le coefficient de
dilatation moyen \(\alpha\paren{T}\) que l'on suppose ne dépendre que
de la température actuelle \(T\). La dilatation thermique linéique
s'écrit alors (voir~\cite{helfer_ecriture_2014})~:
\begin{equation}
  \label{eq:behaviour:epsilonth}
    \Frac{\Delta\,l}{l^{i}}\paren{T}=
    \Frac{1}{1+\alpha\paren{T^{i}}\paren{T^{i}-T^{\alpha}}}\,\left[\alpha\paren{T}\paren{T-T^{\alpha}}-\alpha\paren{T^{i}}\paren{T^{i}-T^{\alpha}}\right]
\end{equation}
où~:
\begin{minipage}[t]{0.8\linewidth}
  \begin{enumerate}[-]
  \item \(T^{\alpha}\) est la température de référence pour
    l'expérience de dilatométrie ayant servie à identifier le
    coefficient de dilatation thermique et qui correspond à une
    dilatation nulle~;
  \item \(T^{i}\) est la température à laquelle la géométrie du corps
    a été mesurée (température de début de calcul)~;
  \item \(\Delta\, l\) est la variation de longueur par rapport à la
    température \(T^{i}\).
  \end{enumerate}
\end{minipage}

La directive \mkey{ComputeThermalExpansion} permet le calcul de la
dilatation thermique par la formule~\eqref{eq:behaviour:epsilonth} en
amont de l'intégration de la loi de comportement.

Elle prend un nom de fichier \mfront{} argument dans le cas isotrope,
ou un tableau de trois noms de fichiers \mfront{} dans le cas
orthotrope (une dilatation pour chaque direction principale). Ces
fichiers doivent décrire le coefficient de dilatation thermique moyen
comme une propriété matériau.

La température \(T^{\alpha}\) est supposée donnée dans le fichier
définissant la dilatation par la définition d'une constante nommée
{\tt Reference\-Temperature} de la directive \mkey{Constant}. Si cette
température n'est pas fournie, \mfront{} prendra la température
ambiante (\(293,15\,K\)).

La température \(T^{i}\) est automatiquement déclarée comme un
paramètre nommé {\tt
  reference\-Temperature\-For\-Thermal\-Expansion}\footnote{La notion
  de paramètres est décrite dans la notice générale de
  \mfront{}~\cite{helfer_generateur_2013-1}.}. Par défaut, sa valeur
est de \(293,15\,K\) (température ambiante).

\paragraph{Cas particuliers} Si le formalisme \og~grandes
transformations~\fg{} choisi se base sur une écriture des lois écrites
dans le formalisme des petites transformations (utilisation d'une
hypothèse de \og~petites déformations, grandes rotations~\fg,
déformations logarithmiques), il est souvent possible de proposer un
traitement des dilatations libres systématique
(voir~\cite{helfer_ecriture_2014}).

\paragraph{Implantation} Pour permettre aux interfaces de proposer des
formalismes \og~grandes transformations\fg{} réutilisant le formalisme
des \og~petites transformation~\fg{}, le calcul de la dilatation doit
être fait en amont de la loi.

Il revient aux interfaces de modifier le chargement du point matériel
de manière adéquate. Par exemple, en petites transformations, la
déformation totale envoyée à la loi de comportement sera la
déformation totale fournie par le code appelant à laquelle on aura
retiré la dilatation.

%\subsubsection{Classes générées}
%
%\begin{figure}[htbp]
%  \centering
%  \includegraphics[height=15cm]{@top_srcdir@/docs/mfront/Images/NortonRK.eps}
%  \caption{Classes générées pour les lois de comportement (diagramme
%    généré par l'outil \doxygen{}).}
%  \label{fig:nortonrk:inheritance}
%\end{figure}
%
%Hormis ce qui est propre aux interfaces, un fichier \mfront{} génère des
%classes \cpp{}. Pour une loi de comportement nommmée \texttt{NortonRK},
%trois classes par loi de comportement sont générées~:
%\begin{itemize}
%  \item une classe \texttt{NortonRKBehaviourData} qui contient la valeur
%  des propriétés matériau, l'état des variables internes et externes en
%  début de pas de temps~;
%  \item une classe \texttt{NortonRKIntegrationData} qui contient la
%  valeur des incréments des variables externes~;
%  \item une classe \texttt{NortonRK} qui contient les variables locales,
%  les différentes variables nécessaires à l'algorithme d'intégration et
%  surtout l'algorithme d'intégration.
%\end{itemize}
%Cette structure est illustrée en figure~\ref{fig:nortonrk:inheritance}.
%
%Cette figure montre également que ces classes sont des classes
%\texttt{template} paramétré par trois arguments~:
%\begin{itemize}
%  \item l'argument \(hypothesis\) désigne l'hypothèse de modélisation~;
%  \item l'argument \(Type\) est le type numérique utilisé~;
%  \item le dernier argument est un booléen qui est généralement faux
%  (choix par défaut). Si il est mis à vrai, par la directive
%  \mkey{UseQt}, les différentes variables sont affectées d'unités, ce
%  qui permet de vérifier à la compilation le respect des unités à la
%  compilation. Cette option est peu usitée car les coefficients des lois
%  de comportements ont souvent des unités exotiques.
%\end{itemize}
%

\subsubsection{Gestion des bornes}

Deux types de bornes sont distinguées dans \mfront{}~:
\begin{itemize}
\item les bornes physiques\footnote{Voir le mot clé
    \mkey{PhysicalBounds}}, qui désignent les plages de valeurs
  acceptables pour une variable donnée. Par exemple, une température
  (en Kelvin) ne peut être négative, une porosité est positive et
  inférieure à 1~;
\item les bornes de validité\footnote{Voir le mot clé \mkey{Bounds}},
  qui désignent les plages de valeurs sur lesquelles la loi de
  comportement a été identifiée.
\end{itemize}

Des bornes peuvent être posées sur~:
\begin{itemize}
  \item les valeurs des variables internes (en début et en fin de pas).
  \item les valeurs des contraintes.
  \item les valeurs des variables externes.
\end{itemize}

\paragraph{Dépassement des bornes de validité} Le traitement d'une
violation d'une borne de validité dépend de l'utilisateur qui peut
choisir parmi trois \og~politiques\fg{}~:
\begin{itemize}
\item ne rien faire~;
\item afficher un message d'avertissement~;
\item arrêter le calcul.
\end{itemize}

La façon de préciser cette politique dépend du code cible (voir
paragraphe~\ref{sec:behaviours:interface}).

% Toutes
%les lois de comportement générées proposent une méthode
%\texttt{set\-Out\-Of\-Bounds\-Policy} qui doit être appelée par les
%interfaces. Cette méthode prend en argument un objet de type énumération
%nommé \tfelb{OutOfBoundsPolicy}{Out\-Of\-Bounds\-Policy} définie dans
%le fichier
%\headerb{TFEL/MaterialLaw/OutOfBoundsPolicy.hxx}{TFEL/\-Material\-Law/\-Out\-Of\-Bounds\-Policy\-.hxx}.
%
%Cette énumération définit trois valeurs valides \texttt{Warning},
%\texttt{Strict} et \texttt{None}.
%
%La politique par défaut est \texttt{None}. Si la politique est fixée à
%\texttt{Warning}, un message est affiché sur la sortie d'erreur en cas
%de violation d'une borne de validité. Si la politique est fixée à
%\texttt{Strict}, une exception du type
%\tfelb{OutOfBoundsException}{Out\-Of\-Bounds\-Exception} est levée.

%
\subsection{Interface aux lois de comportement}
\label{sec:behaviours:interface}

Cette section décrit tout d'abord le rôle des interfaces puis chacune
des deux interfaces disponibles actuellement~:
\begin{itemize}
\item l'interface \umat{} utilisée pour l'adhérence au code aux
  éléments finis \castem{}~\cite{cea_site_2013,helfer_interface_2013}. Cette
  interface est également utilisée par le code d'homogénéisation par
  transformées de \nom{Fourier} rapides
  \tmfft{}~\cite{castelier_specifications_2009,jerome_p._tmfft_2010}~;
\item l'interface \aster{} utilisée pour l'adhérence au code aux
  éléments finis \aster{}~\cite{edf_site_2013,helfer_interface_2013-1}.
\end{itemize}

\subsubsection{Rôle des interfaces aux lois de comportement}

Les interfaces aux lois de comportement ont différentes fonctions~:
\begin{itemize}
\item les lois de comportement choisissent l'une des implantations de
  la loi de comportement en fonction de l'hypothèse de
  modélisation\footnote{Les lois de comportements générées par
    \mfront{} sont représentées par des classes \texttt{template}
    paramétrées par l'hypothèse de modélisation.}~;
  \item assurer la conversion entre la convention utilisée par le code
  appelant pour représenter les tenseurs et la convention utilisée dans
  \TFEL{}~;
\item l'interface doit gérer les lois de comportement orthotropes. Par
  exemple, l'interface peut assurer la rotation des déformations
  totales et de leurs incréments dans le repère propre du matériau
  avant l'appel à la loi de comportement proprement dite et la
  rotation des contraintes dans le repère général après l'appel à la
  loi de comportement.~;
  \item l'interface doit fournir certains éléments nécessaires au calcul
  (matrice d'élasticité) à partir des informations fournies par le
  code~;
\item l'interface doit capter les exceptions \cpp{} et les traduire en
  message d'erreurs adaptés.
\end{itemize}

En fonction des fonctionnalités disponibles ou absentes du code cible,
l'interface peut également~:
\begin{itemize}
  \item assurer le sous-découpage local (au niveau du point de
  \nom{Gauss}) du pas de temps en cas de non convergence~;
  \item gérer le cas des contraintes planes si les lois de
  comportement ne gèrent pas cette hypothèse, ce qui est le cas de la
  plupart des lois générées par \mfront{}~;
  \item permettre l'utilisation des lois écrites pour les petites
  déformations dans un calcul en transformations finies. Pour cela,
  plusieurs pistes semblent intéressantes~:
  \begin{itemize}
  \item le cas des grandes rotations, petites déformations~;
    \item l'utilisation des déformations logarithmiques~;
  \end{itemize}
\end{itemize}

De manière optionnelle, il peut être utile de générer des exemples
d'utilisation de la loi de comportement traitée. Ainsi,
l'interface~\umat{} génère automatiquement un exemple de mise en
données \gibiane{}~\cite{helfer_interface_2013}.

\newpage
\clearpage
\section{Matériaux isotropes à écoulement plastique ou viscoplastique
  incompressible}
\label{sec:mfront:isotropic:analyser}

Nous nous intéressons dans cette section à une famille particulière de
lois de comportement, très utilisée en mécanique, et qui représente la
majorité des comportements de la plate-forme {\pleiades}. Ces lois
s'appliquent à des matériaux {\em isotropes} et décrivent un
comportement {\em plastique} ou {\em viscoplastique}, avec ou sans
{\em écrouissage}, dont les déformations résiduelles sont {\em
  isochores}. Ces particularités permettent d'optimiser les techniques
d'intégration numérique.

Après une présentation des techniques d'intégration adaptées à cette
famille de lois, nous détaillons la syntaxe des quatre analyseurs
{\mfront} qui leurs sont dédiés.

\subsection{Généralités et résolution implicite}

Après une description des particularités des lois de comportement
mécanique traitées dans cette section, les techniques d'intégration
qui leurs sont appliquées dans {\mfront} sont exposées. Celles-ci
s'appuient sur une résolution implicite, une \(\theta\)-méthode
semblable à celles qui sont décrites en section~\ref{sec:Implicite}.

Elles permettent également le calcul de la matrice tangente cohérente.

\subsubsection{Expression des lois de comportement}

Les lois de comportement traitées dans cette section se composent
d'une partie élastique, et de plusieurs mécanismes d'écoulement
plastique ou viscoplastique.

\paragraph{Partition des déformations}
Cette combinaison se traduit par la {\em partition des déformations}:
la déformation totale \(\tepsilonto\) est la somme d'une déformation
élastique \(\tepsilonel\) et d'une déformation inélastique
\(\tepsilonan\):
\begin{subequations}
\label{eq:partition_deformation}
\begin{equation}
  \tepsilonto = \tepsilonel + \tepsilonan
.\end{equation}
Cette dernière se décompose à son tour en plusieurs déformations
relatives aux différents mécanismes d'écoulement, indicés par $i$,
et supposés indépendants:
\begin{equation}
\tepsilonan= \displaystyle\sum_{i}\,\tepsilonan_{i}
.\end{equation}
\end{subequations}

\paragraph{Isotropie}
Les matériaux décrit ici sont supposés {\em isotropes}. Cette
hypothèse sera utilisée pour formuler l'ensemble des mécanismes:
élasticité, écoulements plastiques ou viscoplastiques.

\paragraph{Tenseurs}
Pour simplifier les expressions à venir, il est commode d'introduire
des notations tensorielles, détaillées en
annexe~\ref{sec:oper-tens-dans}~:
\begin{itemize}
\item le produit tensoriel de deux tenseurs $a$, $b$: $a\otimes b$;
\item le produit contracté de deux tenseurs $a$, $b$: $a\colon b$;
\item le tenseur identité d'ordre 2: \(\tenseur{I}\);
\item le tenseur identité d'ordre 4: \(\tenseurq{I}\).
\end{itemize}

\paragraph{Comportement élastique}
Les contraintes \(\tsigma\) se déduisent des déformations élastiques
\(\tepsilonel\) par la loi de \nom{Hooke}. Pour un matériau isotrope,
cette relation peut s'écrire~:
\begin{subequations}
\label{mfront:as:tsigma}
\begin{equation}
\tsigma =
\lambda\,\trace{\tepsilonel}\tenseurq{I}+2\,\mu\,\tepsilonel
\end{equation}
où $\trace{\tepsilonel}$ désigne la trace du tenseur \(\tepsilonel\)
(somme des termes diagonaux). les coefficients de \nom{Lamé}
\(\lambda\) et \(\mu\) du matériau se déduisent du module
d'\nom{Young} et du coefficient de \nom{Poisson}.  Sous forme
tensorielle cette loi s'écrit de manière plus compacte:
\begin{equation}
\tsigma=\tenseurq{D}\colon\tepsilonel
,\quad\text{avec}\quad
\tenseurq{D}=\lambda\,\tenseur{I}\otimes\tenseur{I}+2\,\mu\,\tenseurq{I}
.\end{equation}
La loi est alors résumée par le tenseur élastique $\tenseurq{D}$.
\end{subequations}

\paragraph{Écoulements}
Les écoulements $i$, plastiques ou viscoplastiques, sont les
mécanismes qui créent les déformations inélastiques $\tepsilonan_{i}$
de la partition~\eqref{eq:partition_deformation}. Pour les lois
décrites ici, la direction d'écoulement est supposée proportionnelle
au déviateur des contraintes:
\begin{equation}
\label{eq:direcoul}
\tdepsilonan_{i} \propto \tenseur{s}
,\end{equation}
défini par:
\begin{equation}
\label{eq:mfront:deviateur_contraintes}
\tenseur{s}
=\tsigma-\Frac{1}{3}\,\trace{\tsigma}\,\tenseur{I}
=\tenseurq{K}\colon\tsigma
,\quad\text{avec}\quad
\tenseurq{K}
=\tenseurq{I}-\Frac{1}{3}\tenseur{I}\otimes\tenseur{I}
.\end{equation}
Avec cette hypothèse, les écoulements modélisés vérifient:
\[
\trace{\tepsilonan_{i}} = 0
,\]
c'est-à-dire n'induisent pas de changement de volume: ils sont dits
{\em isochores}.

\paragraph{Contraintes de \nom{Von Mises}}
Les écoulements décrits ici sont supposés dépendre du tenseur des
contraintes $\tsigma$ à travers la norme de son
déviateur~\eqref{eq:mfront:deviateur_contraintes}, appelée contrainte
de \nom{Von Mises}:
\begin{equation}
\label{mfront:sigmaeq} \sigmaeq =
\displaystyle\sqrt{\Frac{3}{2}\,\tenseur{s}\colon\tenseur{s}}
.\end{equation}
La contrainte de \nom{Von Mises} est un des invariants isotropes de la
contrainte. L'hypothèse assure donc le caractère isotrope de la loi
d'écoulement.

\paragraph{Écrouissage}
Pour chaque écoulement $i$, la loi d'écoulement peut également dépendre de son {\em
  écrouissage}, ou {\em déformation inélastique cumulée} $p_i$, défini
par l'équation différentielle:
\begin{subequations}
\label{mfront:as:p}
\begin{equation}
\label{eq.dotpi}
\dot{p}_{i}=\sqrt{\Frac{2}{3}\,\tdepsilonan_{i}\colon\tdepsilonan_{i}}
.\end{equation}
Le facteur \(\pfrac{3}{2}\), introduit par convention dans la
contrainte de \nom{Von Mises}~\eqref{mfront:sigmaeq}, est compensé ici
dans la définition de $\dot{p}_{i}$.  En général, en début de calcul,
au temps $t=t_0$, le matériau est supposé non écroui:
\begin{equation}
p_i(t_0)=0
.\end{equation}
\end{subequations}

Avec l'hypothèse sur la direction d'écoulement~\eqref{eq:direcoul},
la relation~\eqref{eq.dotpi} s'inverse en:
\begin{equation}
\label{mfront:as:epsilonan}
\tdepsilonan_{i} = \dot{p}_{i}\,\tenseur{n}
,\quad\text{avec}\quad
\tenseur{n} = \Frac{3}{2}\Frac{\tenseur{s}}{\sigmaeq}
,\end{equation}
où le tenseur $\tenseur{n}$ est appelé {\em normale} à
l'écoulement. C'est un tenseur déviatorique de norme constante,
qui vérifie:
\begin{equation}
\label{mfront:as:n:2}
\tenseurq{K}\colon\tenseur{n}=\tenseur{n}
,\quad\text{et}\quad
\tenseur{n}\colon\tenseur{n}=\Frac{3}{2}
.\end{equation}

\paragraph{Écoulements supportés}
Il est maintenant possible de formuler les trois types de mécanismes
qui relèvent d'une intégration spécifique dans {\mfront}:
\begin{subequations}
\label{eq:loistype}
\begin{itemize}
\item des écoulements viscoplastiques de la forme suivante, dont
  relève la loi de \nom{Norton}~\cite{chaboche_mecanique_2009}:
\begin{equation}
\label{eq:loistype_a}
  \tdepsilonan_{i} = f_{i}^{\textrm{an}}\paren{\sigmaeq}\,\tenseur{n}
,\quad\text{c'est-à-dire}\quad
\dot{p}_{i} = f_{i}^{\textrm{an}}\paren{\sigmaeq}
;\end{equation}
  \item des écoulements viscoplastiques avec écrouissage de la
  forme suivante, dont relève la loi \nom{Lemaitre}~\cite{chaboche_mecanique_2009}:
\begin{equation}
\label{eq:loistype_b}
  \tdepsilonan_{i} = f_{i}^{\textrm{an}}\paren{\sigmaeq,p_{i}}\,\tenseur{n}
,\quad\text{c'est-à-dire}\quad
\dot{p}_{i} = f_{i}^{\textrm{an}}\paren{\sigmaeq,p_i}
.\end{equation}
\item et des écoulements plastiques qui satisfont une relation du
  type~:
\begin{equation}
\label{eq:loistype_c}
  f_{i}^{\textrm{an}}(\sigmaeq,p_{i})\leq
  0 \quad\quad \dot{p_{i}}\geq 0 \quad\quad
  f_{i}^{\textrm{an}}(\sigmaeq,p_{i})\,\dot{p}_{i} = 0
\end{equation}
Il est classique que la fonction $f_{i}^{\textrm{an}}$ de l'écoulement
plastique ait la dimension d'une contrainte. {\mfront} suppose que
l'utilisateur a respecté cette convention, et divise de ce fait la
fonction $f_{i}^{\textrm{an}}$ par le module d'Young du matériau, pour
que l'ensemble des inéquations~\eqref{eq:loistype_c} aient la
dimension d'une déformation.
\end{itemize}
\end{subequations}
La fonction $f_{i}^{\textrm{an}}$, qui définit l'écoulement, peut
éventuellement faire intervenir des variables externes évoluant
indépendamment de la mécanique (densité de fission, flux de neutrons
rapides, taille de grain, {\it etc.}).

\paragraph{Remarques}
\begin{enumerate}
\item Les écoulements viscoplastiques sans
  écrouissage~\eqref{eq:loistype_a} sont un cas particulier des
  écoulements avec écrouissage~\eqref{eq:loistype_b}. Pour optimiser
  les temps de calculs, il est intéressant de maintenir leur
  distinction lors de l'implantation numérique. Pour les
  développements théoriques, il est préférable de regrouper ces
  écoulements sous leur forme commune~\eqref{eq:loistype_b}.
\item Avec l'hypothèse sur la direction de
  l'écoulement~\eqref{mfront:as:epsilonan}, les écoulements
  viscoplastiques~\eqref{eq:loistype_b} comme les écoulements
  plastiques~\eqref{eq:loistype_c} dépendent de l'écrouissage $p_i$ au
  lieu de dépendre de la déformation $\tepsilonan_{i}$. Cette
  particularité explique l'optimisation possible lors de l'intégration
  des lois~\eqref{eq:loistype}: il suffit d'intégrer des équations
  avec une inconnue scalaire ($p_i$) plutôt que tensorielle
  ($\tepsilonan_{i}$).
\item L'évolution de la déformation
  viscoplastique~\eqref{eq:loistype_b} est donnée par une équation
  différentielle, tandis que la formulation des lois
  plastiques est basée sur une équation~\eqref{eq:loistype_c}: le
  respect de la surface de charge.
\end{enumerate}

\paragraph{Système d'équations}
L'intégration de la loi de comportement proposée par {\mfront}
consiste à charger un point matériel avec la déformation totale
$\tepsilonto$ et à en déduire l'évolution des contraintes
$\tsigma$. Pour cela, la loi s'appuie sur un jeu de variables
internes, conservées en mémoire, et qui évoluent simultanément. Nous
avons décidé de choisir comme variables internes:
\begin{itemize}
  \item la déformation élastique \(\tepsilonel\);
  \item l'écrouissage \(p_{i}\) associé à chaque mécanisme.
\end{itemize}

\begin{subequations}
\label{eq:systemeSpecifique}
Les équations qui permettent de calculer cette évolution sont
rassemblées ici. Il s'agit de:
\begin{enumerate}
\item la partition des déformations~\eqref{eq:partition_deformation}
\begin{equation}
\tepsilonto = \tepsilonel + \tepsilonan
,\quad\text{avec}\quad
\tepsilonan= \displaystyle\sum_{i}\,\tepsilonan_{i}
;\end{equation}

\item la loi d'élasticité~\eqref{mfront:as:tsigma}
\begin{equation}
\tsigma=\tenseurq{D}\colon\tepsilonel
;\end{equation}

\item la décomposition du tenseur des contraintes $\tsigma$ en
  déviateur~\eqref{eq:mfront:deviateur_contraintes} et contrainte de
  \nom{Von Mises}~\eqref{mfront:sigmaeq}:
\begin{equation}
\tenseur{s}= \tenseurq{K}\colon\tsigma
,\quad\text{et}\quad
\sigmaeq=
\displaystyle\sqrt{\Frac{3}{2}\,\tenseur{s}\colon\tenseur{s}}
;\end{equation}

\item la direction des écoulements, orientée suivant la normale à
  l'écoulement~\eqref{mfront:as:epsilonan}:
\begin{equation}
\tdepsilonan_{i} = \dot{p}_{i}\,\tenseur{n}
,\quad\text{avec}\quad
\tenseur{n}= \Frac{3}{2}\Frac{\tenseur{s}}{\sigmaeq}
;\end{equation}

\item les écoulements viscoplastiques~\eqref{eq:loistype_b} ou
  plastiques~\eqref{eq:loistype_c}:
\begin{align}
\dot{p}_{i} &= f_{i}^{\textrm{an}}\paren{\sigmaeq,p_i}
,\\
f_{i}^{\textrm{an}}(\sigmaeq,p_{i})\leq
0,\quad&\dot{p_{i}}\geq 0,\quad
f_{i}^{\textrm{an}}(\sigmaeq,p_{i})\,\dot{p}_{i} = 0
.\end{align}
\end{enumerate}
\end{subequations}

\subsubsection{Méthode d'intégration numérique}
\label{sec:mfront:isotropic-solver}

Dans {\mfront}, la loi de comportement~\eqref{eq:systemeSpecifique}
est intégrée par une méthode implicite, une \(\theta\)-méthode,
semblable à celle présentée au paragraphe~\ref{sec:Implicite}, à ceci
près que le système d'équations peut ici être réduit à une unique
équation scalaire~\cite{proix_integration_2012-1}.

\paragraph{Incrément de temps}
D'après ce qui précède, les variables internes sont la déformation
élastique $\tepsilonel$ et l'écrouissage $p_i$ associé à chaque
mécanisme. Elles sont connues au temps $t$, en début de pas de temps,
et pour une déformation totale $\debutpas{\tepsilonto}$. L'intégration
consiste à calculer leurs nouvelles valeurs induites par un incrément
de déformation totale $\Delta\tepsilonto$ en fin d'un pas de temps
$\Delta t$. Leurs incréments respectifs seront notés
$\Delta\tepsilonel$ et $\Delta p_i$.

La notion de \(\theta\)-méthode, employée ici, consiste à s'intéresser
également aux valeurs de certaines variables au temps intermédiaire
$t+\theta\,\Delta t$, où le paramètre $\theta$, déterminé à l'avance
est compris entre 0 et 1. Une première approximation consiste à
supposer les variables internes linéaires sur le pas de temps
$\Delta t$, donc à poser:
\begin{subequations}
\label{eq:applin}
\begin{align}
\milieupas{\tepsilonel}
&\approx\debutpas{\tepsilonel}+\theta\,\Delta\tepsilonel
,\\
\milieupas{p_i}&\approx
\debutpas{p_i}+\theta\,\Delta p_i
.\end{align}
\end{subequations}

\paragraph{Contraintes}
D'après la loi d'élasticité~\eqref{mfront:as:tsigma} et sa propre
définition~\eqref{eq:mfront:deviateur_contraintes}, le déviateur des
contraintes vaut:
\begin{equation*}
\tenseur{s}=2\,\mu\,\tenseurq{K}\colon\tepsilonel
,\end{equation*}
et d'après les approximations~\eqref{eq:applin}, sa valeur au temps
intermédiaire vaut:
\begin{equation}
\milieupas{\tenseur{s}}=2\,\mu\,\tenseurq{K}\colon
\left(\debutpas{\tepsilonel}+\theta\,\Delta\tepsilonel
\right)
,\label{eq:stthDt}
\end{equation}
et la contrainte de \nom{Von Mises}~\eqref{mfront:sigmaeq} et la
normale à l'écoulement~\eqref{mfront:as:epsilonan} pour cette même
date s'en déduisent:
\begin{equation}
\milieupas{\sigmaeq}=
\displaystyle\sqrt{\Frac{3}{2}\,\milieupas{\tenseur{s}}\colon\milieupas{\tenseur{s}}}
,\quad\text{et}\quad
\milieupas{\tenseur{n}}
=\Frac{3}{2}\Frac{\milieupas{\tenseur{s}}}{\milieupas{\sigmaeq}}
.\label{eq:ntthDt}
\end{equation}

\paragraph{Partition des déformations}
La partition des déformations~\eqref{eq:partition_deformation}, et la
direction des écoulements~\eqref{mfront:as:epsilonan} s'intègrent en:
\begin{equation}
\Delta\tepsilonto
=
\Delta\tepsilonel+
\sum_i\int_t^{t+\Delta t}
\dot{p}_{i}\,\tenseur{n}\,\dtot t
\approx
\Delta\tepsilonel+
\milieupas{\tenseur{n}}\,
\sum_i\Delta p_{i}
.\label{eq:partdef2}
\end{equation}
Cette approximation s'appuie sur la valeur de la normale $\tenseur{n}$
au temps intermédiaire. Avec les relations
précédentes~\eqref{eq:stthDt} et \eqref{eq:ntthDt}, et le type
déviateur~\eqref{mfront:as:n:2} de $\tenseur{n}$, cela permet
d'écrire:
\begin{equation*}
\milieupas{\tenseur{n}}
=\Frac{3\,\mu}{\milieupas{\sigmaeq}}
\left[
\tenseurq{K}\colon
\left(\debutpas{\tepsilonel}+\theta\,\Delta\tepsilonto
\right)
-\milieupas{\tenseur{n}}\,
\theta\,\sum_i\Delta p_{i}
\right]
,\end{equation*}
c'est-à-dire:
\begin{equation}
\label{eq:signB}
\paren{\milieupas{\sigmaeq}+3\mu\theta\sum_i\Delta p_{i}}
\milieupas{\tenseur{n}}
=3\mu\tenseur{B}
,\quad\text{avec}\quad
\tenseur{B}=\tenseurq{K}\colon
\left(\debutpas{\tepsilonel}+\theta\,\Delta\tepsilonto
\right)
.\end{equation}

La norme de cette expression~\eqref{eq:signB}, compte tenu de la
norme~\eqref{mfront:as:n:2} de $\milieupas{\tenseur{n}}$, donne
finalement pour la contrainte équivalente:
\begin{equation}
  \label{eq:IsotropicBehaviour:contrainte_equivalente}
  \milieupas{\sigmaeq}
=\mu\sqrt{6\,\tenseur{B}\colon\tenseur{B}}-3\,\mu\,\theta\,\sum_{i}\,\Delta p_{i}
.\end{equation}
Cette valeur, réinjectée dans la relation~\eqref{eq:signB}, donne la normale:
\begin{equation}
  \label{eq:IsotropicBehaviour:normale}
  \milieupas{\tenseur{n}}=\Frac{3}{\sqrt{6\,\tenseur{B}\colon\tenseur{B}}}\tenseur{B}
.\end{equation}

L'équation~\eqref{eq:IsotropicBehaviour:contrainte_equivalente} va
permettre de former, à partir des équations d'écoulements, un système
non linéaire d'équations dont les inconnues sont les incréments
$\Delta p_{i}$. Il faut maintenant construire ce système en intégrant
les écoulements viscoplastiques, puis plastiques.

\begin{subequations}
\label{eq:fpis}
\paragraph{Écoulement viscoplastique}
L'intégrale sur le pas de temps $\Delta t$ de l'écoulement
viscoplastique~\eqref{eq:loistype_b} s'écrit:
\begin{equation*}
\Delta p_{i}
=\int_t^{t+\Delta t}f_{i}^{\textrm{an}}
\paren{\sigmaeq,p_i}\,\dtot t
\approx\Delta t\,
f_{i}^{\textrm{an}}\paren{\milieupas{\sigmaeq},\milieupas{p_{i}}}
,\end{equation*}
où l'intégrale est approchée par la valeur de $f_{i}^{\textrm{an}}$ au
temps intermédiaire. Avec cette approximation, vérifier l'écoulement
consiste à annuler une fonction à trois paramètres:
\begin{equation}
f_{p_{i}}\paren{\Delta p_{i},\milieupas{\sigmaeq},\milieupas{p_{i}}}
=\Delta p_{i}
-\Delta t\,f^{\textrm{an}}_{i}\paren{\milieupas{\sigmaeq},\milieupas{p_{i}}}
.\label{eq:fpivis}
\end{equation}

\paragraph{Écoulement plastique} L'écoulement plastique consiste à
annuler une fonction de même forme, dépendant des mêmes trois
paramètres:
\begin{equation}
f_{p_{i}}\paren{\Delta p_{i},\milieupas{\sigmaeq},\milieupas{p_{i}}}
=\left\{
\begin{aligned}
  f^{\textrm{an}}_{i}\paren{\milieupas{\sigmaeq},\milieupas{p_{i}}}
  &\quad\text{en cas de charge plastique,}\\
  \Delta p_{i}&\quad\text{sinon.}
\end{aligned}
\right.
\label{eq:fpiplas}
\end{equation}
La condition de charge plastique, qui augmente l'écrouissage $p_i$,
doit respecter les conditions d'écoulement
plastique~\eqref{eq:loistype_c}. Dans la pratique, les incréments
$\Delta p_i$, inconnues d'un système non linéaire, sont calculés de
manière itérative. Chaque nouvelle itération propose de nouvelles
valeurs de ces incréments. Pour ces nouvelles valeurs, la condition de
charge plastique est activée si l'une des conditions suivante est
vérifiée:
\begin{enumerate}
\item \(f^{\textrm{an}}_{i}\paren{\milieupas{\sigmaeq},\milieupas{p_{i}}} >
  \varepsilon\) et \(\Delta p_{i} \geq 0\),
\item \(\Delta p_{i} > \varepsilon\),
\end{enumerate}
où \(\varepsilon\) est un paramètre numérique, qui stabilise les
itérations.

\end{subequations}

\paragraph{Système non linéaire}
Les équations d'écoulement~\eqref{eq:fpis}:
\begin{equation*}
f_{p_{i}}\paren{\Delta p_{i},\milieupas{\sigmaeq},\milieupas{p_{i}}}
=0,\quad\forall\ p_i
,\end{equation*}
modifiées avec les équations~\eqref{eq:applin} et
\eqref{eq:IsotropicBehaviour:contrainte_equivalente}, forment un
système non linéaire d'équations, dont les inconnues sont les
incréments $\Delta
p_i$:
\begin{equation}
f_{p_{i}}\paren{\Delta p_{i},
\mu\sqrt{6\,\tenseur{B}\colon\tenseur{B}}-3\,\mu\,\theta\,\sum_{j}\,\Delta p_{j},
\debutpas{p_i}+\theta\,\Delta p_i}
=0,\quad\forall\ p_i
.\label{eq:sysnonlinDpi}
\end{equation}
Pour le résoudre, une méthode de \nom{Newton-Raphson} est utilisée.
Cette méthode, discutée en détail au paragraphe~\ref{sec:NR},
nécessite les dérivées partielles des fonctions $f_{p_{i}}$ par
rapport aux incréments $\Delta p_i$. Celles-ci s'obtiennent par une
dérivation composée des $f_{p_i}$:
\begin{equation*}
\deriv{f_{p_{i}}}{\Delta p_j}=
\deriv{f_{p_{i}}}{\Delta p_i}
\,\deriv{\Delta p_i}{\Delta p_j}
+\deriv{f_{p_{i}}}{\milieupas{\sigmaeq}}
\,\deriv{\milieupas{\sigmaeq}}{\Delta p_j}
+\deriv{f_{p_{i}}}{\milieupas{p_{i}}}
\,\deriv{\milieupas{p_{i}}}{\Delta p_j}
,\end{equation*}
avec, d'après les équations~\eqref{eq:applin} et
\eqref{eq:IsotropicBehaviour:contrainte_equivalente}:
\begin{equation*}
\deriv{\Delta p_i}{\Delta p_j}=\delta_{ij}
,\quad
\deriv{\milieupas{\sigmaeq}}{\Delta p_j}=-3\,\mu\,\theta
,\quad\text{et}\quad
\deriv{\milieupas{p_{i}}}{\Delta p_j}=\theta\,\delta_{ij}
,\quad\text{avec}\quad
\delta_{ij}=
\left\{
\begin{array}{ll}
1&\text{si $i=j$,}\\0&\text{sinon.}
\end{array}
\right.
\end{equation*}

Appliquée aux écoulements, ces relations différentielles donnent:
\begin{itemize}
\item pour les écoulement viscoplastiques~\eqref{eq:fpivis}:
\[
\deriv{f_{p_{i}}}{\Delta p_{j}}=
\left\{
\begin{aligned}
1-3\,\mu\,\theta\,\Delta t\,\left[\deriv{f^{\textrm{an}}_{i}}{\sigmaeq}\paren{\milieupas{\sigmaeq},\milieupas{p_{i}}}\right]-\theta\,\Delta t\,\left[\deriv{f^{\textrm{an}}_{i}}{p_{i}}\paren{\milieupas{\sigmaeq},\milieupas{p_{i}}}\right] & \text{ si } i=j, \\
-3\,\mu\,\theta\,\Delta t\,\left[\deriv{f^{\textrm{an}}_{i}}{\sigmaeq}\paren{\milieupas{\sigmaeq},\milieupas{p_{i}}}\right] & \text{ si } i\neq j;
\end{aligned}
\right.
\]

\item pour les écoulement plastiques~\eqref{eq:fpiplas}, en cas de
  décharge plastique:
\[
\deriv{f_{p_{i}}}{\Delta p_{j}}= \left\{
\begin{aligned}
  \theta\,\paren{3\,\mu\,\left[\deriv{f^{\textrm{an}}_{i}}{\sigmaeq}\paren{\milieupas{\sigmaeq},\milieupas{p_{i}}}\right]+
    \left[\deriv{f^{\textrm{an}}_{i}}{p_{i}}\paren{\milieupas{\sigmaeq},\milieupas{p_{i}}}\right]} & \text{ si } i=j, \\
  -3\,\mu\,\theta\,\deriv{f^{\textrm{an}}_{i}}{\sigmaeq}\paren{\milieupas{\sigmaeq},\milieupas{p_{i}}}
  & \text{ si } i\neq j.
\end{aligned}
\right.
\]
et en absence de charge plastique:
\[
\deriv{f_{p_{i}}}{\Delta p_{j}}= \left\{
\begin{aligned}
  1 & \text{ si } i=j, \\
  0 & \text{ si } i\neq j.
\end{aligned}
\right.
\]
\end{itemize}

\paragraph{Critère d'arrêt} L'algorithme de \nom{Newton} s'arrête quand
la différence entre deux estimations des incréments des déformations
viscoplastiques cumulées est inférieur à un certain critère
\(\varepsilon\)~:
\begin{equation}
  \label{eq:mfront:isotropic:convergence_criteria}
  \Frac{1}{N}\displaystyle\sum_{i=1}^{N}\left|\Delta_{n}\Delta
  p_{i}\right|<\varepsilon
\end{equation}
où\(\Delta_{n}\Delta p_{i}\) désigne la différence entre l'estimation
de l'incrément de la déformation viscoplastique cumulée \(\Delta\,
p_{i}\) à l'étape \(n+1\) et de l'estimation à l'étape \(n\).

L'algorithme échoue si le nombre d'itérations dépasse une
borne maximale.

\paragraph{Étapes}
Finalement les différentes étapes de l'intégration de la loi de
comportement mécanique par la $\theta$-méthode sont les suivantes:
\begin{enumerate}
\item calcul du tenseur $\tenseur{B}$~\eqref{eq:signB}:
\begin{equation*}
\tenseur{B}=\tenseurq{K}\colon
\left(\debutpas{\tepsilonel}+\theta\,\Delta\tepsilonto
\right)
~;\end{equation*}

\item calcul des incréments $\Delta p_i$, par résolution du système
  non linéaire~\eqref{eq:sysnonlinDpi}:
\begin{equation*}
f_{p_{i}}\paren{\Delta p_{i},
\mu\sqrt{6\,\tenseur{B}\colon\tenseur{B}}-3\,\mu\,\theta\,\sum_{j}\,\Delta p_{j},
\debutpas{p_i}+\theta\,\Delta p_i}
=0,\quad\forall\ p_i
,\end{equation*}
en utilisant la méthode de \nom{Newton-Raphson}, et calcul des
nouvelles valeurs d'écrouissage:
\begin{equation*}
\finpas{p_i}=\debutpas{p_i}+\Delta p_i
~;\end{equation*}

\item calcul de l'incrément de déformation élastique, grâce à la
  partition des déformations~\eqref{eq:partdef2}, de la
  normale~\eqref{eq:IsotropicBehaviour:normale}:
\begin{equation*}
\Delta\tepsilonel
=
\Delta\tepsilonto
-\milieupas{\tenseur{n}}\,
\sum_i\Delta p_{i}
,\quad\text{avec}\quad
\milieupas{\tenseur{n}}
=\Frac{3}{\sqrt{6\,\tenseur{B}\colon\tenseur{B}}}\tenseur{B}
,\end{equation*}
et de la nouvelle valeur de déformation élastique:
\begin{equation*}
\finpas{\tepsilonel}=\debutpas{\tepsilonel}+\Delta\tepsilonel
~;\end{equation*}
\item calcul de la contrainte en fin de pas de temps, par la loi
  d'élasticité~\eqref{mfront:as:tsigma}:
\begin{equation*}
\finpas{\tsigma}=\tenseurq{D}\colon\finpas{\tepsilonel}
.\end{equation*}
\end{enumerate}

\subsubsection{Matrice tangente cohérente}

La méthode de résolution implicite permet le calcul de la matrice
tangente cohérente. Cette matrice tangente cohérente permet une
convergence quadratique du calcul de structure, ce qui peut nettement
accélérer les calculs.

\paragraph{Définition}
L'intégration numérique décrite au paragraphe précédent, permet de
calculer l'évolution des con\-train\-tes $\Delta\tsigma$ induite par un
incrément de déformation $\Delta\tepsilonto$ sur un pas de temps
$\Delta t$. La matrice tangente cohérente est définie comme le
tenseur:
\begin{equation}
\tenseurq{L}^{\mathit{tc}}
=\deriv{\Delta\tsigma}{\Delta\tepsilonto}
.\label{eq:defTangent}
\end{equation}
La loi d'élasticité~\eqref{mfront:as:tsigma} et la partition des
déformations~\eqref{eq:partdef2} permettent d'exprimer l'incrément de
contrainte:
\begin{equation*}
\Delta\tsigma
=\tenseurq{D}\colon\Delta\tepsilonel
=\tenseurq{D}\colon\paren{\Delta\tepsilonto
-\milieupas{\tenseur{n}}\,
\sum_i\Delta p_{i}
}
\end{equation*}
et d'en déduire par dérivation~\eqref{eq:defTangent} la matrice
tangente cohérente~:
\begin{equation}
\tenseurq{L}^{\mathit{tc}}
=\tenseurq{D}-\tenseurq{D}\colon\paren{\milieupas{\tenseur{n}}\otimes
\sum_i\deriv{\Delta p_{i}}{\Delta\tepsilonto}
+\sum_i\Delta p_{i}\deriv{\milieupas{\tenseur{n}}}{\Delta\tepsilonto}}
.\label{eq:tang1}
\end{equation}
Il faut maintenant calculer les dérivées de la normale
$\milieupas{\tenseur{n}}$ et des incréments $\Delta p_i$ pour chaque écoulement.

\paragraph{Normale}
La dérivée du tenseur~\eqref{eq:signB} $\tenseur{B}$ et de la
normale~\eqref{eq:IsotropicBehaviour:normale} valent respectivement:
\begin{subequations}
\begin{align}
\deriv{\tenseur{B}}{\Delta\tepsilonto}
&=\theta\,\tenseurq{K}
,\label{eq:dBdDeto}
\\
\deriv{\milieupas{\tenseur{n}}}{\Delta\tepsilonto}
&=\Frac{3}{\sqrt{6\,\tenseur{B}\colon\tenseur{B}}}
\,\left(
\deriv{\tenseur{B}}{\Delta\tepsilonto}
-\frac{2}{3}
\deriv{\tenseur{B}}{\Delta\tepsilonto}\colon
\milieupas{\tenseur{n}}\otimes\milieupas{\tenseur{n}}
\right)
,\notag\\
&=\Frac{3\,\theta}{\sqrt{6\,\tenseur{B}\colon\tenseur{B}}}
\,\left(\tenseurq{K}
-\frac{2}{3}\milieupas{\tenseur{n}}\otimes\milieupas{\tenseur{n}}
\right)
.\label{eq:dndDeto}\end{align}
\end{subequations}
Ce calcul utilise la propriété~\eqref{mfront:as:n:2} du tenseur
$\milieupas{\tenseur{n}}$.

\paragraph{Incréments}
Les incréments d'écrouissage $\Delta p_i$ sont issus de la résolution
du système non linéaire, formé par les équations
d'écoulement~\eqref{eq:fpis}:
\begin{equation*}
f_{p_{i}}\paren{\Delta p_{i},\milieupas{\sigmaeq},\milieupas{p_{i}}}
=0,\quad\forall\ p_i
.\end{equation*}
La dérivation de ces équations conduit à:
\begin{equation}
\label{eq:dfpidDetoa}
\deriv{f_{p_{i}}}{\Delta p_i}
\,\deriv{\Delta p_i}{\Delta\tepsilonto}
+\deriv{f_{p_{i}}}{\milieupas{\sigmaeq}}
\,\deriv{\milieupas{\sigmaeq}}{\Delta\tepsilonto}
+\deriv{f_{p_{i}}}{\milieupas{p_{i}}}
\,\deriv{\milieupas{p_{i}}}{\Delta\tepsilonto}
=0,\quad\forall\ p_i
.\end{equation}
La dérivation des équations~\eqref{eq:applin} et
\eqref{eq:IsotropicBehaviour:contrainte_equivalente}, en utilisant
les relations~\eqref{mfront:as:n:2} et \eqref{eq:dBdDeto}:
\begin{align*}
\deriv{\milieupas{\sigmaeq}}{\Delta\tepsilonto}
&=2\,\mu\,\theta\,\milieupas{\tenseur{n}}
-3\,\mu\,\theta\,\sum_j\,\deriv{\Delta p_j}{\Delta\tepsilonto}
,\\
\deriv{\milieupas{p_{i}}}{\Delta\tepsilonto}
&=\theta\,\deriv{\Delta p_i}{\Delta\tepsilonto}
,\end{align*}
permet de modifier les équations~\eqref{eq:dfpidDetoa}:
\begin{equation}
\label{eq:dfpidDetob}
\left(
\deriv{f_{p_{i}}}{\Delta p_i}
+\theta\,\deriv{f_{p_{i}}}{\milieupas{p_{i}}}
\right)
\,\deriv{\Delta p_i}{\Delta\tepsilonto}
-3\,\mu\,\theta\,\deriv{f_{p_{i}}}{\milieupas{\sigmaeq}}
\,\sum_j\,\deriv{\Delta p_j}{\Delta\tepsilonto}
=-\deriv{f_{p_{i}}}{\milieupas{\sigmaeq}}
\,2\,\mu\,\theta\,\milieupas{\tenseur{n}}
,\quad\forall\ p_i
.\end{equation}
C'est un jeu d'équations linéaires, dont les inconnues sont les
tenseurs dérivées recherchés $\deriv{\Delta p_i}{\Delta\tepsilonto}$.

\paragraph{Solutions scalaires}
Le produit contracté des équations~\eqref{eq:dfpidDetob} par
$\milieupas{\tenseur{n}}$ donne un système d'équations:
\begin{subequations}
\label{eq:dfpidDetoc}
\begin{equation}
\left(
\deriv{f_{p_{i}}}{\Delta p_i}
+\theta\,\deriv{f_{p_{i}}}{\milieupas{p_{i}}}
\right)\,d_i
-3\,\mu\,\theta\,\deriv{f_{p_{i}}}{\milieupas{\sigmaeq}}
\,\sum_j\,d_j
=-\deriv{f_{p_{i}}}{\milieupas{\sigmaeq}}
\,3\,\mu\,\theta
,\quad\forall\ p_i
.\end{equation}
avec pour inconnues, les produits contractés:
\begin{equation}
d_i=\deriv{\Delta p_i}{\Delta\tepsilonto}\colon\milieupas{\tenseur{n}}
.\end{equation}
\end{subequations}

Le produit contracté des équations~\eqref{eq:dfpidDetob} par tout
tenseur $\tenseur{x}$ orthogonal à $\milieupas{\tenseur{n}}$ donne un
système d'équations:
\begin{equation*}
\left(
\deriv{f_{p_{i}}}{\Delta p_i}
+\theta\,\deriv{f_{p_{i}}}{\milieupas{p_{i}}}
\right)\,\deriv{\Delta p_i}{\Delta\tepsilonto}\colon\tenseur{x}
-3\,\mu\,\theta\,\deriv{f_{p_{i}}}{\milieupas{\sigmaeq}}
\,\sum_j\,\deriv{\Delta p_j}{\Delta\tepsilonto}\colon\tenseur{x}
=0
,\quad\forall\ p_i
,\end{equation*}
de second membre nul. Ses inconnues sont donc nulles:
\begin{equation}
\deriv{\Delta p_i}{\Delta\tepsilonto}\colon\tenseur{x}
=0
,\quad\forall\ p_i
,\end{equation}
ce qui montre que les tenseurs $\deriv{\Delta p_i}{\Delta\tepsilonto}$
sont tous colinéaires au tenseur $\milieupas{\tenseur{n}}$. Ils
s'écrivent donc:
\begin{equation}
\deriv{\Delta p_i}{\Delta\tepsilonto}
=\Frac{\deriv{\Delta p_i}{\Delta\tepsilonto}:\milieupas{\tenseur{n}}}
{\milieupas{\tenseur{n}}\colon\milieupas{\tenseur{n}}}
\,\milieupas{\tenseur{n}}
=\frac{2\,d_i}{3}\,\milieupas{\tenseur{n}}
,\label{eq:dDpdDeto}
\end{equation}
où les composantes $d_i$ sont les solutions du
système~\eqref{eq:dfpidDetoc}.

\paragraph{Matrice tangente cohérente}
Finalement, la matrice tangente cohérente peut être reconstituée à
partir des expressions~\eqref{eq:tang1}, \eqref{eq:dndDeto} et
\eqref{eq:dDpdDeto}, et en notant que:
\begin{equation*}
\tenseurq{D}\colon\tenseurq{K}=2\,\mu\,\tenseurq{K}
,\quad\text\quad
\tenseurq{D}\colon\milieupas{\tenseur{n}}=2\,\mu\,\milieupas{\tenseur{n}}
.\end{equation*}
Elle vaut ainsi:
\begin{equation}
\tenseurq{L}^{\mathit{tc}}
=\tenseurq{D}-2\,\mu\,\left(
\frac{3\,\theta\,\sum_i\Delta p_i}
{\sqrt{6\,\tenseur{B}\colon\tenseur{B}}}
\,\tenseurq{K}
+\left(\frac{2\,\sum_id_i}{3}
-\frac{2\,\theta\,\sum_i\Delta p_i}
{\sqrt{6\,\tenseur{B}\colon\tenseur{B}}}
\right)
\,\milieupas{\tenseur{n}}\otimes\milieupas{\tenseur{n}}
\right)
,\label{eq:tang2}
\end{equation}
avec les composantes $d_i$ solutions du système~\eqref{eq:dfpidDetoc}.

\subsubsection{Expression de la matrice tangente cohérente pour un
  mécanisme unique}

Dans le cas où le système se limite à un mécanisme unique,
l'expression de la matrice tangente cohérente se simplifie.

\paragraph{Écoulement viscoplastique}
Dans le cas d'un seul écoulement viscoplastique~\eqref{eq:fpivis},
d'incrément $\Delta p$, le système~\eqref{eq:dfpidDetoc} se réduit à
une équation d'inconnue $d$:
\begin{equation*}
  \left(
    1
    -\Delta t\,\theta\,\deriv{f^{\textrm{an}}}{p}
    +3\,\mu\,\Delta t\,\theta\,\deriv{f^{\textrm{an}}}{\sigmaeq}
  \right)\,d
  =\deriv{f^{\textrm{an}}}{\sigmaeq}
  \,3\,\mu\,\Delta t\,\theta
  .\end{equation*}
L'expression~\eqref{eq:tang2} de la matrice tangente cohérente se simplifie
alors en:
\begin{equation}
\label{eq:mfront:Dt:visco}
\tenseurq{L}^{\mathit{tc}}
=\tenseurq{D}-2\,\mu\,\theta\,\left(
\Frac{3\,\Delta p}
{\sqrt{6\,\tenseur{B}\colon\tenseur{B}}}
\,\tenseurq{K}
+\left(
\Frac{2\,\mu\,\Delta t
\,\deriv{f^{\textrm{an}}}{\sigmaeq}}
{1+\theta\,\Delta t\,\left(
3\,\mu\,\deriv{f^{\textrm{an}}}{\sigmaeq}
-\deriv{f^{\textrm{an}}}{p}\right)
}
-\Frac{2\,\Delta p}
{\sqrt{6\,\tenseur{B}\colon\tenseur{B}}}
\right)
\,\milieupas{\tenseur{n}}\otimes\milieupas{\tenseur{n}}
\right)
\end{equation}

\paragraph{Écoulement plastique}
Dans le cas d'un seul écoulement viscoplastique,
d'incrément $\Delta p$, le système~\eqref{eq:dfpidDetoc} se réduit à
une équation d'inconnue $d$. Si l'écoulement plastique est
activé~\eqref{eq:fpiplas}, cette équation devient:
\begin{equation*}
  \left(
    \theta\,\deriv{f^{\textrm{an}}}{p}
    -3\,\mu\,\theta\,\deriv{f^{\textrm{an}}}{\sigmaeq}
  \right)\,d
  =-\deriv{f^{\textrm{an}}}{\sigmaeq}
  \,3\,\mu\,\theta
  .\end{equation*}
L'expression~\eqref{eq:tang2} de la matrice tangente cohérente se simplifie
alors en:
\begin{equation}
\label{eq:mfront:Dt:plasti}
\tenseurq{L}^{\mathit{tc}}
=\tenseurq{D}-2\,\mu\,\theta\,\left(
\Frac{3\,\Delta p}
{\sqrt{6\,\tenseur{B}\colon\tenseur{B}}}
\,\tenseurq{K}
+\left(
\Frac{2\,\mu\,\deriv{f^{\textrm{an}}}{\sigmaeq}}
{3\,\mu\,\theta\,\deriv{f^{\textrm{an}}}{\sigmaeq}
-\theta\,\deriv{f^{\textrm{an}}}{p}
}
-\Frac{2\,\Delta p}
{\sqrt{6\,\tenseur{B}\colon\tenseur{B}}}
\right)
\,\milieupas{\tenseur{n}}\otimes\milieupas{\tenseur{n}}
\right)
\end{equation}

Lorsque l'écoulement plastique n'est pas activé, l'équation d'écoulement~\eqref{eq:fpiplas} et
le système~\eqref{eq:dfpidDetoc} se réduisent respectivement à:
\begin{equation*}
\Delta p=0
,\quad\text{et}\quad
d=0
,\end{equation*}
et la matrice tangente cohérente au tenseur d'élasticité:
\begin{equation*}
\tenseurq{L}^{\mathit{tc}}
=\tenseurq{D}
.\end{equation*}

\subsubsection{Matrice tangente}

La matrice tangente définie par la relation en vitesse~:
\[
\tdsigma = \tenseurq{L}\,\colon\,\tdepsilonto
\]

Sans chercher à justifier cette affirmation, la matrice tangente peut
être calculée comme la limite de la matrice tangente cohérente lorsque
le pas de temps \(\Delta \,t\) et l'incrément de déformation plastique
cumulée \(\Delta \,p\) tendent vers \(0\) et en posant \(\theta=1\)~:
\[
\tenseurq{L} = \lim_{
  \tiny
  \begin{aligned}
    \Delta\,t\,&\rightarrow\,0\\
    \Delta\,p\,&\rightarrow\,0\\
    \theta&=1
  \end{aligned}
}\tenseurq{L}^{\mathit{tc}}
\]

Les expressions précédentes peuvent donc être réutilisées en prenant
garde à l'instant où l'on veut calculer la matrice tangente (en début
ou fin de pas de temps).

\subsubsection{Expression de la matrice tangente pour un mécanisme
  unique}

Dans le cas où le système se limite à un mécanisme unique,
l'expression de la matrice tangente cohérente se simplifie.

\paragraph{Écoulement viscoplastique}
Dans le cas d'un seul écoulement viscoplastique,
l'expression~\eqref{eq:mfront:Dt:visco} de la matrice tangente
cohérente se réduit à la matrice d'élasticité.

\paragraph{Écoulement plastique}
Dans le cas d'un écoulement plastique, deux cas se présentent suivant
que l'on soit en charge ou pas.

En cas de charge plastique, l'expression~\eqref{eq:mfront:Dt:plasti}
de la matrice tangente cohérente conduit à~:
\[
\tenseurq{L}
=\tenseurq{D}-\Frac{4\,\mu^{2}\,\deriv{f^{\textrm{an}}}{\sigmaeq}}
{3\,\mu\,\deriv{f^{\textrm{an}}}{\sigmaeq}
  -\deriv{f^{\textrm{an}}}{p}
}
\tenseur{n}\otimes\tenseur{n}
\]

Dans le domaine élastique ou en cas de décharge, la matrice tangente
se réduit à la matrice d'élasticité.

\subsection{Utilisation des analyseurs spécifiques}

Nous présentons dans ce paragraphe les analyseurs \mfront{} dédiés aux
lois de comportement plastique et viscoplastique incompressible des
matériaux isotropes. Ils sont au nombre de \(4\)~:
\begin{itemize}
  \item l'analyseur \texttt{IsotropicMisesCreep} gère
  exclusivement les lois de comportement viscoplastique isotrope de la
  forme~:
  \[
  \dot{p} = f\paren{\sigmaeq}
  \]
  \item l'analyseur \texttt{IsotropicStrainHardeningMisesCreep}
  gère exclusivement les lois de comportement viscoplastique
  isotrope de la forme~:
  \[
  \dot{p} = f\paren{\sigmaeq,p}
  \]
  \item l'analyseur \texttt{IsotropicPlasticMisesFlow} gère
  exclusivement les lois de comportement plastique isotrope de la
  forme~:
  \[
  f(\sigmaeq,p)\leq 0 \quad\quad \dot{p}\geq 0
  \quad\quad f(\sigmaeq,p)\,\dot{p} = 0
  \]
  \item l'analyseur \texttt{MultipleIsotropicMisesFlows} gère
  une combinaison arbitraire d'écoulements des trois types précédents.
  Les différents écoulements sont supposés non couplés.
\end{itemize}

\subsubsection{La directive \texttt{@FlowRule}} La directive \mkey{FlowRule}
permet de définir un écoulement.

\paragraph{Cas des analyseurs \texttt{Isotropic\-Mises\-Creep},
  \texttt{Isotropic\-Strain\-Hardening\-Mises\-Creep} et
  \texttt{Isotropic\-Plastic\-Mises\-Flow}}

Pour les analyseurs \texttt{Isotropic\-Mises\-Creep},
\texttt{Isotropic\-Strain\-Hardening\-Mises\-Creep} et
\texttt{Isotropic\-Plastic\-Mises\-Flow}, une variable interne nommée
\(p\) et représentant la déformation inélastique cumulée est
automatiquement définie. Aucune autre variable interne ne peut être
définie.

L'évolution de cette variable est introduite par la directive
\mkey{FlowRule}. Cette directive est suivie d'un bloc définissant
l'écoulement. Ce bloc doit renseigner la valeur de la fonction
\texttt{f}, dont la définition a été donnée plus haut en fonction de
l'écoulement traité, et sa dérivée par rapport à la contrainte
équivalente \texttt{df\textunderscore{}dseq}. Pour les analyseurs
\texttt{Isotropic\-Strain\-Hardening\-Mises\-Creep} et
\texttt{Isotropic\-Plastic\-Mises\-Flow} il est également nécessaire
de donner la dérivée de \texttt{f} par rapport à la déformation
cumulée \texttt{df\textunderscore{}dp}.

Dans le bloc suivant la directive \mkey{FlowRule}, les variables
\texttt{f}, \texttt{df\textunderscore{}dseq} et éventuellement
\texttt{df\textunderscore{}dp} sont automatiquement définies. La
contrainte équivalente actualisée en \(t+\theta\,\Delta\,t\) est
accessible par la variable \texttt{seq}. Pour les analyseurs
\texttt{Isotropic\-Strain\-Hardening\-Mises\-Creep} et
\texttt{Isotropic\-Plastic\-Mises\-Flow}. Si nécessaire, la
déformation équivalente actualisée en \(t+\theta\,\Delta\,t\) est
accessible par la variable \texttt{p}.

\paragraph{Cas de l'analyseur
  \texttt{Multiple\-Isotropic\-Mises\-Flows}} Plusieurs blocs
\mkey{FlowRule} peuvent être définis dans le cas de l'analyseur
\texttt{Multiple\-Isotropic\-Mises\-Flows}. La directive \mkey{FlowRule}
est suivie d'un des trois types d'écoulement supportés, respectivement
\texttt{Creep}, \texttt{StrainHardeningCreep} et \texttt{Plasticity}. Le
bloc suivant décrit l'écoulement en suivant les règles données au
paragraphe précédent. Notons que la déclaration d'un nouvel écoulement
déclare automatiquement la déformation viscoplastique cumulée associée
sous le nom \texttt{pi} où \texttt{i} est le nombre d'écoulement défini
jusque là. Il n'est possible d'associer des noms de glossaire ou des
bornes à ces variables qu'après la définition du bloc.

\paragraph{Transformation du code dans les blocs \texttt{@FlowRule}} Le
code dans les blocs suivant la directive \mkey{FlowRule} est modifié
ainsi~:
\begin{itemize}
  \item les variables externes sont remplacées par leurs valeurs en
  \(t+\theta\,\Delta\,t\)~;
  \item la déformation inélastique cumulée \(p\) est remplacée par sa
  valeur en \(t+\theta\,\Delta\,t\).
\end{itemize}

\subsubsection{Mise à jour des variables auxiliaires}

La directive \mkey{UpdateAuxiliaryStateVariables} permet de mettre à
jour les variables auxiliaires.

\paragraph{Déformations inélastiques} Les déformations inélastiques
\(\tepsilonan_{i}\) ne sont pas des variables internes. Pour y avoir
accès (pour des posttraitements), il est possible de définir des
variables internes auxiliaires.

Ces variables sont mises à jour après les variables internes et les
contraintes, mais avant les variables externes (incluant la température)
ou la déformation totale\footnote{Les variables externes ne sont pas
  mises à jour car cela est du ressort de l'interface en cas de
  sous-découpage. Si aucun sous-découpage n'a lieu ou si l'on a réalisé
  le dernier pas de temps, on économise l'opération de mise à jour.}.

\paragraph{Exemple de la déformation inélastique totale} Supposons la
déformation inélastique totale représentée par une variable auxiliaire
tensorielle nommée \texttt{evp} (déclarée par la directive
\mkeyb{AuxiliaryStateVariable}{Auxiliary\-State\-Variable}). Cette
variable peut être calculée, dans le bloc suivant la directive
\mkeyb{UpdateAuxiliaryStateVariables}{Update\-Auxiliary\-State\-Variables}
ainsi~:
\begin{center}
  \texttt{evp += deto-deel;}
\end{center}
Une autre manière de calculer cette variable est~:
\begin{center}
  \texttt{evp = eto+deto-eel;}
\end{center}
qui montre que la déformation élastique a été mise à jour et non la
déformation totale.

\subsubsection{Paramètres numériques automatiquement définis}

Le paramètre \(\theta\) de la \(\theta\)-méthode vaut par défaut \(1\)
pour l'analyseur \texttt{Isotropic\-Plastic\-Mises\-Flow} et
\(\pfrac{1}{2}\) pour les autres. Cette valeur par défaut peut être
modifiée par la directive \mkey{Theta}. Cette valeur est également un
paramètre de la loi, nommé \varcpp{theta}, qui peut être modifié à
l'exécution.

La valeur du critère d'arrêt est par défaut de \(10^{-8}\). Cette
valeur par défaut peut être modifiée par la directive \mkey{Epsilon}.
Cette valeur est également un paramètre de la loi, nommé
\varcpp{epsilon}, qui peut être modifié à l'exécution.

Le nombre maximum d'itération de l'algorithme est de \(100\) par
défaut. Cette valeur par défaut peut être modifiée par la directive
\mkey{IterMax}. Cette valeur est également un paramètre de la loi, nommé
\varcpp{iterMax}, qui peut être modifié à l'exécution.

\subsubsection{Noms réservés}

Les noms réservés par cette analyseur sont décrits en
annexe~\ref{sec:noms-de-variables}.

\subsubsection{Comportement viscoplastique du \sic{}}

Nous voulons décrire le comportement viscoplastique du \sic{} Le
\sic{} est supposé avoir un comportement viscoplastique isotrope
donnée par~:
\[
\dot{p} = f\paren{\sigmaeq}
\]
où\begin{minipage}[t]{0.75\linewidth}
\begin{itemize}
  \item \(p\) est la déformation viscoplastique équivalente~;
  \item \(\sigmaeq\) est la contrainte de \nom{Von Mises}.
\end{itemize}
\end{minipage}

La fonction d'écoulement est donnée par~:
\begin{equation}
  \label{eq:sicBehaviour}
  f\paren{\sigmaeq}=\paren{A\,\exp\paren{-\Frac{B}{T}}+a\,\phi}\,\sigmaeq
\end{equation}
où~:
\begin{minipage}[t]{0.75\linewidth}
  \begin{itemize}
  \item \(A\), \(B\) et \(a\) sont des coefficients~;
  \item \(T\) est la température~;
  \item \(\phi\) est le flux de neutrons rapides~;
  \end{itemize}
\end{minipage}

\begin{figure}[htbp]
  \centering
  \code{{\ttfamily \input{@abs_top_srcdir@/docs/mfront/mfront/SICCreepBehaviour.tex}}}
  \caption{Implantation de la loi de comportement viscoplastique
    du \sic{} en \mfront{}.}
  \label{fig:mfront_sic_behaviour}
\end{figure}

Cette loi de comportement dépend d'une variable externe, le flux de
neutron rapide \(\phi\).

Pour implanter cette loi de comportement, nous utilisons l'analyseur
\texttt{IsotropicMisesCreep}. Le code source est donné en
figure~\ref{fig:mfront_sic_behaviour}.

\paragraph{Le mot clé \mkey{Parser}} La première ligne, commençant
par le mot clé \mkey{Parser}, décrit le type d'analyseur utilisé,
ici \texttt{IsotropicMisesCreep}.

\paragraph{Le mot clé \mkey{Behaviour}} La seconde ligne, commençant
par le mot clé \mkey{Behaviour}, donne le nom de la loi de
comportement.

\paragraph{Les mots clé \mkey{Author} et \mkey{Date}} La
troisième et la quatrième ligne renseignent respectivement l'auteur du
fichier et la date de création à l'aide des mots clés \mkey{Author}
et \mkey{Date}.

\paragraph{Le mot clé \mkey{Description}} Le mot clé
\mkey{Description} permet de donner les références bibliographiques
d'où la loi est extraite.

\paragraph{Le mot clé \texttt{@External\-State\-Variable}} Le mot clé
\mkey{External\-State\-Variable} définit une variable externe
scalaire (\texttt{real}) nommée \texttt{flux}. Cette variable, et son
incrément \texttt{dflux}, sont dès lors accessibles.

\paragraph{Le mot clé \mkey{StaticVariable}} Le mot clé
\mkey{Static\-Variable} sert à définir les constantes utilisées par
la loi. Ces constantes sont des scalaires (\texttt{real}).

\paragraph{Le mot clé \mkey{LocalVariable}} Le mot clé
\mkey{Local\-Variable} sert à définir des variables de travail
locales. Le rôle des variables \texttt{AF1} et \texttt{AF3} sera
explicité dans la suite.

\paragraph{Le mot clé \mkey{InitLocalVariables}} Le mot clé
\mkey{InitLocalVars} permet d'écrire du code appelé avant tout
calcul. Nous y initialisons les valeurs de coefficients de la loi. En
effet, la méthode d'intégration utilisée évalue la fonction \(f\) au
temps \(t+\theta\,\Delta t\) (intégration implicite).  Le coefficient
dépendant de la température est donc connu et il est avantageux de
l'évaluer ici plutôt qu'au cours des itérations de convergence de
l'algorithme (nous économisons des appels à la fonction
exponentielle). Nous utilisons ici le fait que la température au temps
\(t+\theta\,\Delta t\) est égale à \(T+\theta\,\Delta T\). De même la
valeur du flux de neutrons rapides \(\phi\) est égal à
\(\phi+\theta\,\Delta \phi\).

\paragraph{Le mot clé \mkey{FlowRule}} Le mot clé
\mkey{FlowRule} permet de renseigner la fonction d'écoulement
\(f\) et sa dérivée \(\derivtot{f}{\sigmaeq}\).

\clearpage
\newpage
\section{Intégration des lois de comportement, intégrateur par défaut}
\label{sec:defaultparser}

L'intégrateur par défaut est essentiellement utilisé à des fins de
test ou pour certaines lois spécifiques (loi d'endommagement pilotée
en déformations).

\subsection{Exemple de l'élasticité orthotrope}

\begin{figure}[htbp]
  \centering
  \begin{minipage}[htbp]{0.9\linewidth}
    \shorthandoff{:}
    \code{
      \small
      \input{@abs_top_srcdir@/docs/mfront/mfront/OrthotropicElastic.tex}
    }
    \shorthandon{:}
  \end{minipage}
  \caption{Implantation d'une loi orthotrope élastique}
  \label{fig:OrthoElastique}
\end{figure}

La figure~\ref{fig:OrthoElastique} décrit l'implantation d'une loi de
comportement orthotrope élastique. Cette loi ne nécessitant pas
d'algorithme d'intégration, l'analyseur par défaut est utilisé. La
commande \mkey{Orthotropic\-Behaviour} permet le support des lois
orthotropes. La commande \mkey{Require\-Stiff\-ness\-Tensor} demande à
ce que l'interface au code appelant mette à disposition une variable
{\tt D} contenant la matrice d'élasticité\footnote{Pour ce faire,
  l'interface déclare des propriétés matériaux supplémentaires. Leur
  nombre et l'ordre dans lequel il est nécessaire de les passer sont
  décrits dans les documentations des
  interfaces~\cite{helfer_interface_2013,helfer_interface_2013-1}.}.

Le calcul des contraintes en fin de pas de temps (à l'instant
\(t+\Delta\,t\)) s'écrit de manière similaire à son expression
mathématique~:
\[
\tsigma_{t+\Delta\,t} = \tenseurq{D}\colon\paren{\tepsilonto_{t}+\Delta\,\tepsilonto}
\]

Dans ce cas, le calcul de la matrice tangente cohérente est immédiat.

\subsection{Noms des variables et des méthodes \cpp{} utilisés en interne}

Pour éviter les conflits avec les noms de variables de travail ou des
méthodes \cpp{} utilisés par les classes générées, certains noms de
variables ne peuvent être utilisés. Leur liste est donnée en
annexe~\ref{sec:noms-communs-tous}.

\subsection{La directive \texttt{@Integrator}}

La directive \mkey{Integrator} permet d'intégrer la loi de comportement.

\paragraph{Conventions spécifiques} Les
conventions suivantes s'appliquent~:
\begin{itemize}
  \item {\tt Dt} représente la matrice tangente cohérente qu'il
  faut calculer~;
  \item {\tt sig} représente la contrainte qu'il faut calculer~;
  \item {\tt eto} représente la déformation totale en début de
  pas~;
  \item {\tt deto} représente l'incrément de déformation totale
  (constante sur le pas)~;
  \item {\tt T} représente la valeur de la température en début
  de pas~;
  \item {\tt dT} représente l'incrément de changement de
  température (constante sur le pas)~;
  \item pour toute variable interne {\tt Y}, {\tt Y} représente
  sa valeur en début de pas~;
  \item pour toute variable interne {\tt Y}, {\tt dY} représente
  l'incrément de cette variable sur le pas, incrément qu'il faut
  calculer~;
  \item pour toute variable externe {\tt V}, {\tt V} représente
  sa valeur en début de pas~;
  \item pour toute variable externe {\tt V}, {\tt dV} représente
    son incrément de variation sur le pas de temps (constante sur le pas).
\end{itemize}

\paragraph{Calcul d'une matrice de raideur} La loi de comportement
peut éventuellement fournir une matrice de raideur pour réaliser les
itérations de la résolution globale. Pour que les interfaces puissent
gérer cette possibilité, l'analyseur {\tt Default\-Parser} fournit
deux mots clés~:
\begin{itemize}
\item {\tt @Provides\-Tangent\-Operator} indique que la loi fournit
  une matrice tangente non symétrique~;
\item {\tt @Provides\-Symmetric\-Tangent\-Operator} indique que la loi
  fournit une matrice tangente symétrique~;
\end{itemize}

Si le code appelant demande le calcul de la matrice tangente
cohérente, la variable booléenne
\varcpp{compute\-Tangent\-Operator\textunderscore{}} est vraie.

Le type de matrice demandé est stocké dans la variable {\tt smt}
(\varcpp{stiffness matrix type}). L'utilisateur doit tester sa valeur
et effectuer le calcul le cas échéant. Sa valeur peut être~:
\begin{itemize}
\item {\tt ELASTIC}, pour la matrice d'élasticité (matrice
  d'élasticité)~;
\item {\tt SECANTOPERATOR}, pour la matrice sécante (matrice
  d'élasticité endommagée)~;
\item {\tt TANGENTOPERATOR},  pour la matrice tangente~;
\item {\tt CONSISTENTTANGENTOPERATOR}, pour la matrice tangente
  cohérente~;
\end{itemize}

\paragraph{Matrice de prédiction} Il n'est aujourd'hui pas possible de
préciser une matrice de prédiction avec cet analyseur.

\clearpage
\newpage
\section{Intégration des lois de comportement par
  une méthode explicite (algorithme de \nom{Runge-Kutta})}
\label{sec:RK}

Cette section décrit l'analyseur \texttt{RungeKutta} qui permet
l'intégration des lois de comportement, formulées en vitesse, par une
méthode explicite.

Contrairement aux analyseurs spécifiques décrits dans la
section~\ref{sec:mfront:isotropic:analyser}, aucune hypothèse n'est
faite sur les lois de comportement présentées~: l'analyseur
\texttt{RungeKutta} est dit générique.

Nous supposons que l'utilisateur a réussi à exprimer la loi de
comportement sous la forme d'un système différentiel~:
\begin{equation}
  \label{eq:systeme_diff}
  \dot{Y}=G\paren{Y,t}
\end{equation}
où \(G\) est une fonction {\em a priori} non linéaire et que nous
supposerons {\em a minima} continûment dérivable. Cette fonction se
construit en rassemblant les équations régissant les différentes
variables d'états.

\(Y\) un vecteur regroupant les différentes variables internes~:
\[
Y=
\begin{pmatrix}
  y_{1} \\
  \vdots \\
  y_{i} \\
  \vdots \\
  y_{n} \\
\end{pmatrix}
\]
Cette écriture est {\em symbolique} et chaque terme \(y_{j}\) peut
représenter une variable interne qui peut être soit un scalaire soit
un tenseur symétrique d'ordre \(2\). L'analyseur \texttt{RungeKutta}
impose que le premier terme de ce vecteur soit la déformation
élastique, dont la variable associée est \texttt{eel}, qui est
automatiquement déclarée. Le nom de glossaire de cette variable est
\texttt{Elastic\-Strain}.

La dépendance en temps de la fonction \(G\) qui apparaît dans
l'équation~\eqref{eq:systeme_diff} désigne en réalité une dépendance à
la variation de certaines variables {\em externes} qui influencent la
loi de comportement. Des exemples de telles variables externes sont~:
\begin{minipage}[t]{0.5\linewidth}
  \begin{itemize}
    \item le taux de combustion~;
    \item la taille de grain~;
    \item la densité de fission~;
    \item etc\ldots
  \end{itemize}
\end{minipage}

Nous détaillons les méthodes d'intégration proposées par l'analyseur
\texttt{RungeKutta}. Nous donnons en même temps les directives
\mfront{} qui permettent de programmer le système différentiel et de
modifier le comportement des algorithmes.

\subsection{Les algorithmes de \nom{Runge-Kutta}}
\label{sec:resol-dun-syst}

Les méthodes de \nom{Runge-Kutta} désignent une famille d'algorithmes
telles que~:
\[
\finpas{Y}=\debutpas{Y}+\displaystyle\sum_{i=1}^{n}b_{i}k_{i}
\]
où~:
\begin{minipage}[t]{0.9\linewidth}
  \(
  \left\{
  \begin{aligned}
    k_{1} &= \Delta\,t\,G\paren{\debutpas{Y},t}\\
    k_{2} &= \Delta\,t\,G\paren{\debutpas{Y}+a_{21}\,k_{1},t+c_{2}\,\Delta\,t}\\
    k_{3} &= \Delta\,t\,G\paren{\debutpas{Y}+a_{31}\,k_{1}+a_{32}\,k_{2},t+c_{3}\,\Delta\,t}\\
    &\vdots\\
    k_{n} &= \Delta\,t\,G\paren{\debutpas{Y}+a_{n1}\,k_{1}+\ldots+a_{n,n-1}\,k_{n-1},t+c_{n}\,\Delta\,t}\\
  \end{aligned}
  \right.
  \)
\end{minipage}

Une méthode particulière est caractérisée par la donnée du nombre
d'étapes \(n\) et des différents coefficients \(a_{ij}\), \(b_{i}\) et
\(c_{i}\).

L'un des intérêts des méthodes de \nom{Runge-Kutta} est que les
coefficients \(k_{i}\) d'une méthode d'ordre élevée peuvent parfois
être utilisés pour construire une méthode d'ordre moins élevée. On
obtient ainsi deux estimations de la valeur de fin de pas qui peuvent
être comparées~: si la comparaison entre les deux estimations n'est
pas satisfaisante, on en déduit que le pas de temps utilisé est trop
grand et il est possible de mettre en place une stratégie de
redécoupage du pas de temps.

\paragraph{Pas de temps minimal} En cas de convergence
forcée~\cite{pascal_notice_2005}, le code aux éléments finis \castem{}
envoie à la loi de comportement un pas de temps nul qui ne permet pas
d'écrire les lois de comportement sous la forme d'un système
différentiel. Il est nécessaire de préciser un pas de temps minimal
par la directive \mkeyb{MinimalTimeStep}{Minimal\-Time\-Step} pour
gérer ce cas. À l'exécution, la valeur choisie peut être modifiée à
l'aide du paramètre \varcpp{dtmin}.

\subsection{Algorithmes disponibles}


\begin{table}[htbp]
  \centering
  \begin{tabular}{|c|p{9cm}|}
    \hline
    Argument & Algorithme associé \\
    \hline
    \hline
    \texttt{Euler} & algorithme d'\nom{Euler} \\
    \hline
    \texttt{rk2} & algorithme de \nom{Runge-Kutta} d'ordre \(2\) \\
    \hline
    \texttt{rk4} & algorithme de \nom{Runge-Kutta} d'ordre \(4\) \\
    \hline
    \texttt{rk42} & algorithme avec contrôle de l'erreur et adaptation
    du pas de temps basé sur des méthodes de \nom{Runge-Kutta} d'ordre
    \(4\) et d'ordre \(2\) \\
    \hline
    \texttt{rk54} & algorithme avec contrôle de l'erreur et adaptation
    du pas de temps basé sur des méthodes de \nom{Runge-Kutta} d'ordre
    \(5\) et d'ordre \(4\) \\
    \hline
    \texttt{rkCastem} & algorithme avec contrôle de l'erreur et adaptation
    du pas de temps extrait du code \castem{} \\
    \hline
  \end{tabular}
  \caption{Arguments possibles de la directive \texttt{@Algorithm}.}
  \label{tab:Algorithm}
\end{table}

Différents algorithmes de \nom{Runge-Kutta} sont disponibles.
L'utilisateur peut choisir l'algorithme utilisé par la directive
\mkey{Algorithm}. Les arguments possibles de cette directive sont
regroupés au tableau~\ref{tab:Algorithm}.

\subsubsection{Algorithmes sans contrôle du pas de temps}

\mfront{} permet d'utiliser un certain nombre d'algorithmes de
\nom{Runge-Kutta} sans contrôle du pas de temps. L'utilisation de ces
algorithmes n'est pas recommandée car leurs prédictions sont
généralement peu satisfaisantes~: il est préférable d'utiliser les
algorithmes avec contrôle de l'erreur et pas de temps adaptatif.

\paragraph{Méthode d'\nom{Euler} explicite} L'estimation la plus simple
consiste à remplacer la dérivée \(\dot{Y}\) par la différence finie
\(\Frac{\finpas{Y}-\debutpas{Y}}{\Delta\, t}\), ce qui conduit à~:
\begin{equation}
  \label{eq:explicite}
  \finpas{Y}=\debutpas{Y}+\Delta\,t\,G\paren{\debutpas{Y},t}
\end{equation}
Cette estimation, dite {\em explicite}, est connue pour être à la fois
peu précise et peu robuste. L'erreur de cette méthode est
(asymptotiquement) proportionnelle au pas de temps~: elle est dite
d'ordre \(1\).

\paragraph{Méthode de \nom{Runge-Kutta} d'ordre 2}
Les méthodes de \nom{Runge-Kutta} permettent d'améliorer la précision
du schéma d'\nom{Euler}. Pour obtenir une méthode d'ordre \(2\), la
valeur \(\finpas{Y}\) est recherchée de la forme~:
\[
\finpas{Y} = \debutpas{Y}+a_{1}\Delta\,t\,G\paren{\debutpas{Y},t}+a_{2}\Delta\,t\,G\paren{\debutpas{Y}+a_{4}\Delta\, t,t+a_{3}\Delta\,t}
\]
Les quatre coefficients \(\left.a_{i}\right|_{i \in [0:4]}\) sont
choisis pour que l'expression précédente coïncide avec le
développement de \nom{Taylor} à l'ordre \(2\) de la fonction \(G\).
Cette contrainte conduisant à un système de \(3\) équations à \(4\)
inconnues, plusieurs choix sont possibles.

Le choix fait dans \mfront{} est dite \og~du point milieu~\fg{}. Ce
choix correspond à la solution \(a_{1}=0\), \(a_{2}=1\),
\(a_{3}=\pfrac{1}{2}\) et \(a_{4}=\pfrac{G\paren{Y_{i},t}}{2}\). Cette
méthode s'écrit ainsi~:
\begin{equation}
  \label{eq:RK2}
  \finpas{Y} = \debutpas{Y}+\Delta\,t\,G\paren{\debutpas{Y}+\Frac{\Delta\,
      t}{2}G\paren{\debutpas{Y},t},t+\Frac{\Delta\, t}{2}}
\end{equation}

\paragraph{Méthode de \nom{Runge-Kutta} d'ordre $4$}

La méthode de \nom{Runge-Kutta} d'ordre \(4\) retenue dans \mfront{} est
donnée par~:
\begin{equation}
  \label{eq:RK4}
  \finpas{Y} = \debutpas{Y}+\Frac{1}{6}\paren{k_{1}+2\,k_{2}+2\,k_{3}+k_{4}}
\end{equation}
où~:
\begin{minipage}[t]{0.9\linewidth}
  \(
  \left\{
  \begin{aligned}
    k_{1} &= \Delta\,t\,G\paren{\debutpas{Y},t}\\
    k_{2} &= \Delta\,t\,G\paren{\debutpas{Y}+\Frac{1}{2}\,k_{1},t+\Frac{1}{2}\,\Delta\,t}\\
    k_{3} &= \Delta\,t\,G\paren{\debutpas{Y}+\Frac{1}{2}\,k_{2},t+\Frac{1}{2}\,\Delta\,t}\\
    k_{4} &= \Delta\,t\,G\paren{\debutpas{Y}+k_{3},t+\Delta\,t}
  \end{aligned}
  \right.
  \)
\end{minipage}

\subsubsection{Algorithmes avec contrôle du pas de temps (correcteur/prédicteur)}

\mfront{} propose trois algorithmes de \nom{Runge-Kutta} avec
sous-découpage du pas de temps pour garantir la qualité de la solution
trouvée. Pour cela, on introduit un temps local \(t^{l}\), compris
entre\(t\) et \(t+\Delta\,t\), et un pas de temps local
\(\delta\,t^{l}\) variable permettant d'augmenter \(t^{l}\).

\paragraph{Méthode de \nom{Runge-Kutta} d'ordre $\pfrac{4}{2}$} Les
algorithmes~\eqref{eq:RK2} et~\eqref{eq:RK4} partagent les mêmes
coefficients \(k_{1}\) et \(k_{2}\)~: si l'on utilise une méthode
d'ordre \(4\), on a \og~gratuitement\fg{} une estimation d'ordre \(2\).
La différence entre ces deux estimations est en \(\bigO 3\) et peut
servir à contrôler le pas de temps.

\paragraph{Méthode de \nom{Runge-Kutta} d'ordre $\pfrac{5}{4}$ (méthode
  de \nom{Fehlberg})} Pour décrire cette méthode, utilisée par défaut,
nous reprenons la présentation d'\nom{A. Fortin}~\cite{fortin_analyse_2001}.
Il est utile d'introduire les \(6\) constantes suivantes~:
\[
\begin{aligned}
  k_{1} &= \Delta\,t\,\cdot\,G\paren{\debutpas{Y_{i}},t} \\
  k_{2} &=
  \Delta\,t\,\cdot\,G\paren{\debutpas{Y_{i}}+\Frac{1}{4}\,k_{1},t+\Frac{1}{4}\,\Delta\,t}\\
 k_{3} &=
  \Delta\,t\,\cdot\,G\paren{\debutpas{Y_{i}}+\Frac{3}{32}\,k_{1}+\Frac{9}{32}\,k_{2},t+\Frac{3}{8}\,\Delta\,t}
 \\
  k_{4} &=
  \Delta\,t\,\cdot\,G\paren{\debutpas{Y_{i}}+\Frac{1932}{2197}\,k_{1}-\Frac{7200}{2197}\,k_{2}+\Frac{7296}{2197}\,k_{3},t+\Frac{12}{13}\,\Delta\,t}\\
 k_{5} &=
  \Delta\,t\,\cdot\,G\paren{\debutpas{Y_{i}}+\Frac{439}{216}\,k_{1}-8\,k_{2}+\Frac{3680}{513}\,k_{3}-\Frac{845}{4104}\,k_{4},t+\Delta\,t}\\
 k_{6} &=
  \Delta\,t\,\cdot\,G\paren{\debutpas{Y_{i}}-\Frac{8}{27}\,k_{1}+2\,k_{2}-\Frac{3544}{2565}\,k_{3}+\Frac{1859}{4104}\,k_{4}-\Frac{11}{40}\,k_{5},t+\Frac{1}{2}\Delta\,t}\\
\end{aligned}
\]

Ces coefficients permettent de construire une méthode de
\nom{Runge-Kutta} d'ordre \(4\)~:
\[
\finpas{Y_{i}}=\debutpas{Y_{i}}+\paren{\Frac{25}{216}\,k_{1}+\Frac{1408}{2565}\,k_{3}+\Frac{2197}{4104}\,k_{4}-\Frac{1}{5}\,k_{5}}
\]
et une autre d'ordre \(5\)~:
\[
\finpas{Y_{i}}=\debutpas{Y_{i}}+\paren{\Frac{16}{135}\,k_{1}+\Frac{6656}{12\,825}\,k_{3}+\Frac{28\,561}{56\,430}\,k_{4}-\Frac{9}{50}\,k_{5}+\Frac{2}{55}\,k_{6}}
\]

La différence entre ces deux estimations est en \(\bigO 5\) et peut
servir à contrôler le pas de temps.

\paragraph{Méthode \castem{}} Une méthode de \nom{Runge-Kutta}
particulière est disponible dans \castem{}. Cette méthode se distingue
par sa gestion du pas de temps et son critère d'arrêt qui sont
uniquement basés sur la différence des contraintes entre deux ordres
d'approximations (et non sur les variables internes). Cette méthode est
décrite en annexe~\ref{sec:algoRKcastem}.

\subsubsection{Gestion du pas de temps}

Les méthodes de \nom{Runge-Kutta} d'ordre $\pfrac{4}{2}$ et d'ordre
$\pfrac{5}{4}$ permettent de comparer les estimations fournies par
deux méthodes d'ordre différents. La différence entre ces deux
estimations mesure l'erreur d'intégration et est comparée à une
certaine tolérance \(\varepsilon\)~:
\begin{equation}
  \label{eq:RK:epsilon}
  \norm{Y^{(a)}-Y^{(b)}}_{\text{RK}}<\varepsilon\quad\text{avec}\quad
\end{equation}

Deux mesures de l'erreur sont utilisées. Pour un nombre restreint de
variables internes\footnote{Cette mesure est utilisée si la taille du
  vecteur \(Y\) est inférieure à \(20\) en \(1D\).}, la mesure
utilisée est la suivante~:
\begin{equation}
  \norm{Y^{(a)}-Y^{(b)}}_{\text{RK}}<\varepsilon\quad\text{avec}\quad
\norm{Y^{(a)}-Y^{(b)}}_{\text{RK}}=
  \Frac{1}{\displaystyle\sum_{i=1}^{N}N_{y_{i}}}\displaystyle\sum_{i=1}^{N}\Frac{\left|
    y_{i}^{(a)}-y_{i}^{(b)}\right|}{E_{y_{i}}}
\end{equation}

Pour des systèmes différentiels de plus grandes tailles, la norme
suivante est utilisée~:
\begin{equation}
  \norm{Y^{(a)}-Y^{(b)}}_{\text{RK}}=
  \Frac{1}{\displaystyle\sum_{i=1}^{N}N_{y_{i}}}\displaystyle\sum_{i=1}^{N}\Frac{\left|
    y_{i}^{(a)}-y_{i}^{(b)}\right|}{E_{y_{i}}}
\end{equation}

Dans ces deux expressions, nous avons noté~:\\
\begin{minipage}[t]{0.9\linewidth}
  \begin{itemize}
    \item \(N\) est le nombre de variables internes~;
    \item \(N_{y_{i}}\) le nombre de composantes de la variable interne
    \(y_{i}\)~;
    \item \(y_{i}^{(a)}\) est l'estimation de la variable interne en fin
    de pas de temps obtenue par la première méthode~;
    \item \(y_{i}^{(b)}\) est l'estimation de la variable interne en fin
    de pas de temps obtenue par la seconde méthode~;
    \item \(E_{y_{i}}\) est un paramètre de normalisation associé à la
    variable \(E_{y_{i}}\), si l'utilisateur en a définit une par la
    méthode \texttt{set\-Error\-Norma\-lisa\-tion\-Factor}. Ce paramètre
    permet de normaliser toutes les variables à quelque chose de l'ordre
    des déformations. Par exemple, si la variable \(y_{i}\) est de
    l'ordre des contraintes, il est possible de la normaliser ainsi~:
    \begin{center}
      \texttt{yi.setErrorNormalisationFactor(young);}
    \end{center}
    L'argument de la méthode
    \texttt{set\-Error\-Norma\-lisa\-tion\-Factor} peut être soit une
    valeur scalaire soit une propriété matériau.
  \end{itemize}
\end{minipage}

Si l'inégalité~\eqref{eq:RK:epsilon} n'est pas respectée, le pas de
temps \(\delta\,t^{l}\) est sous-découpé automatiquement. Inversement,
si elle est respectée, le pas de temps \(\delta\,t^{l}\) est augmenté.

Supposons que l'erreur soit en \(\bigO n\), un facteur \(\beta\)
d'augmentation ou de réduction du pas de temps peut être proposé ainsi~:
\[
\beta = \paren{\Frac{\varepsilon}{\norm{Y^{(a)}-Y^{(b)}}_{\text{RK}}}}^{\Frac{1}{n}}
\]

L'augmentation ou la réduction du pas de temps local est limité à un
facteur \(10\). Naturellement, le pas de temps local \(\delta\,t^{l}\)
ne peut être supérieur à \(t+\Delta\, t-t^{l}\).

\paragraph{Gestion par l'utilisateur} Si le code
suivant la directive \mkey{Derivative}, présenté dans la suite et qui
sert à programmer le système différentiel, renvoie la valeur booléenne
\texttt{false}, le pas de temps est automatiquement divisé par \(2\).

Cette possibilité est en pratique assez importante, en particulier
pour des comportements très raides. Il est par exemple fréquent de
rencontrer des lois viscoplastique de type \nom{Norton} dont
l'exposant est grand. Dans ces cas, les premières estimations des
solutions conduisent à des contraintes équivalentes assez importantes
qui conduisent à des vitesses d'écoulement infinies\footnote{Ceci est
  particulièrement vrai quand on utilise le système international où
  les contraintes s'expriment en Pascal. Une autre manière de résoudre
  le problème est de réécrire la loi de \nom{Norton} différemment. Au
  lieu d'écrire \(\dot{p}=A\,\sigmaeq^{n}\), il faudrait mieux écrire
  \(\dot{p}=\dot{p}_{0}\,\paren{\Frac{\sigmaeq}{\sigma_{0}}}^{n}\)
  avec \(\dot{p}_{0}=A\,\sigmaeq^{0}\).}.  Informatiquement, ces
vitesses se traduisent par l'apparition de la valeur flottante
\varcpp{NaN}\footnote{Acronyme de \og~Not a Number~\fg{}}. Les
algorithmes de \nom{Runge-Kutta} captent ce cas et sous-découpent
alors automatiquement le pas de temps, mais cette solution présente
l'inconvénient de ne pas être compatible avec certaines techniques de
déverminage extrêmement utiles. Il est dès lors préférable, et
également plus efficace puisque des opérations lourdes (calculs de
puissance) peuvent être évitées, que l'utilisateur introduise des
tests appropriés.

\paragraph{Choix de la tolérance} La tolérance
\(\varepsilon\) est fixée par défaut à \(10^{-8}\). L'utilisateur peut
modifier cette valeur par la directive \mkey{Epsilon}. À l'exécution,
l'utilisateur peut modifier cette valeur à l'aide du paramètre
\varcpp{epsilon}.

\subsubsection{Critère d'arrêt} L'intégration
s'arrête quand l'une des conditions suivantes est vérifiée~:
\begin{itemize}
  \item le temps local devient trop petit, ce qui marque l'échec
  de l'algorithme. Le pas de temps minimal admissible est égal au pas de
  temps total \(\Delta\,t\) multiplié par 100 fois la précision machine
  du type réel considéré~;
  \item le temps local est égal au temps de fin de pas, ce qui
  marque la réussite de l'intégration. Ceci est vérifié par~:
  \[
  \left|t+\Delta\,
  t-t^{l}\right|<\Frac{\delta\,t^{l}}{2}
  \]
\end{itemize}

\subsection{La directive \texttt{@ComputeStress}}

La directive \mkey{ComputeStress} permet de calculer les contraintes aux
différentes étapes de la méthode et en fin de calcul. Les contraintes
sont supposées fonction des valeurs actuelles des variables internes.

\paragraph{Conventions spécifiques aux différentes étapes de la méthode}
Dans la méthode \mkey{ComputeStress}, les conventions suivantes
s'appliquent~:
\begin{itemize}
  \item {\tt sig} représente la contrainte qu'il faut calculer~;
  \item {\tt eto} représente la déformation totale actualisée~;
  \item {\tt deto} représente la vitesse de déformation totale
  (constante sur le pas)~;
  \item {\tt T} représente la valeur actualisée de température~;
  \item {\tt dT} représente la vitesse de changement de température
  (constante sur le pas)~;
  \item pour toute variable interne {\tt Y}, {\tt Y} représente sa
  valeur actualisée~;
  \item pour toute variable interne auxiliaire {\tt Y}, {\tt Y}
  représente sa valeur en {\em début} de pas de temps~;
  \item pour toute variable externe {\tt V}, {\tt V} représente sa
  valeur actualisée~;
  \item pour toute variable externe {\tt V}, {\tt dV} représente sa
  vitesse de variation sur le pas de temps (constante sur le pas).
\end{itemize}

\paragraph{Conventions spécifiques en fin d'intégration} Dans la méthode
\mkey{ComputeStress}, les conventions suivantes s'appliquent~:
\begin{itemize}
\item {\tt sig} représente la contrainte qu'il faut calculer~;
\item {\tt eto} représente la déformation totale en fin de pas~;
\item {\tt deto} représente la vitesse de déformation totale
  (constante sur le pas)~;
\item {\tt T} représente la valeur de la température en fin de pas~;
\item {\tt dT} représente la vitesse de changement de température
  (constante sur le pas)~;
\item pour toute variable interne {\tt Y}, {\tt Y} représente sa
  valeur en fin de pas~;
\item pour toute variable interne auxiliaire {\tt Y}, {\tt Y}
  représente sa valeur en {\em début} de pas de temps~;
\item pour toute variable externe {\tt V}, {\tt V} représente sa
  valeur en fin de pas~;
\item pour toute variable externe {\tt V}, {\tt dV} représente sa
  vitesse de variation sur le pas de temps (constante sur le pas).
\end{itemize}

\subsection{La directive \texttt{@Derivative}}

La directive \mkey{Derivative} permet de préciser les équations
différentielles associées aux différentes variables internes.

Le code fourni est appelé après celui fourni par la méthode
\mkey{ComputeStress}.

\paragraph{Conventions spécifiques} Dans la méthode \mkey{Derivative},
les conventions suivantes s'appliquent~:
\begin{itemize}
  \item {\tt sig} représente la valeur actualisée de la contrainte~;
  \item {\tt eto} représente la déformation totale actualisée~;
  \item {\tt deto} représente la vitesse de déformation totale
  (constante sur le pas)~;
  \item {\tt T} représente la valeur actualisée de température~;
  \item {\tt dT} représente la vitesse de changement de température
  (constante sur le pas)~;
  \item pour toute variable interne {\tt Y}, {\tt Y} représente sa
  valeur actualisée~;
  \item pour toute variable interne auxiliaire {\tt Y}, {\tt Y}
  représente sa valeur en {\em début} de pas de temps~;
  \item pour toute variable externe {\tt V}, {\tt V} représente sa
  valeur actualisée~;
  \item pour toute variable externe {\tt V}, {\tt dV} représente sa
  vitesse de variation sur le pas de temps (constante sur le pas).
\end{itemize}

Pour toute variable interne {\tt Y}, {\tt dY} représente sa vitesse de
variation~: c'est cette variable que l'utilisateur doit renseigner.

%\subsection{La directive \texttt{@TangentOperator}}
%
%La directive \mkey{TangentOperator} permet de préciser une matrice
%tangente.

%Le code fourni est appelé après celui fourni par la méthode
%\mkey{ComputeStress}.
%
%\paragraph{Conventions spécifiques} Dans la méthode \mkey{Derivative},
%les conventions suivantes s'appliquent~:
%\begin{itemize}
%  \item {\tt sig} représente la valeur actualisée de la contrainte~;
%  \item {\tt eto} représente la déformation totale actualisée~;
%  \item {\tt deto} représente la vitesse de déformation totale
%  (constante sur le pas)~;
%  \item {\tt T} représente la valeur actualisée de température~;
%  \item {\tt dT} représente la vitesse de changement de température
%  (constante sur le pas)~;
%  \item pour toute variable interne {\tt Y}, {\tt Y} représente sa
%  valeur actualisée~;
%  \item pour toute variable interne auxiliaire {\tt Y}, {\tt Y}
%  représente sa valeur en {\em début} de pas de temps~;
%  \item pour toute variable externe {\tt V}, {\tt V} représente sa
%  valeur actualisée~;
%  \item pour toute variable externe {\tt V}, {\tt dV} représente sa
%  vitesse de variation sur le pas de temps (constante sur le pas).
%\end{itemize}
%
%Pour toute variable interne {\tt Y}, {\tt dY} représente sa vitesse de
%variation~: c'est cette variable que l'utilisateur doit renseigner.

\subsection{Mise à jour des variables auxiliaires}

Les variables auxiliaires sont mises à jour à chaque fois qu'un pas est
jugé valide au sens du critère d'arrêt~\eqref{eq:RK:epsilon} après la
mise à jour des variables internes.

Cette méthode peut donc être appelée plusieurs fois si un algorithme à
pas de temps adaptatif est utilisé. Une variable spécifique, nommée
\texttt{dt\textunderscore{}}, permet de récupérer l'incrément de temps
effectivement utilisé par l'algorithme. Le pas de temps total est
donné par la variable \texttt{dt}.

\subsection{Paramètres automatiquement définis}

L'analyseur implicite définit automatiquement différents
paramètres~:
\begin{itemize}
  \item \varcpp{epsilon}, la valeur du critère d'arrêt de
  l'algorithme de \nom{Runge-Kutta}. La valeur par défaut de ce
  paramètre est de \(10^{-8}\) et peut être modifiée par la directive
  \mkey{Epsilon}. Ce paramètre doit avoir une valeur strictement
  positive~;
  \item \varcpp{dtmin}, la valeur minimale du pas de temps
  utilisé pour l'intégration\footnote{Ce paramètre a été introduit
    essentiellement pour contourner les difficultés liées à l'algorithme
    de convergence forcée du code \castem{} qui envoie, lorsqu'il
    s'active, un incrément de temps nul à la loi de comportement.}. Par
  défaut, ce paramètre n'existe que si la directive
  \mkeyb{MinimalTimeStep}{Minimal\-Time\-Step} est utilisée. Ce
  paramètre doit avoir une valeur strictement positive.
\end{itemize}

\subsection{Notes sur l'utilisation des tableaux de variables internes}

L'utilisation des tableaux de variables internes dans le
cas des algorithmes de \nom{Runge-Kutta} mérite quelques précisions.

Pour réduire le temps de développement de cette fonctionnalité et
surtout réduire les sources d'erreurs, nous avons utilisé des classes
livrés avec TFEL (la librairie mathématique sur laquelle se base
\mfront{}). En fonction de la taille des tableaux, deux types de
vecteurs sont utilisés~:
\begin{itemize}
\item pour des tableaux de taille inférieure à \(10\), des vecteurs de
  taille finie, alloués sur la pile, sont utilisés~;
\item pour des tableaux plus grands, des vecteurs alloués
  dynamiquement ont été utilisés.
\end{itemize}
Dans les deux cas, les opérations mathématiques usuelles ont été
définies~: l'addition d'un vecteur de tenseur est possible et
optimisée.

Ainsi, la génération des différents algorithmes de \textsc{Runge-Kutta}
est restée inchangée, ce qui est satisfaisant.

\subsection{Intégration d'une loi d'écoulement viscoplastique
  orthotrope par un algorithme de \nom{Runge-Kutta}}
\label{sec:integration-dune-loi}

Nous nous intéressons maintenant à l'écriture d'une loi de
comportement viscoplastique orthotrope. La forme retenue est proche
d'une loi de \nom{Norton}, la vitesse de l'écoulement étant donnée
par~:
\[
\tdepsilonvis = A\,\paren{\sigmaH}^{E}\,\tenseur{n}
\]
où \(sigmaH\) désigne la contrainte de \nom{Hill} et \(\tenseur{n}\)
le tenseur normal. Ils sont définis en annexe au
paragraphe~\ref{sec:critere-de-nomhill}.

Il est d'usage d'introduire une déformation viscoplastique cumulée
\(p\) dont la vitesse est donnée par~:
\[
\dot{p} = A\paren{\sigmaH}^{E}
\]

Le système d'équations complet est donc~:
\begin{equation}
  \label{eq:law}
  \begin{aligned}
    \dot{p}       &= A\paren{\sigmaH}^{E} \\
    \tdepsilonvis &= \dot{p}\,\tenseur{n}   \\
    \tdepsilonel  &= \tdepsilonto - \tdepsilonvis
  \end{aligned}
\end{equation}

\begin{figure}[htbp]
  \centering
  \begin{minipage}[htbp]{0.9\linewidth}
    \shorthandoff{:}
    \code{
      \small
      \input{@abs_top_srcdir@/docs/mfront/mfront/OrthotropicCreep.tex}
    }
    \shorthandon{:}
  \end{minipage}
  \caption{Implantation d'une loi d'écoulement viscoplastique
    orthotrope par une mé\-tho\-de de \nom{Runge-Kutta}}
  \label{fig:OrthoViscoRK}
\end{figure}

\paragraph{Implantation} La figure~\ref{fig:OrthoViscoRK} présente
l'implantation de ce système. Le calcul du tenseur de \nom{Hill} a été
présenté au paragraphe~\ref{sec:calcul-du-tenseur}. Rappelons
simplement qu'il est nécessaire d'inclure le fichier d'entête {\tt
  TFEL/\-Material\-Law/\-Hill.hxx} dans la directive {\tt
  @Includes}. Le calcul du tenseur normal est immédiat d'après les
formules données au paragraphe~\ref{sec:critere-de-nomhill}.

Deux remarques peuvent être faites~:
\begin{itemize}
  \item les intégrations de type \nom{Runge-Kutta} conduisent souvent à
  des prédictions initiales élevées des contraintes. Nous testons ici si
  la contrainte de \nom{Hill} est supérieure à \(10^{9}\, MPa\)~: dans
  ce cas, nous renvoyons la valeur {\tt false} qui entraîne un découpage
  du pas de temps~;
  \item cette implantation tient compte des limitations du code aux
  éléments finis \castem{} qui ne permet pas d'avoir une définition
  homogène des directions d'orthotropie quelque soit l'hypothèse de
  modélisation. Les conventions utilisées ici sont conformes à celles
  décrites en annexe~\ref{sec:annexe:orthotropie} pour les tubes.
\end{itemize}

\clearpage
\newpage
\section{Intégration des lois de comportement par une méthode implicite}
\label{sec:Implicite}

Nous décrivons dans cette section l'intégration des lois de
comportement par une méthode implicite. Ces méthodes sont très
efficaces, stables et peuvent être utilisées pour intégrer toutes les
lois de comportement.

Elles présentent de nombreux avantages par rapport aux autres
méthodes. Les plus importants sont que~:
\begin{itemize}
  \item ces méthodes introduisent naturellement la matrice
  tangente cohérente. De plus, dans la plupart des cas, cette matrice
  est pratiquement \og{}~gratuite~\fg{} (voir
  paragraphe~\ref{sec:mfront:implicit:tangentOperator})~;
  \item elles permettent de traiter naturellement les lois de
  comportement indépendantes du temps. Un exemple est donné par le
  traitement de la plasticité dans les analyseurs spécifiques (voir
  paragraphe~\ref{sec:mfront:isotropic:analyser}).
\end{itemize}

Ces méthodes supposent que l'on soit capable de résumer l'intégration
de la loi de comportement à la résolution d'un système non-linéaire.
Nous avons détaillé comment procéder pour les analyseurs spécifiques
(voir paragraphe~\ref{sec:mfront:isotropic:analyser}). Nous montrons
dans le paragraphe suivant comment procéder pour transformer un
système différentiel en système d'équations non linéaires et nous en
profitons pour introduire les méthodes implicites en se basant sur les
méthodes de \nom{Runge-Kutta}.

Deux analyseurs, nommés \texttt{Implicit} et \texttt{ImplicitII}, sont
disponibles pour intégrer les lois de comportement par une méthode
implicite. Ils se distinguent par le fait que l'analyseur
\texttt{Implicit} déclare automatiquement la déformation élastique
comme première variable interne\footnote{Nous verrons plus loin une
  justification de ce choix~: il permet un calcul quasi-automatique de
  la matrice tangente cohérente.}.

Différentes méthodes pour résoudre le système non-linéaire
introduit par la méthode implicite sont disponibles~:
\begin{itemize}
  \item \varcpp{NewtonRaphson}, l'algorithme classique de
  \nom{Newton-Raphson}~;
  \item \varcpp{NewtonRaphson\textunderscore{}NumericalJocabian},
  l'algorithme classique de \nom{Newton-Raphson} avec une approximation
  numérique du jacobien~;
  \item \varcpp{Broyden}, premier algorithme de \nom{Broyden}~;
  \item \varcpp{BroydenII}, second algorithme de \nom{Broyden}.
\end{itemize}
L'utilisateur choisi l'algorithme par la directive
\mkey{Algorithm}.

Outre la présentation des bases mathématiques de ces méthodes, nous
profitons de ce paragraphe pour justifier certaines initialisations
implicites faites dans \mfront{} qui rendent l'écriture des lois de
comportement plus concise.

\subsection{Résolution d'un système différentiel par une méthode implicite}
\label{sec:resol-dun-syst-implicite}

\subsubsection{Généralités} Soit à résoudre le système différentiel
suivant~:
\begin{equation}
  \label{eq:systeme_diff-implicite}
  \dot{Y}=G\paren{Y,t}
\end{equation}
où \(G\) est une fonction {\em a priori} non linéaire et que nous
supposerons {\em a minima} continûment dérivable.

La dépendance en temps évoquée ici dans la fonction \(G\) désigne en
réalité une dépendance à la variation de certaines variables {\em
  externes} qui influencent la loi de comportement. Des exemple de
telles variables externes sont~:
\begin{minipage}[t]{0.5\linewidth}
  \begin{itemize}
    \item le taux de combustion~;
    \item la taille de grain~;
    \item la densité de fission~;
    \item etc\ldots
  \end{itemize}
\end{minipage}

\(Y\) un vecteur regroupant les différentes variables internes~:
\[
Y=
\begin{pmatrix}
  y_{1} \\
  \vdots \\
  y_{i} \\
  \vdots \\
  y_{n} \\
\end{pmatrix}
\]
Cette écriture est {\em symbolique} et chaque terme \(y_{j}\) peut
représenter une variable interne qui peut être soit un scalaire soit un
tenseur symétrique d'ordre \(2\).

L'analyseur \texttt{Implicit} impose que le premier terme de ce
vecteur soit la déformation élastique, dont la variable associée est
\texttt{eel}, qui est automatiquement déclarée\footnote{Le nom de
  glossaire de cette variable est \texttt{Elastic\-Strain}. La notion
  de nom de glossaire est décrite dans la première partie de la
  documentation de \mfront{}~\cite{helfer_generateur_2013-1}.}.

La résolution numérique de cette équation consiste à estimer, à partir
de la valeur \(\debutpas{Y}\) à un instant \(t\), la valeur
\(\finpas{Y}\) à un instant \(t+\Delta\, t\).

\paragraph{Méthode d'\nom{Euler}}
L'estimation la plus simple consiste à remplacer la dérivée
\(\dot{Y}\) par la différence finie
\(\Frac{\finpas{Y}-\debutpas{Y}}{\Delta\, t}\), ce qui conduit à~:
\begin{equation}
  \label{eq:explicite-implicite}
  \finpas{Y}=\debutpas{Y}+\Delta\,t\,G\paren{\debutpas{Y},t}
\end{equation}
Cette estimation, dite {\em explicite}, est connue pour être à la fois
peu précise et peu robuste. L'erreur de cette méthode est
(asymptotiquement) proportionnelle au pas de temps~: elle est dite
d'ordre \(1\).

\paragraph{Méthode de \nom{Runge-Kutta} (explicite) d'ordre 2}
Les méthodes de \nom{Runge-Kutta} permettent d'améliorer la précision
du schéma d'\nom{Euler}. Pour obtenir une méthode d'ordre \(2\), la
valeur \(\finpas{Y}\) est recherchée de la forme~:
\[
\finpas{Y} = \debutpas{Y}+a_{1}\Delta\,t\,G\paren{\debutpas{Y},t}+a_{2}\Delta\,t\,G\paren{\debutpas{Y}+a_{4}\Delta\, t,t+a_{3}\Delta\,t}
\]
Les quatre coefficients \(\left.a_{i}\right|_{i \in [0:4]}\) sont
choisis pour que l'expression précédente coïncide avec le
développement de \nom{Taylor} à l'ordre \(2\) de la fonction \(G\).
Cette contrainte conduisant à un système de \(3\) équations à \(4\)
inconnues, plusieurs choix sont possibles.

\paragraph{Méthode du point milieu} Ce choix correspond à la solution
\(a_{1}=0\), \(a_{2}=1\), \(a_{3}=\pfrac{1}{2}\) et
\(a_{4}=\pfrac{G\paren{Y_{i},t}}{2}\). Cette méthode s'écrit ainsi~:
\begin{equation}
  \label{eq:mfront:implicite:midpoint}
  \finpas{Y} = \debutpas{Y}+\Delta\,t\,G\paren{\debutpas{Y}+\Frac{\Delta\, t}{2}G\paren{\debutpas{Y},t},t+\Frac{\Delta\, t}{2}}
\end{equation}

\paragraph{Méthode d'\nom{Euler} modifiée} Ce choix correspond à la solution
\(a_{1}=a_{2}=\pfrac{1}{2}\), \(a_{3}=1\) et
\(a_{4}=G\paren{Y_{i},t}\). Cette méthode s'écrit ainsi~:
\begin{equation}
  \label{eq:mfront:implicite:euler}
  \finpas{Y} = \debutpas{Y}+\Frac{\Delta\,t}{2}\left[G\paren{\debutpas{Y},t}+G\paren{\debutpas{Y}+\Delta\, t\,G\paren{\debutpas{Y},t},t+\Delta\,t }\right]
\end{equation}

\paragraph{Estimation implicite}
Des estimations dites {\em implicites} de la solution consistent à
introduire, à la place de l'estimation \(\debutpas{Y}+\Delta\, t\,
G\paren{\debutpas{Y},t}\), la valeur recherchée \(\finpas{Y}\).  Deux
schémas implicites, se déduisent alors des schémas de
\nom{Runge-Kutta} précédent.

À partir de la méthode du point
milieu~\eqref{eq:mfront:implicite:midpoint}, nous obtenons~:
\[
\finpas{Y}
=\debutpas{Y}+\Delta\,t\,G\paren{\Frac{1}{2}\left(\debutpas{Y}+\finpas{Y}\right),t+\Frac{\Delta\,
 t}{2}}
\]

À partir de la méthode d'\nom{Euler}
modifiée~\eqref{eq:mfront:implicite:euler}, nous obtenons~:
\[
\finpas{Y} =\debutpas{Y}+\Frac{\Delta\,t}{2}\,\left[G\paren{\debutpas{Y},t}
+\,G\paren{\finpas{Y},t+\Delta\,t}\right]
\]

\paragraph{\(\theta\)-méthodes}
Les schémas implicites précédents sont généralement présentés en
introduisant un paramètre \(\theta\), compris entre \(0\) et \(1\),
dans les expressions précédentes afin de rendre les apports des termes
\(\debutpas{Y}\) et \(\finpas{Y}\) dissymétriques.

À partir de la méthode du point
milieu~\eqref{eq:mfront:implicite:midpoint}, nous obtenons~:
\begin{equation}
  \label{eq:implicite}
  \finpas{Y} =\debutpas{Y}+\Delta\,t\,G\paren{\paren{1-\theta}\debutpas{Y}+\theta\,\finpas{Y},t+\theta\,\Delta\, t} \\
\end{equation}

La seconde forme, obtenue à partir de la méthode d'\nom{Euler}
modifiée~\eqref{eq:mfront:implicite:euler}, nous obtenons~:
\[
\finpas{Y} =\debutpas{Y}+\Delta\,t\,\left[\paren{1-\theta}G\paren{\debutpas{Y},t}
+\theta\,G\paren{\finpas{Y},t+\Delta\,t}\right]
\]
Cette seconde forme d'équation implicite est rarement utilisée, au moins
pour l'intégration des lois de comportement mécaniques.

\subsection{Méthode implicite utilisée dans \mfront{}}

\mfront{} tend, suivant l'usage, à supporter naturellement la
formulation~\eqref{eq:implicite}, mais rien n'interdit d'utiliser la
première. Pour simplifier la présentation, nous ne décrirons le
principe de la résolution que pour la première forme.

L'équation~\eqref{eq:implicite} peut s'écrire~:
\begin{equation}
  \label{eq:F_implicite} F\paren{\Delta\, Y} = \Delta Y -
  \Delta\,t\,G\paren{\debutpas{Y}+\theta\Delta{Y},t+\theta\,\Delta t} =
  0
\end{equation}
où nous avons introduit l'incrément \(\Delta\,Y\) de la variable \(Y\)
au cours du pas de temps \(\Delta\,t\), variable qui apparaît comme
l'inconnue \og~naturelle~\fg{} de cette équation. Nous avons donc ramené
la résolution du système différentiel~\eqref{eq:systeme_diff-implicite}
à la recherche d'un zéro de la fonction \(F\paren{\Delta\, Y}\).

\paragraph{Convention d'écriture} Il est d'usage d'introduire des
notations (et des abus de langage) pratiques.

Si dans le système différentiel initial, une fonction \(f\) des
variables internes \(f\paren{Y}\) est présente, alors l'expression
\(f\paren{Y+\theta\,\Delta\,Y}\) apparaîtra dans le système
d'équations.

Par abus de langage, nous dirons que \(Y+\theta\,\Delta\,Y\) est la
valeur en milieu de pas de \(Y\), que nous noterons \(\milieupas{Y}\)~:
\[
\milieupas{Y} = Y+\theta\,\Delta\,Y
\]
Cette notation s'étend naturellement aux composantes de \(Y\)~:
\[
\milieupas{y_{i}} = y_{i}+\theta\,\Delta\,y_{i}
\]
En particulier, ces notations donnent un sens à l'expression
\(\deriv{\milieupas{y_{i}}}{\Delta\,y_{i}}\) qui est égale à \(\theta\).

De même, nous dirons que \(f\paren{Y+\theta\,\Delta\,Y}\) est la
valeur en milieu de pas de \(f\) et nous la noterons \(\milieupas{f}\)~:
\[
\milieupas{f} = f\paren{Y+\theta\,\Delta\,Y}
\]

\paragraph{Décomposition du vecteur $F$} L'écriture
générique~\eqref{eq:systeme_diff-implicite} résulte de la concaténation
des équations différentielles régissant chaque variable interne \(Y\).
En pratique, le vecteur \(F\) est décomposé en termes associés à chaque
variable interne~:
\[
F=
\begin{pmatrix}
  f_{y_{1}} \\
  \vdots \\
  f_{y_{i}} \\
  \vdots \\
  f_{y_{n}} \\
\end{pmatrix}
\]
Cette écriture est là aussi {\em symbolique} et chaque terme
\(f{y_{j}}\) peut être soit un scalaire soit un tenseur symétrique
d'ordre \(2\).

\mfront{} définit automatiquement les différents termes
\(f_{y_{j}}\). Ces termes se manipulent comme un scalaire si \(y_{i}\)
est un scalaire, et comme une tenseur d'ordre \(2\) si \(y_{i}\) l'est.

Pour des raisons de performances, la variable \(f_{y_{j}}\) est l'image
d'une portion du vecteur \(F\) et toute modification de cette variable
modifie le vecteur \(F\)

Si \(y_{i}\) désigne un tableau de variables internes, \(f_{y_{j}}\)
est une méthode qui prend en argument un indice entier.

\paragraph{Choix du paramètre $\theta$} L'optimum en terme
d'efficacité numérique est obtenu pour \(\theta=\pfrac{1}{2}\) (schéma
dit \og~semi-implicite\fg{}), ce qui est le choix par défaut dans
\mfront{}. Le traitement de mécanismes indépendants du temps
(plasticité ou endommagement) nécessite, pour être physiquement
satisfaisant, d'utiliser la valeur \(1\) (schéma dit \og~purement
implicite\fg{}). La valeur par défaut peut être modifiée par la
directive \mkey{Theta}. À l'exécution, elle peut être modifiée par le
paramètre \varcpp{theta}.

\paragraph{Méthodes mixtes} Le choix du paramètre
\(\theta\) dépend du phénomène traité. Pour des lois mêlant différents
types de phénomènes, certains auteurs choisissent de traiter les
écoulements viscoplastiques par une méthode semi-implicite et les
phénomènes indépendants du temps par une méthode purement implicite. Il
est possible de faire de même dans \mfront{} en prenant certaines
précautions sur le calcul des contraintes.

\paragraph{Initialisation du vecteur $F$}
L'équation~\eqref{eq:F_implicite} montre que la fonction \(f_{y_{i}}\)
comporte comme premier terme l'incrément \(\Delta\,y_{i}\). \mfront{}
initialise donc le vecteur \(F\) à cette valeur avant chaque itération
de la méthode de résolution.

\subsubsection{Initialisation}

Les méthodes de résolution convergent si l'initialisation de la
recherche \(\Delta\, Y^{0}\) est \og~suffisamment proche~\fg{} de la
solution. Dans le cas présent, l'équation~\eqref{eq:explicite-implicite}
conduit à supposer que \(\Delta\, Y\) ne sera qu'une faible correction
de \(Y_{t}\) (pour des pas de temps suffisamment petits), de sorte
qu'une estimation initiale nulle est généralement satisfaisante. Il
s'agit du choix fait par défaut par \mfront{} qui peut être modifié par
l'utilisateur dans un bloc \mkey{Predictor}. Un choix classique est
d'utiliser l'estimation explicite de la solution~:
\[
\Delta\, Y^{0} = \Delta\,t\,G\paren{\debutpas{Y},t}
\]

\subsubsection{Critère d'arrêt}

Les algorithmes présentés s'arrêtent quand la différence entre deux
estimations des incréments des déformations viscoplastiques cumulées
est inférieure à un certain critère \(\varepsilon\)~:
\[
\Frac{1}{\displaystyle\sum_{i=1}^{N}N_{i}}\displaystyle\sum_{i=1}^{N}\left|\Delta_{n}\Delta\,
y_{i}\right|<\varepsilon
\]
où \(\Delta_{n}\Delta\,y_{i}\) désigne la différence entre
l'estimation de l'incrément de la i\ieme{} variable à l'étape \(n+1\) et
de son estimation à l'étape \(n\) et \(N_{i}\) est le nombre de
composantes de la variable interne \(y_{i}\)\footnote{Pour les tenseurs,
  \(\left|s\right|\) représente la somme des valeurs absolues des termes
  du tenseur (les composantes extra-diagonales ne sont comptées qu'une
  fois).}.

\paragraph{Notes} Cette écriture du critère d'arrêt
montre qu'il n'est valide que si les différentes variables internes sont
du même ordre de grandeur. Ceci est également nécessaire pour le bon
comportement numérique de la méthode de \nom{Newton}. Dans le cas
contraire, il est possible de normaliser les variables par la méthode
{\tt setNormalisationFactor}.

\subsubsection{Avantages et inconvénients des méthodes implicites} Les
méthodes implicites sont généralement opposées aux méthodes
d'intégration plus ou moins basées sur des algorithmes de type
\nom{Runge-Kutta}. Ces dernières méthodes présentent plusieurs
avantages~:
\begin{itemize}
\item un ordre de convergence qui peut être élevé~;
\item une maîtrise de la précision des résultats par
  \og~correction/prédiction\fg{}. Cette maîtrise permet de
  sous-découper localement le pas de temps~;
\item une facilité d'implantation par rapport aux méthodes implicites,
  ces méthodes ne nécessitant pas le calcul de matrice jacobienne
  notamment.
\end{itemize}

\mfront{} propose un analyseur générique basé sur les algorithmes de
type \nom{Runge-Kutta} (par défaut un correcteur-prédicteur \(5/4\) à
pas de temps adaptatif qui allie un ordre de convergence élevé à une
efficacité numérique acceptable) décrit à la section~\ref{sec:RK}.

Le principal défaut de ces méthodes est le nombre d'appels nécessaires
à la fonction \(G\) du système~\eqref{eq:systeme_diff-implicite}, en
particulier si une stratégie à pas de temps adaptatif est utilisée.

En comparaison, les méthodes implicites sont réputées beaucoup plus
stables numériquement et ne nécessitent qu'un nombre faible d'appels à
la fonction \(G\). De plus, les méthodes implicites sont
particulièrement adaptées aux lois de comportement indépendantes du
temps (endommagement ou plasticité) où le système différentiel peut être
remplacé par la nullité du critère plastique en fin de pas de temps, ce
qui est précisément la condition {\em physique} à vérifier\footnote{Nous
  avons traité le cas de la plasticité isotrope au
  paragraphe~\ref{sec:mfront:isotropic-solver}.}. Cette particularité
explique pourquoi les méthodes implicites sont toujours à préférer aux
méthodes de type \nom{Runge-Kutta} dans le cas de lois de comportement
indépendantes du temps, même si elles présentent sur le papier des
ordres de convergence moins élevés.

L'application des méthodes implicites aux lois de comportement
dépendantes du temps (écoulements viscoplastiques notamment) est plus
délicate en raison justement de ces ordres de convergences plus
faibles (que les méthodes de type \nom{Runge-Kutta})\footnote{Il ne
  faut pas oublier qu'il ne sert généralement à rien d'utiliser des
  pas de temps trop grands. En effet, l'algorithme de résolution
  global (en quasi-statique) utilise une résolution implicite~: sur un
  pas de temps, il faut que l'hypothèse qui permet la linéarisation de
  l'évolution mécanique de la structure soit vérifiée.}.

En pratique, nous avons {\em toujours} constaté des résultats
extrêmement satisfaisants tant en termes de précision qu'en termes de
performance numérique avec les méthodes implicites, quel que soit le
type de loi de comportement.

\begin{table}[htbp]
  \centering
  \begin{tabular}{|c|c|c|}
    \hline Méthode d'intégration & Temps CPU & Nombre
    d'itérations \\
    \hline \hline Implicite & 5 m 26 s & 1817 itérations \\
    \hline \nom{Runge-Kutta} \(5/4\) & 32 m 21 s & 17493
    itérations \\
    \hline
  \end{tabular}
  \caption{Temps de calculs obtenus sur une éprouvette
    entaillée en traction modélisée en $2D\paren{r,z}$. Calcul réalisé
    par \nom{J.-M. Proix} sur la version $11$ de code \aster{}. La loi
    de comportement locale est la loi de
    \nom{HayHurst}~\cite{proix_comportement_2012}.
    Cette comparaison a été fournie par J.-M. Proix (EDF/AMA).}
  \label{tab:mfront:NR:comparison:methode}
\end{table}

Enfin, les méthodes implicites introduisent naturellement la
matrice tangente cohérente, ce qui permet de réduire considérablement le
nombre d'itérations nécessaires à la méthode globale. Ceci est illustré
sur la figure~\ref{tab:mfront:NR:comparison:methode}.

Ceci explique que la plupart des analyseurs syntaxiques
disponibles dans \mfront{} soient basés sur ces méthodes.

Les deux principales faiblesses des méthodes implicites restent~:
\begin{itemize}
  \item la nécessité de calculer la matrice jacobienne ou d'utiliser des
  méthodes alternatives moins efficaces~;
  \item l'absence de contrôle de l'erreur.
\end{itemize}

\subsection{Méthode de \nom{Newton-Raphson}}
\label{sec:NR}

Une méthode généralement utilisée pour trouver un zéro d'une fonction
est la méthode itérative de \nom{Newton-Raphson}.

Connaissant une estimation \(\Delta\,Y^{n}\) de la solution, la fonction
\(F\) peut être approximée dans un voisinage \(\Delta\,Y^{n}\) par~:
\[
F\paren{\Delta\,Y}=F\paren{\Delta\,Y^{n}}+\left.\derivtot{F}{\Delta\,Y}\right|_{\Delta
 Y=\Delta\,Y^{n}}\,\paren{\Delta\,Y-\Delta\,Y^{n}}
\]
où \(\derivtot{F}{\Delta\,Y}\) est la matrice jacobienne de la fonction
\(F\), qui sera notée \(J\) dans la suite.

Une nouvelle estimation \(\Delta\,Y^{n+1}\) est choisie en annulant
cette approximation de \(F\), ce qui conduit à l'expression suivante~:
\begin{equation}
  \label{eq:NR}
  \Delta\, Y^{n+1}=\Delta\,Y^{n}-\left[J^{-1}\paren{\Delta\,Y^{n}}\right]\,F\paren{Y^{n}}
\end{equation}

Dans cette équation, l'écriture \(\left[J^{-1}\paren{\Delta\,Y^{n}}\right]\,F\paren{Y^{n}}\) représente la
solution \(\Delta^{2}\,Y\) du système linéaire~:
\[
\left[J\paren{\Delta\,Y^{n}}\right]\,\Delta^{2}\,Y = F\paren{Y^{n}}
\]
%
%\begin{figure}[htbp]
%  \centering
%  \input{@top_srcdir@/docs/mfront/Images/Methode_Newton.pstricks}
%  \caption[Illustration de la méthode de \nom{Newton} en $1D$ utilisée
%  pour la résolution de l'équation $f\paren{y}=0$.]{Illustration de la
%    méthode de \nom{Newton} en $1D$ utilisée pour la résolution de
%    l'équation $f\paren{y}=0$. Les estimations successives de la
%    solution sont données par~:
%    $y_{n+1}=y_{n}-\Frac{f\paren{y_{n}}}{f\primeexp\paren{y_{n}}}$.}
%  \label{fig:illustration_newton_raphson}
%\end{figure}
%
%\paragraph{Illustration en \(1D\)} La
%figure~\ref{fig:illustration_newton_raphson} illustre la méthode de
%\nom{Newton} en $1D$.

\paragraph{Matrice jacobienne}
La matrice jacobienne s'exprime ainsi~:
\begin{equation}
  \label{eq:J_implicite}
  J = \derivtot{F}{\Delta\,Y} = I-\theta\Delta\,t\,\derivtot{G}{Y}\paren{\debutpas{Y}+\theta\Delta{Y},t+\theta\,\Delta t}
\end{equation}
où \(I\) est la matrice identité.

{\em Le calcul de la matrice jacobienne est généralement le point le
  plus difficile de la méthode implicite.}

\paragraph{Calcul de la matrice jacobienne par blocs} En pratique, la
matrice \(J\) est calculée par blocs et prend la forme suivante~:
\[
J = \deriv{F}{\Delta\,Y} =
\begin{pmatrix}
  \deriv{f_{y_{1}}}{\Delta\,y_{1}} & \ldots & \ldots & \ldots & \ldots \\
  \vdots & \vdots & \vdots & \vdots & \vdots \\
  \vdots & \vdots & \deriv{f_{y_{i}}}{\Delta\,y_{j}} & \vdots & \vdots \\
  \vdots & \vdots & \vdots & \vdots & \vdots \\
  \ldots & \ldots & \ldots & \ldots & \deriv{f_{y_{N}}}{\Delta\,y_{N}} \\
\end{pmatrix}
\]

Les différents termes \(\deriv{f_{y_{i}}}{\Delta\,y_{j}}\) sont
automatiquement définies par \mfront{}. La nature mathématique de ces
termes dépend de la nature des variables \(y_{i}\) et \(y_{j}\). Elles
peuvent être manipulées par les opérations mathématiques
usuelles. Pour des raisons de performances, ces variables sont des
images de la matrice jacobienne~: en les manipulant, on modifie
directement la matrice jacobienne.

Ainsi, la dérivée d'une variable \varcpp{a} par rapport à la
variable \varcpp{b} sera accessible par \varcpp{dfa\_ddb}.

Dans le cas particulier des tableaux de variables internes,
\varcpp{dfa\_ddb} est une méthode qui prend~:
\begin{itemize}
\item un unique argument entier si soit \varcpp{a} ou soit \varcpp{b}
  désigne un tableau de variables internes, mais pas les deux en mêmes
  temps~;
  \item deux arguments entiers si \varcpp{a} et \varcpp{b}
  désignent deux tableaux de variables internes~;
\end{itemize}
Ces méthodes renvoient des objets spéciaux, qui agissent comme des
scalaires ou des tenseurs d'ordre \(2\) ou \(4\) suivant les cas. Ces
objets permettent un accès aux valeurs de la jacobienne qui peut être,
en fonction des capacités du compilateur utilisé, moins performant que
dans le cas des variables ordinaires.

\paragraph{Comparaison de la matrice jacobienne
  fournie à une approximation numérique} Le calcul de la matrice
jacobienne est, à juste titre, jugé délicat et une erreur peut avoir un
impact significatif sur les performances ou la robustesse de
l'intégration et dans le meilleur des cas une non
convergence\footnote{En effet, dans ce cas, l'erreur apparaît alors que
  sinon le seul symptôme est une convergence lente qui peut passer
  inaperçue si l'on y prend garde.}.

Pour détecter les erreurs dans le calcul du jacobien, nous avons
introduit le mot clé \mkey{Compare\-To\-Numerical\-Jacobian}. Ce mot
clé déclenche un processus de vérification par comparaison du jacobien
fournit par l'utilisateur à un jacobien évalué numériquement par
perturbations. Cette vérification se fait par étapes~:
\begin{itemize}
  \item la première étape est le calcul d'un jacobien numérique
  par perturbations. Pour gagner en précision, une différence finie
  centrée est utilisée\footnote{En \(1D\), la dérivée numérique d'une
    fonction \(f\paren{x}\) au point \(x\) est évaluée ainsi~:
    \[
    \Frac{f\paren{x+\epsilon}-f\paren{x-\epsilon}}{2\,\epsilon}
    \]
    La perturbation est choisie égale à la valeur du critère
    d'arrêt de l'algorithme implicite. }~;
  \item la seconde étape est de comparer les termes du jacobien.
  La comparaison se fait bloc par bloc~: la différence entre la dérivée
  \(\deriv{f_{y_{i}}}{\Delta\,y_{j}}\) calculée par l'utilisateur et son
  approximation numérique \(\Frac{\Delta\,f_{y_{i}}}{\Delta\,\Delta\,y_{j}}\)
  est évaluée ainsi~:
  \[
  \Frac{1}{N_{y_{i}}\,N_{y_{j}}}\norm{\deriv{f_{y_{i}}}{\Delta\,y_{j}}-\Frac{\Delta\,f_{y_{i}}}{\Delta\,\Delta\,y_{j}}}
  \]
  où \(N_{y_{i}}\) est le nombre de composantes de la variables
  \(y_{i}\) et où \(\norm{a}\) est la somme des valeurs absolues des
  composantes de \(a\). Cette différence est comparée à un certain
  critère. Ce critère est par défaut égal au critère de convergence de
  l'algorithme implicite, mais l'utilisateur peut en préciser en autre
  en utilisant le mot clé {\tt @Jacobian\-Comparison\-Criterium}.
  L'utilisateur peut également utiliser le paramètre
  \varcpp{jacobian\-Comparison\-Criterium} pour modifier ce critère à
  l'exécution. Si la différence est supérieure au critère, un message
  donnant la valeur de la différence, le bloc considéré et son
  approximation numérique est affiché sur la sortie standard.
\end{itemize}

\begin{table}
  \centering
  \begin{tabular}{|c|c|c|}
    \cline{2-3}
    \multicolumn{1}{c|}{} &
    \begin{minipage}{4cm}
      Nombre de cycles
    \end{minipage} &
    \begin{minipage}{4cm}
      Rapport à la référence
    \end{minipage} \\
    \hline
    Référence             & 4\,776\,120. & 1 \\
    \hline
    \begin{minipage}[p]{4cm}
    Calcul du jacobien numérique et comparaison
    \end{minipage}
    & 26\,944\,564 & 5,65\\
    \hline
  \end{tabular}
  \caption{Impact du calcul de jacobien numérique et de sa comparaison à
    celui fourni par l'utilisateur pour une loi de comportement simple.}
  \label{tab:NR:nJ}
\end{table}

Le coût du calcul du jacobien numérique est souvent important. Nous
pouvons le constater au tableau~\ref{tab:NR:nJ} où nous avons analysé
le nombre de cycles processeur pour réaliser l'intégration d'une loi
de \textsc{Norton} lors d'un essai de traction uniaxial\footnote{Cet
  essai est décrit plus en détails en
  annexe~\ref{QNR:ComparaisonNumeriques}.}. Pour le cas traité, le
coût de l'intégration est multiplié par \(5,65\). Il faut donc prendre
garde à retirer cette option pour les versions de production des lois
de comportement.

\paragraph{Initialisation de la matrice jacobienne}
L'expression~\eqref{eq:J_implicite} montre que le premier terme de la
matrice jacobienne est la matrice identité. Cette écriture nous a
conduit, dans \mfront{} à initialiser, à chaque itération, la matrice
jacobienne à la matrice identité.

\subsection{Méthodes alternatives}

Insistons, le calcul de la matrice jacobienne est la
principale difficulté de la méthode de \nom{Newton-Raphson}.

Afin de dispenser l'utilisateur de fournir la matrice jacobienne, nous
avons envisagé d'approximer le jacobien par différentiation
numérique. Un jacobien correctement évalué (par une différence finie
d'ordre \(2\)) conduit à une bonne convergence de l'algorithme, mais a
un coût numérique important\footnote{Voir, en
  annexe~\ref{QNR:ComparaisonNumeriques}, les tableaux~\ref{tab:QNR:1}
  et~\ref{tab:QNR:2}}. Cette méthode est accessible par l'algorithme
\texttt{Newton\-Raphson\textunderscore{}\-Numerical\-Jacobian}.

Ce constat est général et de nombreux auteurs ont proposé des
algorithmes qui ne se basent pas sur une différentiation numérique
(sauf éventuellement pour une estimation initiale ou en cas de non
convergence).

L'idée de ces méthodes est de modifier l'algorithme~\eqref{eq:NR} en
substituant à la matrice jacobienne une matrice \(\underset{\sim}{J}\)
bien choisie~:
\begin{equation}
  \label{eq:QNR}
  Y_{n}=Y_{n-1}-\underset{\sim}{J}^{-1}\,.\,F\paren{Y_{n-1}}
\end{equation}

Parmi ces méthodes, les algorithmes de \textsc{Broyden} sont les plus utilisés.

\subsubsection{Les algorithmes de \textsc{Broyden}}

Nous décrivons dans ce paragraphe deux algorithmes dû à \textsc{Broyden}.

\paragraph{Premier algorithme de \textsc{Broyden}} Un des algorithmes
les plus utilisés est celui de \textsc{Broyden} qui met à jour au cours
des itérations une approximation \(\underset{\sim}{J}_{n}\) de la
matrice jacobienne à partir de son expression
\(\underset{\sim}{J}_{n-1}\) à l'itération précédente~:
\begin{equation}
  \label{eq:Broyden} \underset{\sim}{J}_{n} =
  \underset{\sim}{J}_{n-1}+\Frac{\Delta\,F_{n}-\underset{\sim}{J}_{n-1}\,.\,\Delta\,Y_{n}}{\norm{\Delta\,Y_{n}}^{2}}\otimes\Delta\,Y_{n}
\end{equation}
où~:
\begin{minipage}[t]{0.8\linewidth}
  \begin{itemize}
    \item \(\Delta\,Y_{n}=Y_{n}-Y_{n-1}\)~;
    \item
    \(\Delta\,F_{n}=F\paren{Y_{n}}-F\paren{Y_{n-1}}\)~;
    \item \(v_{1}\,\otimes\,v_{2}\) est le produit tensoriel de deux
    vecteurs~;
    \item \(\norm{v}\) est la norme euclidienne du vecteur \(v\).
  \end{itemize}
\end{minipage}

\paragraph{Calcul exact de certains blocs}
Il est intéressant de remarquer qu'il est toujours possible de
modifier la matrice jacobienne par blocs~: si certains termes sont
connus et/ou faciles à calculer, il peut être intéressant de les
injecter dans la matrice jacobienne en écrasant les modifications
fournies par l'algorithme de \textsc{Broyden}.

\paragraph{Initialisation du jacobien} Il est nécessaire de fournir
une estimation initiale \(\underset{\sim}{J}_{0}\) du jacobien. Deux
choix ont été regardés dans cette note~:
\begin{itemize}
  \item \(\underset{\sim}{J}_{0}\) est choisi égal à l'identité.
  L'expression~\eqref{eq:J_implicite} montre que pour des pas de temps
  raisonnablement petits, l'identité constitue une bonne approximation
  de la matrice jacobienne~;
  \item \(\underset{\sim}{J}_{0}\) est calculé par perturbations du
  système initial.
\end{itemize}

\paragraph{Second algorithme de \textsc{Broyden}}
Une variante de l'algorithme de \textsc{Broyden} met à jour une
approximation \(\underset{\sim}{J}_{n}^{-1}\) de l'inverse de la matrice
jacobienne par la relation suivante~:
\begin{equation}
  \label{eq:Broyden2}
  \underset{\sim}{J}_{n}^{-1} =
  \underset{\sim}{J}_{n-1}^{-1}+\Frac{\Delta\,Y_{n}-\underset{\sim}{J}_{n-1}^{-1}\,.\,\Delta\,F_{n}}{\Delta\,Y_{n}\,\mid\,\underset{\sim}{J}_{n-1}^{-1}\,.\,\Delta\,F_{n}}\,\otimes\,\paren{\Delta\,Y_{n}\,\mid\,\underset{\sim}{J}_{n-1}^{-1}}
\end{equation}

Cette variante est disponible dans \mfront{} sous le nom de
\texttt{Broyden2}.

Dans ce cas, il n'est pas possible de fournir des blocs du jacobien.

\subsection{Amélioration de la robustesse des algorithmes implicites}

Différentes méthodes permettant d'améliorer la robustesse des
algorithmes ont été testées. Nous les décrivons maintenant bien que leur
utilisation conduise à des résultats mitigés.

Pour des raisons d'efficacité du code généré, il a été décidé que le
support d'une méthode d'accélération était une décision prise au
moment de la génération du code. Les méthodes présentées présentent
des paramètres permettant de fixer un nombre minimal d'itérations
avant que la méthode s'active. En donnant une valeur suffisamment
grande à ce paramètre, la méthode ne s'activera jamais, ce qui ne veut
pas dire qu'elle n'aura aucun coût : par exemple, les méthodes
présentées utilisent des valeurs des itérées précédentes qui sont
toujours stockées.

Enfin, les paramètres associés à ces méthode ne seront disponibles que
si elles sont activées.

\subsubsection{Méthode d'accélération utilisé par
  les algorithmes de résolution globaux de \castem{}}

L'idée de cette méthode d'accélération est de trouver une nouvelle
estimation de la solution à partir des itérations précédentes. Soient
\(F_{n-2}\), \(F_{n-1}\) et \(F_{n}\) les résidus associés aux trois
dernières estimations \(\Delta\, Y_{n-2}\), \(\Delta\, Y_{n-1}\) et
\(\Delta\, Y_{n}\) de la solution.

Il est possible de construire à partir de ces trois résidus un
hyperplan. L'idée est d'exprimer la projection du vecteur nul (qui
correspond à la solution recherchée) sur cet hyperplan comme une
combinaison linéaire de ces trois résidus et de déduire une nouvelle
estimation de la solution en appliquant la même combinaison linéaire aux
estimations \(\Delta\, Y_{n-2}\), \(\Delta\, Y_{n-1}\) et \(\Delta\,
Y_{n}\).

Dans certains cas difficiles, cette méthode peut stabiliser les
calculs, mais elle empêche la convergence quadratique de la méthode de
\nom{Newton}. Son intérêt apparaît aujourd'hui très limité.

La directive \mkey{UseAcceleration} est utilisée pour activer cette
méthode d'accélération. Elle est suivie soit du mot clé \varcpp{true}
soit du mot clé \varcpp{false} et d'un point-virgule.

La directive \mkey{AccelerationTrigger} permet de préciser à partir de
quelle itération la méthode d'accélération est utilisée. La valeur
fournie doit être supérieure à \(3\) pour que la méthode ait un sens. La
valeur fournie peut être modifiée par le paramètre
\varcpp{acceleration\-Trigger}.

La directive \mkey{AccelerationPeriod} permet de préciser le nombre
d'itérations séparant deux emplois successifs de cette méthode. La
valeur fournie peut être modifiée par le paramètre
\varcpp{acceleration\-Period}.

\subsubsection{Méthode de relaxation}

Une méthode de relaxation a été introduite pour améliorer la
convergence de certains cas difficiles où des oscillations lors de la
recherche de solution peuvent apparaître. Soient \(\Delta\, Y_{n-1}\) et
\(\Delta\, Y_{n}\) les deux dernières estimations de la solution. La
relaxation consiste à prendre comme nouvelle estimation~:
\[
\paren{1-w}\,\Delta\, Y_{n-1}+w\,\Delta\, Y_{n}
\]
où \(w\) est un coefficient de relaxation.

La directive \mkey{RelaxationTrigger} permet de préciser à partir de
quelle itération la méthode d'accélération est utilisée. La valeur
fournie doit être supérieure à \(3\) pour que la méthode ait un sens. La
valeur fournie peut être modifiée par le paramètre
\varcpp{relaxation\-Trigger}.

La directive \mkey{RelaxationCoefficient} permet de préciser le
coefficient \(w\). La valeur fournie peut être modifiée à l'exécution
par le paramètre \varcpp{relaxation\-Coefficient}.

% \subsubsection{Autres pistes}

% D'autres pistes d'amélioration peuvent être suivies. Nous renvoyons à
% l'annexe~\ref{sec:mfront:implicit:robustness} pour la description d'un
% algorithme dû à \nom{Powell} qui semble particulièrement intéressant.

\subsection{Calcul de la matrice tangente cohérente}

L'utilisateur peut calculer la matrice tangente cohérente
dans un bloc introduit par la directive \mkey{TangentOperator}.

Cette directive est appelée si l'intégration de la loi a été un
succès, avant la mise à jour des variables internes et la mise à jour de
contraintes.

\paragraph{Symétrie de la matrice tangente} Par défaut, la matrice
tangente est supposée non symétrique. L'utilisateur peut explicitement
dire que la matrice est symétrique par la directive
\mkeyb{IsTangentOperatorSymmetric}{Is\-Tangent\-Operator\-Symmetric}
qui est suivie d'une des valeurs booléennes \varcpp{true} ou
\varcpp{false}\footnote{Si l'opérateur tangent est symétrique,
  l'interface \aster{} fait l'économie d'une transposition de la
  matrice (cette transposition est nécessaire pour la compatibilité
  avec le {\tt fortran}.}.

\subsubsection{Une façon générique de calculer de la matrice tangente
  cohérente}
\label{sec:mfront:implicit:tangentOperator}

La matrice tangente cohérente est la dérivée d'incrément de
contraintes \(\tsigma\) par rapport à l'incrément de déformations
totales \(\tepsilonto\). Cet incrément était considéré jusqu'à présent
comme un paramètre du système d'équations implicites. Il est considéré
dans ce paragraphe comme la variable principale.

Pour commencer, nous supposerons que le tenseur d'élasticité est une
des variables internes de la loi. Plus loin, nous imposerons que cette
variable soit la première dans le système implicite. Ces conditions sont
imposées si on utilise l'analyseur \texttt{Implicit}.

Nous supposerons ensuite que la contrainte ne dépende que du
tenseur d'élasticité (éventuellement de manière non linéaire). Nous
reviendrons sur cette hypothèse plus tard. Nous en tirons que~:
\[
\deriv{\Delta\,\tsigma}{\Delta\,\tepsilonto}
=
\left.\deriv{\tsigma}{\tepsilonel}\right|_{\tepsilonel+\Delta\,\tepsilonel}\,\colon\,\deriv{\Delta\,\tepsilonel}{\Delta\,\tepsilonto}
\]

Le problème du calcul de la matrice tangente cohérente se ramène donc
au calcul de la dérivée
\(\deriv{\Delta\,\tepsilonel}{\Delta\,\tepsilonto}\).

Nous montrons dans ce paragraphe que si la matrice jacobienne du
système est connue, alors il existe une façon générique de calculer ce
terme dans de très nombreux cas.

\paragraph{Rappels mathématiques} Avant de rentrer
dans les détails de la méthode, il est utile de faire quelques rappels
mathématiques. Les incréments des variables internes \(\Delta\,Y\) sont
reliées à l'incrément de déformations totales par la relation
implicite~:
\[
F\paren{\Delta\,Y\paren{\Delta\tepsilonto},\Delta\,\tepsilonto}
= 0
\]
Si \(\Delta\tepsilonto\) varie, cette relation reste vraie. Par
différentiation, nous obtenons~:
\[
\dtot\,F = \deriv{F}{\Delta\,Y}\,\dtot\,\Delta\,Y+\deriv{F}{\Delta\,\tepsilonto}\,\dtot\,\Delta\,\tepsilonto=0
\]

Le terme \(\deriv{F}{\Delta\,Y}\) est la jacobienne \(J\) du système
qui est connue après la résolution.

Pour déduire \(\deriv{\Delta\,\tsigma}{\Delta\,\tepsilonto}\) de cette
équation, nous allons faire une hypothèse qui permet de calculer
explicitement \(\deriv{F}{\Delta\,\tepsilonto}\).

\paragraph{Hypothèse} Nous supposons que l'incrément de déformation
totales \(\Delta\tepsilonto\) n'apparaît que dans l'équation relative
aux déformations élastiques par la relation de partition des
déformations, qui, une fois discrétisée, donne~:
\[
\Delta\,\tepsilonel+\sum_{i}\Delta\,\tepsilonan_{i}-\Delta\,\tepsilonto=0
\]

La dérivée de \(\deriv{F}{\Delta\,\tepsilonto}\) est alors évidente et
le vecteur
\(\deriv{F}{\Delta\,\tepsilonto}\,\dtot\,\Delta\,\tepsilonto\) est égal
à~:
\[
\deriv{F}{\Delta\,\tepsilonto}\,\dtot\,\Delta\,\tepsilonto = -
\begin{pmatrix}
  \dtot\,\Delta\tepsilonto \\
  \vdots \\
  0 \\
  \vdots \\
  0\\
\end{pmatrix}
\]

Nous en déduisons que~:
\[
\dtot\,\Delta\,\tepsilonel=J^{-1}_{\tepsilonel}\,\colon\,\dtot\,\Delta\tepsilonto
\]
où \(J^{-1}_{\tepsilonel}\) est la partie supérieure gauche de l'inverse de
la jacobienne.

Finalement, nous obtenons~:
\[
\deriv{\Delta\,\tsigma}{\Delta\,\tepsilonto} =  \left.\deriv{\tsigma}{\tepsilonel}\right|_{\tepsilonel+\Delta\,\tepsilonel}\,\colon\, J^{-1}_{\tepsilonel}
\]

La matrice \(J^{-1}_{\tepsilonel}\) peut être calculée par la méthode
\varcpp{getPartialJacobianInvert} à qui l'on doit passer un tenseur
d'ordre \(4\). Cette méthode ne doit être utilisée que dans le bloc
\mkey{TangentOperator}.

Ce cas explique pourquoi la déformation élastique est toujours définie
comme la première variable interne dans le cas de l'analyseur
\texttt{Implicit}.

\paragraph{Généralisation} La méthode
précédente se généralise dans le cas où la contrainte dépend de
plusieurs variables internes et éventuellement de la déformation totale.
En effet, il suffit d'écrire~:
\[
\deriv{\Delta\,\tsigma}{\Delta\,\tepsilonto}
=
\sum_{i}\left.\deriv{\tsigma}{y_{i}}\right|_{y_{i}+\Delta\,y_{i}}\,\colon\,\deriv{\Delta\,y_{i}}{\Delta\,\tepsilonto}
+
\left.\deriv{\tsigma}{\tepsilonto}\right|_{\tepsilonto+\Delta\,\tepsilonto}
\]
Les mêmes considérations que précédemment montrent que les
différentes dérivées \(\deriv{\Delta\,y_{i}}{\Delta\,\tepsilonto}\) se
lisent sur les \(N_{\tepsilonto}\) premières colonnes de l'inverse de la
matrice jacobienne, où \(N_{\tepsilonto}\) est le nombre de composantes
du tenseur \(\tepsilonto\).

Il est possible de récupérer ces dérivées par la méthode
\varcpp{get\-Partial\-Jacobian\-Invert} à qui l'on doit passer
successivement des variables qui vont contenir les dérivées des
variables d'intérêt dans {\em l'ordre de déclaration des
  variables}. Ces variables doivent être des tenseurs d'ordre \(2\)
dans le cas de variables scalaires et des tenseurs d'ordre \(4\) dans
le cas de variables tensorielles.

Dans le cas des tableaux de variables internes, il faut passer des
vecteurs de tailles finies du type adéquat\footnote{On utilisera la
  classe \tfel{tvector} de \TFEL{} qui paramétrée par sa taille et le
  type d'objet contenu.}.

Insistons~: si l'on a besoin que des dérivées de deux
premières variables internes (par exemple la déformation élastique et la
variable d'endommagement) alors il suffit de passer deux arguments à la
méthode \varcpp{get\-Partial\-Jacobian\-Invert}, on ne calculera donc
pas les autres dérivées pour rien. À l'utilisateur de déclarer les
variables dans le bon ordre !

\subsubsection{Une méthode alternative de calcul de la matrice
  tangente cohérente}

Il existe une autre méthode de calcul de la matrice tangente cohérente
basée sur un partitionnement de la jacobienne et non de son
inverse~\cite{proix_integration_2012}.

Pour cela, notons \(z_{i}\) les variables internes autres que les
déformations élastiques. Nous avons, avec les mêmes hypothèses que
précédemment~:
\[
J\,.\,
\begin{pmatrix}
  \dtot\,\Delta\,\tepsilonel \\
  \dtot\,\Delta\, Z
\end{pmatrix}
-
\begin{pmatrix}
  \dtot\,\Delta\tepsilonto \\
  0 \\
\end{pmatrix}
= 0
\]

En découpant la matrice jacobienne en quatres blocs~:
\[
J=
\begin{pmatrix}
  J_{0} & J_{1} \\
  J_{2} & J_{3} \\
\end{pmatrix}
\]
on obtient deux systèmes~:
\[
\left\{
\begin{aligned}
J_{0}\,\dtot\,\Delta\,\tepsilonel+J_{1}\, \dtot\,\Delta\, Z &= \dtot\,\Delta\tepsilonto \\
J_{2}\,\dtot\,\Delta\,\tepsilonel+J_{3}\, \dtot\,\Delta\, Z &= 0\\
\end{aligned}
\right.
\]
Le second système conduit à~:
\[
\dtot\,\Delta\, Z = -J_{3}^{-1}\,J_{2}\,\dtot\,\Delta\,\tepsilonel
\]
Par report dans le premier système, nous obtenons~:
\[
\dtot\,\Delta\,\tepsilonel = \paren{J_{0}-J_{1}\,J_{3}^{-1}\,J_{2}}^{-1}\,\dtot\,\Delta\tepsilonto
\]
et finalement~:
\[
\deriv{\Delta\,\tepsilonel}{\Delta\tepsilonto}=\paren{J_{0}-J_{1}\,J_{3}^{-1}\,J_{2}}^{-1}
\]

Cette méthode peut être avantageuse si le nombre de variables
internes est petit.

\subsection{La directive \texttt{@ComputeStress}}

La directive \mkey{ComputeStress} permet de calculer les contraintes aux
différentes itérations de la méthode et en fin de calcul. Les contraintes
sont supposées fonction des valeurs actuelles des variables internes.

\paragraph{Conventions spécifiques aux différentes itérations de la
  méthode} Dans la méthode \mkey{ComputeStress}, les conventions
suivantes s'appliquent~:
\begin{itemize}
  \item {\tt sig} représente la contrainte qu'il faut calculer~;
  \item {\tt eto} représente la déformation totale en milieu de pas (en
  \(t+\theta\,\Delta\,t\))~;
  \item {\tt deto} représente l'incrément de la déformation totale
  (constante sur le pas)~;
  \item {\tt T} représente la valeur de la température en milieu de
  pas (en
  \(t+\theta\,\Delta\,t\))~;
  \item {\tt dT} représente la variation de température (constante sur
  le pas)~;
  \item pour toute variable interne {\tt Y}, {\tt Y} représente son
  estimation actuelle en milieu de pas (en
  \(t+\theta\,\Delta\,t\))~;
  \item pour toute variable interne auxiliaire {\tt Y}, {\tt Y}
  représente sa valeur en {\em début} de pas de temps~;
  \item pour toute variable externe {\tt V}, {\tt V} représente sa
  valeur en milieu de pas (en
  \(t+\theta\,\Delta\,t\))~;
  \item pour toute variable externe {\tt V}, {\tt dV} représente son
  incrément sur le pas de temps (constante sur le pas).
\end{itemize}

\paragraph{Conventions spécifiques en fin d'intégration} Dans la méthode
\mkey{ComputeStress}, les conventions suivantes s'appliquent~:
\begin{itemize}
  \item {\tt sig} représente la contrainte qu'il faut calculer~;
  \item {\tt eto} représente la déformation totale en fin de pas~;
  \item {\tt deto} représente la vitesse de déformation totale
  (constante sur le pas)~;
  \item {\tt T} représente la valeur de la température en fin de pas~;
  \item {\tt dT} représente la vitesse de changement de température
  (constante sur le pas)~;
  \item pour toute variable interne {\tt Y}, {\tt Y} représente sa
  valeur en fin de pas~;
  \item pour toute variable interne auxiliaire {\tt Y}, {\tt Y}
  représente sa valeur en {\em début} de pas de temps~;
  \item pour toute variable externe {\tt V}, {\tt V} représente sa
  valeur en fin de pas~;
  \item pour toute variable externe {\tt V}, {\tt dV} représente sa
  vitesse de variation sur le pas de temps (constante sur le pas).
\end{itemize}

\subsection{La directive \texttt{@Integrator}}

La directive \mkey{Integrator} permet de construire le système implicite
à résoudre.

\paragraph{Conventions spécifiques} Les conventions suivantes
s'appliquent~:
\begin{itemize}
  \item {\tt sig} représente la contrainte actualisée (calculée par le
  code fourni après la directive
  \mkeyb{ComputeStress}{Compute\-Stress})~;
  \item {\tt eto} représente la déformation totale en début de pas~;
  \item {\tt deto} représente la vitesse de déformation totale
  (constante sur le pas)~;
  \item {\tt T} représente la valeur de la température en début de pas~;
  \item {\tt dT} représente la vitesse de changement de température
  (constante sur le pas)~;
  \item pour toute variable interne {\tt Y}, {\tt Y} représente sa
  valeur en début de pas~;
  \item pour toute variable interne {\tt Y}, {\tt Y} représente
  l'estimation courante de l'incrément de cette variable sur le pas~;
  \item pour toute variable interne {\tt Y}, {\tt fY} représente
  l'équation implicite associée à cette variable~;
  \item pour toute variable interne auxiliaire {\tt Y}, {\tt Y}
  représente sa valeur en début de pas~;
  \item pour toute variable externe {\tt V}, {\tt V} représente sa
  valeur en début de pas~;
  \item pour toute variable externe {\tt V}, {\tt dV} représente sa
  vitesse de variation sur le pas de temps (constante sur le pas).
\end{itemize}

Si l'algorithme de \nom{Newton} ou le premier algorithme de
\nom{Broyden} est utilisé, pour tout couple de variables internes {\tt
  Yi} et {\tt Yj}, {\tt dfYi\textunderscore{}ddYj} représente le terme du jacobien
\(\deriv{f_{Y_{i}}}{Y_{j}}\)~;

\subsection{Mise à jour des variables auxiliaires}

Les variables auxiliaires sont mises à jour après la mise à jour des
variables internes et la mise à jour des contraintes, une fois
l'algorithme convergé.

À ce moment, les valeurs externes (déformations totales et températures
incluses) n'ont pas été mises à jour.

\subsection{Ordre d'évaluation des blocs définis par l'utilisateur}


\begin{figure}[htbp]
  \centering
  \includegraphics{@top_srcdir@/docs/mfront/Images/ImplicitDSL.eps}
  \caption{Ordre d'évaluation des différents blocs de code fournis par
    l'utilisateur dans le cas d'une intégration implicite.}
  \label{fig:ImplicitOrder}
\end{figure}

La figure~\ref{fig:ImplicitOrder} montre l'ordre d'appel des
différents blocs de code définis par l'utilisateur.

\subsection{Paramètres automatiquement définis}

L'analyseur implicite définit automatiquement différents
paramètres~:
\begin{itemize}
  \item \varcpp{epsilon}, la valeur du critère d'arrêt de
  l'algorithme implicite. La valeur par défaut de ce paramètre est de
  \(10^{-8}\) et peut être modifiée par la directive \mkey{Epsilon}. Ce
  paramètre doit avoir une valeur strictement positive~;
  \item \varcpp{theta}, la valeur de \(\theta\) utilisé par
  l'algorithme implicite. La valeur par défaut de ce paramètre est de
  \(0.5\) et peut être modifiée par la directive \mkey{Theta}. Les
  valeurs admissibles sont comprises entre \(0\) et \(1\)~;
  \item \varcpp{iterMax}, le nombre maximum d'itérations autorisées.
\end{itemize}

Si la méthode d'accélération de \castem{} est utilisée, les paramètres
suivants sont déclarés~:
\begin{itemize}
  \item \varcpp{accelerationTrigger}, le nombre d'itérations
  minimum pour utiliser la méthode d'accélération~;
  \item \varcpp{accelerationPeriod}, le nombre d'itérations
  entre deux utilisations de la méthode d'accélération~;
\end{itemize}

Si la méthode de relaxation de \castem{} est utilisée, les paramètres
suivants sont déclarés~:
\begin{itemize}
  \item \varcpp{relaxationTrigger}, le nombre d'itérations minimum
  pour utiliser la méthode de relaxation~;
  \item \varcpp{relaxationCoefficient}, le coefficient de
  relaxation.
\end{itemize}


\subsection{Intégration d'une loi d'écoulement viscoplastique
  orthotrope par un algorithme implicite}

Nous reprenons ici l'exemple précédent en modifiant l'algorithme
d'intégration utilisé par un algorithme implicite. Nous traitons ici la
généralisation de la loi de \nom{Norton} au cas orthotrope.

Les algorithmes implicites (\(\theta\)-méthodes) conduisent à réécrire
le système différentiel~\eqref{eq:law} ainsi~:
\begin{equation}
  \begin{aligned}
    \Delta\,\tepsilonel + \Delta\,p\,\tenseur{n} - \Delta\,\tepsilonto &= 0 \\
    \Delta{p}   - A\,\paren{\sigmaH}^{E}\,\Delta\, t &= 0 \\
  \end{aligned}
\end{equation}
où~:
\begin{minipage}[t]{0.8\linewidth}
  \begin{itemize}
  \item \(\Delta\,\tepsilonto\) est l'incrément sur le pas de temps de
    la déformation totale (imposée)~;
  \item \(\Delta\,\tepsilonel\) est l'incrément sur le pas de temps de
    la déformation élastique~;
  \item \(\Delta\,p\) est l'incrément sur le pas de temps de la
    déformation viscoplastique équivalente~;
  \item \(\tenseur{n}\) est la normale à la surface de charge évaluée
    à l'instant \(t+\theta\,\Delta\,t\)~;
  \item \(\sigmaH\) est la contrainte équivalente de \nom{Hill}
    évaluée à l'instant \(t+\theta\,\Delta\,t\)~;
  \item \(\Delta\,t\) est l'incrément de temps.
  \item \(\theta\) est un paramètre d'intégration compris entre \(0\)
    et \(1\).
  \end{itemize}
\end{minipage}

Ce système non-linéaire d'inconnues \(\Delta\,\tepsilonel\) et
\(\Delta\,p\) se décompose naturellement en deux équations
\(f_{\tepsilonel}\) et \(f_{p}\).

Ce système est résolu par une méthode de \nom{Newton-Raphson}, ce qui
nécessite de calculer sa matrice jacobienne. Cette matrice se décompose
en quatre parties notées respectivement
\(\deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel}\),
\(\deriv{f_{\tepsilonel}}{\Delta\,p}\),
\(\deriv{f_{p}}{\Delta\,\tepsilonel}\) et \(\deriv{f_{p}}{\Delta\,p}\).

Ces différents termes se calculent ainsi~:
\[
\left\{
\begin{aligned}
  \deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel} &= \tenseurq{I}+\Delta\,p\deriv{n}{\tsigma}\deriv{\tsigma}{\Delta\,\tepsilonel} \\
                                            &= \tenseurq{I} +\Delta\,p\deriv{n}{\tsigma}\deriv{\tsigma}{\mts{\tepsilonel}}\deriv{\mts{\tepsilonel}}{\Delta\,\tepsilonel} \\
                                            &= \tenseurq{I} +2\,\mu\,\theta\,\Delta\,p\,\deriv{n}{\tsigma} \\
  &= \tenseurq{I}+\Frac{2\,\mu\,\theta\,\Delta\,p}{\sigmaeq}\,\paren{\tenseurq{H}-\tenseur{n}\otimes\tenseur{n}} \\
  \deriv{f_{\tepsilonel}}{\Delta\,p} &= \tenseur{n} \\
  \deriv{f_{p}}{\Delta\,\tepsilonel} &= - A\,E\,\Delta\, t\,\paren{\sigmaH}^{E-1}\,\deriv{\sigmaH}{\Delta\,\tepsilonel} \\
  &= - A\,E\,\theta\,\Delta\, t\,\paren{\sigmaH}^{E-1}\,\tenseur{n}\,\colon\,\tenseurq{D} \\
  \deriv{f_{p}}{\Delta\,p} &= 1
\end{aligned}
\right.
\]
où \(\tenseur{H}\) est le tenseur de \nom{Hill} du matériau déjà
rencontré au paragraphe~\ref{sec:integration-dune-loi}.

\begin{figure}[htbp]
  \centering
  \begin{minipage}[htbp]{0.9\linewidth}
    \shorthandoff{:}
    \code{
      \small
      \input{@abs_top_srcdir@/docs/mfront/mfront/ImplicitOrthotropicCreep.tex}
    }
    \shorthandon{:}
  \end{minipage}
  \caption{Implantation d'une loi d'écoulement viscoplastique
    orthotrope par un algorithme implicite}
  \label{fig:OrthoViscoImplicit}
\end{figure}

\paragraph{Implantation} La figure~\ref{fig:OrthoViscoImplicit}
présente l'implantation de cette loi. La définition du système
implicite se fait dans la partie \mkey{Integrator}. Les variables {\tt
  feel}, {\tt fp}, {\tt dfeel\_ddeel}, {\tt dfeel\_ddp}, {\tt
  dfp\_ddeel}, {\tt dfp\_ddp} sont automatiquement définies. Les termes
du système sont initialisés à la valeur des incréments des variables
correspondantes ({\tt feel} est initialisé à {\tt deel}) et la matrice
jacobienne est initialisée à l'identité.

Cette implantation tient compte des limitations du code aux éléments
finis \castem{} qui ne permet pas d'avoir une définition homogène des
directions d'orthotropie quelle que soit l'hypothèse de
modélisation. Les conventions utilisées ici sont conformes à celles
décrites en annexe~\ref{sec:annexe:orthotropie} pour les tubes.

\clearpage
\newpage
\section{Conclusions}

Cette note a décrit les possibilités de \mfront{} dans le domaine de
l'intégration des lois de comportement mécanique. Elle complète la
présentation générale de \mfront{}~\cite{helfer_generateur_2013-1}.

Pour des raisons pédagogiques, nous avons choisi de ne présenter dans
cette note que des exemples élémentaires. Des exemples plus complexes
peuvent être trouvés dans les sources de la librairie \TFEL{}. Pour
les applications de la plate-forme \pleiades{}, des lois de
comportement des matériaux usuels du nucléaire peuvent être trouvés
dans la base de données \sirius{}. Les utilisateurs \cea{} de
l'application \licos{} peuvent également accéder à un corpus de
connaissances matériau précompilées via le projet
\mfm{}~\cite{helfer_connaissances_2013}.

\subsection{Une solution flexible, mature et performante}

\mfront{} constitue aujourd'hui une solution flexible et mature au
c\oe{}ur de la gestion des matériaux de la plate-forme \pleiades{}. Il
a été évalué en interne ou par les partenaires industriels du
\cea{}~\cite{olagnon_analysis_2013}.

\paragraph{L'évaluation par l'équipe \aster{}}
Pour illustrer ce point, nous pouvons nous appuyer sur l'évaluation
faite par l'équipe de code \aster{} qui est particulièrement
complète~\cite{proix_integration_2013}.

Différentes classes de lois de comportement ont été testées~:
\begin{itemize}
\item loi de \nom{Mazars} décrivant l'endommagement fragile du
  béton~\cite{hamon_modeendommagement_2013}~;
\item l'endommagement sous fluage tertiaire de structure en acier par
  une loi dite de \nom{Hayhurst}~\cite{proix_comportement_2012}~;
\item lois plastiques et viscoplastiques de type \nom{Chaboche} avec
  écrouissage cinématique (linéaire ou non linéaire)~;
\item comportement mono-cristallin et poly-cristallins pour des
  structures CFC.
\end{itemize}
De plus, l'équipe \aster{} a reversé de nombreux tests de loi de
comportement dans la librairie \TFEL{}. Ces tests ont été écrits à
l'aide de l'utilitaire \mtest{}~\cite{helfer_mtest_2013}.

\begin{table}[htbp]
  \centering
  \begin{tabular}[htbp]{|>{\centering}m{3.cm}|>{\centering}m{3.cm}|>{\centering}m{3.cm}|m{1mm}c|}
    \cline{2-5}
    \multicolumn{1}{c|}{}&
    \mfront{} implicite (\(19\) équations) &
    {\tt VISC\_CIN1\_CHAB} implicite aster optimisé \\(1 équation) & \vrule width 0mm height 1.4cm&
    \multirow{3}{*}{\includegraphics[width=6.cm]{@top_srcdir@/docs/mfront/Images/hsnv125d.eps}}\\
    \cline{1-3}
    Nombre de pas de temps                 &76          &76      & \vrule width 0mm height 1.4cm&\\
    \cline{1-3}
    Nombre d'itérations de \nom{Newton}    &160         &151     & \vrule width 0mm height 1.4cm&\\
    \cline{1-3}
    Temps CPU                              &4,67s      &4,19s    & \vrule width 0mm height 1.4cm&\\
    \hline
  \end{tabular}
  \caption{Résultat du test {\tt HSNV125}~\cite{proix_hsnv125_2011}.}
  \label{fig:mfront:behaviours:hsnv125}
\end{table}

\paragraph{Test de robustesse} Le test {\tt HSNV125} est extrait de la
base de cas test de \aster{}~\cite{proix_hsnv125_2011}. Il s'agit d'un
test particulièrement sollicitant pour les algorithmes d'intégration
des loi de comportement. Un exemple de résultat de ce test, qui permet
d'apprécier la complexité du chargement imposé, est donné au
tableau~\ref{fig:mfront:behaviours:hsnv125} pour une loi de type
\nom{Chaboche} avec deux variables cinématiques non linéaires. Dans
cet exemple, la loi \mfront{} est légèrement moins efficace que la loi
optimisée d'\aster{}\footnote{Une version non optimisée de la loi est
  également disponible dans \aster{} (\(21\) équations résolues par
  une méthode implicite). Cette version conduit à une divergence du
  calcul pour le test proposé.}.

\begin{table}[htbp]
  \centering
  \begin{tabular}[htbp]{|>{\centering}m{2.5cm}|>{\centering}m{2.5cm}|>{\centering}m{2.5cm}|>{\centering}m{2.5cm}|c|}
    \cline{2-5}
    \multicolumn{1}{c|}{}&
    MFRONT Implicite&
    MFRONT implicite avec jacobienne numérique&
    Aster avec jacobienne numérique &
    \multirow{6}{*}{\includegraphics[width=3.25cm]{@top_srcdir@/docs/mfront/Images/hayhurst.eps}}\\
    \cline{1-4}
    \cline{1-4}
    Nombre de pas de temps                 &11          &11            &11& \\
    \cline{1-4}
    Nombre d'itérations de \nom{Newton}    &37          &36            &44& \\
    \cline{1-4}
    Temps CPU                              &8mn11s      &7mn58s    &17mn43s& \\
    \cline{1-4}
    Temps d'intégration de la loi de comportement &2mn26s &2mn21s  &10mn51s& \\
    \cline{1-4}
    Temps d'inversion de la matrice               &5mn11s &5mn5s   &6mn11s& \\
    \hline
  \end{tabular}
  \caption{Résultat des test de fluage sur éprouvette entaillée $3D$,
    loi de \nom{Hayhurst}~\cite{proix_integration_2013}.}
  \label{tab:mfront:behaviours:hayhurst}
\end{table}

\paragraph{Tests de performances} Plusieurs comparaisons avec les lois
implantées dans \aster{} ont été menées. Dans l'ensemble, les lois
générées par \mfront{} se sont avérées efficaces, autant que les lois
\aster{} et parfois meilleures.

Un premier exemple est donnée par le
tableau~\ref{tab:mfront:behaviours:hayhurst} où l'on reporte les
résultats obtenus sur une éprouvette entaillée subissant un essai de
fluage. Le calcul est effectué en 3D en grandes déformations grâce aux
formalisme des déformations logarithmiques qui permet d'utiliser des
lois écrites en petites
déformations~\cite{bargellini_modelisation_2013,miehe_anisotropic_2002}.

\begin{table}[htbp]
  \centering
  \begin{tabular}[htbp]{|>{\centering}m{2.cm}|>{\centering}m{2.cm}|>{\centering}m{2.cm}|m{1mm}m{8.cm}|}
    \cline{2-5}
    \multicolumn{1}{c|}{}&
    \mfront{} (RK42)  &
    \aster{}  & \vrule width 0mm height 2.75cm&
    \multirow{2}{*}{\includegraphics[width=7.cm]{@top_srcdir@/docs/mfront/Images/polybz.eps}}\\
    \cline{1-3}
    Temps CPU                              &90s      &136s    & \vrule width 0mm height 2.75cm&\\
    \hline
  \end{tabular}
  \caption{Comparaison des temps de calculs pour une loi issue de l'homogénéisation d'un poly-cristal constitué de $100$ grains.}
  \label{fig:mfront:behaviours:polybz}
\end{table}

Un second test, basé sur une loi d'homogénéisation de lois de
plasticité cristalline, démontre les capacités de \mfront{} à traiter
un grand nombre de variables internes (plus de \(4\,200\) ici) en
conservant des performances numériques intéressantes, reportées au
tableau~\ref{fig:mfront:behaviours:polybz}. Le fichier \mfront{} fait
une centaine de lignes.

\subsection{Perspectives}

L'utilisation de \mfront{} hors de la plate-forme \pleiades{}, au DMN
et à EDF notamment, offre des possibilités de qualification de l'outil
extrêmement intéressantes et des voies de collaboration extrêmement
riches et enrichissantes.

\subsection{Développement}

Dans le cadre \pleiades{} différentes perspectives ont été
identifiées~:
\begin{itemize}
\item un meilleur support des contraintes planes~;
\item la possibilité d'écrire des lois de comportements poutre~;
\item l'amélioration de la robustesse de algorithmes implicites, qui
  peuvent diverger si l'estimation initiale est mauvaise. La méthode
  de \nom{Powell} semble particulièrement
  intéressante~\cite{chen_modification_1981}.
\end{itemize}

Pour une utilisation dépassant le cadre combustible, différents
développements sont envisagés. Ces développements se résument
généralement à la possibilité d'utiliser des lois de comportement qui
ne soient pas limitées au cas volumique petites déformations. Ceci
inclut~:
\begin{itemize}
\item le cas des poutres qui utilisent des variables généralisées (on
  retrouve le besoin combustible)~;
\item les grandes transformations~;
\item le cas des modèles de zones cohésives qui travaillent en
  traction/saut de déplacement~;
\item les lois à gradients (de déformations ou à variables internes).
\end{itemize}

À moyen terme, ces développements intéressent également nos
applications combustibles~:
\begin{itemize}
\item les grandes transformations semblent intéressantes dans des
  contextes accidentels (RIA, APRP)~;
\item les modèles de zones cohésives sont utilisés au DMN/SEMI pour la
  modélisation de la propagation d'une fissure par corrosion sous
  contrainte dans la gaine~;
\item les lois à gradients permettraient une meilleure représentation
  des phénomènes d'endommagement dans le cadre des travaux des
  thésards du SESC (modélisation de l'essai de flexion des barreaux
  \uod{}).
\end{itemize}

La structure de \mfront{} fait que tous ces développements sont
informatiquement connexes et peuvent être regroupés dans une même
action.

\clearpage
\newpage
\referencecea
\listetableaux
\listefigures

\appendix

\section{Résolution d'un problème
  mécanique quasi-statique par la méthode des éléments finis}
\label{sec:mfront:mechanical_equilibrium}

Nous décrivons succinctement la résolution par la méthode des éléments
finis utilisées pour décrire l'évolution quasi-statique d'un système
mécanique. Ceci nous permet de bien préciser la place de la loi de
comportement et de faire apparaître la notion de matrice tangente
cohérente.

La présentation faite est suffisamment générale pour notre propos mais
nous avons passé sous silence de nombreux points (transformations
finies, gestion du contact, etc...). Chaque code aux éléments finis
présentera sa propre variation, plus au moins
sophistiquée~\cite{besson_mecanique_2001,pascal_notice_2005,abbas_algorithme_2013}.

\subsection{Discrétisation}

La méthode des éléments finis repose sur une approximation de l'espace
des solutions par un espace de fonctions de dimension finie.

\paragraph{Maillage de la discrétisation} Un espace
d'approximation adapté à de tels espaces est construit en approchant le
volume $\Omega$ d'intérêt par un volume $\Omega^{h}$ appelé le {\em
  maillage} obtenu par agrégation d'un ensemble de polygones convexes
appelés {\em mailles} ou {\em éléments}. Ces polygones définissent~:
\begin{itemize}
  \item un ensemble de points particuliers incluant en
  particulier leur sommet où sont estimées les inconnues~;
  \item des fonctions d'interpolation qui permettent des
  construire les solutions approchées.
\end{itemize}

\paragraph{Notion de champ nodal} La discrétisation
éléments finis permet de définir la notion de {\em champ nodal}~: à
chaque n{\oe}ud est associé des valeurs qui permettent de construire une
approximation d'une fonction, dont la solution du problème mécanique, le
champ de déplacements.

\subsection{Principe des travaux virtuels}

La méthode présentée repose sur un découpage temporel. Connaissant le
champ de déplacement et l'état mécanique des matériaux, caractérisé
par des variables internes \(y_{i}\), à l'instant \(t\), la méthode
des éléments finis propose de rechercher sur une géométrie
discrétisée, un champ nodal \(\Delta\discret{\vec{u}}\) telle que la
fonction associée \(\Delta\vec{u}^{h}\) vérifie une approximation du
principe des travaux virtuels~:
\[
\forall \vec{v}^{h}\; C.A.^{h},\quad
\underbrace{\int_{\Omega^{h}}\tsigma_{t+\Delta
    t}\paren{\Delta\tepsilonto,t,\Delta\,
    t}\colon\tepsilonto\paren{\vec{v}^{h}}\;\dtot V}_{\text{Travaux
    intérieurs}}=\underbrace{\int_{\Omega^{h}}\vec{f}_{t+\Delta
    t}.\vec{v}^{h}\;\dtot
  V+\int_{\partial\Omega^{h}_{\vec{T}}}\vec{T}_{t+\Delta
    t}.\vec{v}^{h}\;\dtot S}_{\text{Travaux extérieurs}}
\]
où nous avons noté \(\Delta\tepsilonto\) la déformation totale
associée à la fonction \(\Delta\vec{u}^{h}\)
\(\paren{\Delta\tepsilonto=\tepsilonto\paren{\Delta\vec{u}^{h}}}\). La
dépendance au \(t\) symbolise la dépendance de la réponse du matériau à
des variables externes (fluence, flux neutronique, etc...). L'espace
\(C.A.^{h}\) est l'ensemble des fonctions \(\vec{v}^{h}\) correspondant
aux champs nodaux vérifiant les conditions aux limites cinématiques.
L'espace \(C.A.^{h}\) étant fini, il existe deux champs nodaux
\(\forceintEF\paren{\Delta\discret{\vec{u}}}\) et \(\forceextEF\) tels
que~:
\[
\begin{aligned}
  \forall \vec{v}^{h}\; C.A.^{h},\quad &
  \forceintEF\paren{\Delta\discret{\vec{u}}}.\champEF&=&
  \int_{\Omega^{h}}\tsigma_{t+\Delta t}\paren{\Delta\tepsilonto,\Delta\,
    t}\colon\tepsilonto\paren{\vec{v}}\;\dtot V &\\
  &\forceextEF.\champEF&=& \int_{\Omega^{h}}\vec{f}_{t+\Delta
    t}.\vec{v}^{h}\;\dtot
  V+\int_{\partial\Omega^{h}_{\vec{T}}}\vec{T}_{t+\Delta
    t}.\vec{v}^{h}\;\dtot S &
\end{aligned}
\]
Les deux champs nodaux \(\forceintEF\paren{\Delta\discret{\vec{u}}}\)
et \(\forceextEF\) sont appelés respectivement les {\em forces nodales
  intérieures} et {\em extérieures}. Le principe des travaux virtuels
est alors équivalent à l'égalité des forces nodales intérieures et
extérieures~:
\begin{equation}
  \label{eq:bibliographie:equilibre_forces_interieures_forces_exterieures}
 \boxed{\forceintEF\paren{\Delta\discret{\vec{u}}}=\forceextEF}
\end{equation}

\paragraph{Calcul des forces intérieures} Les forces
nodales intérieures peuvent être calculées à partir de la contribution
de chacun des éléments finis. La reconstruction du vecteur
\(\forceintEF\) à partir des contributions élémentaires
\(\forceintElem\) est l'étape d'{\em assemblage}. Sur chaque élément
fini, cette contribution élémentaire est calculée par intégration du
champ de contraintes~:
\[
\forceintElem=\int_{V^{e}}\tsigma_{t+\Delta
  t}\paren{\Delta\,\tepsilonto,\Delta\, t}\colon\tenseur{B}\;\dtot V
\]
où \(V^{e}\) est le volume de l'élément fini et \(\tenseur{B}\) une
matrice propre à l'élément fini permettant de définir la déformation de
la fonction interpolée à partir des valeurs du champ nodal aux n{\oe}uds
de l'élément. La matrice \(\tenseur{B}\) est directement reliée aux
gradients des fonctions d'interpolation.

\paragraph{Intégration numérique, points de
  \nom{Gauss}} Numériquement, l'intégrale
\(\displaystyle\int_{v^{e}}\tsigma_{t+\Delta t}\colon\tenseur{B}\;\dtot
V\) est évaluée par quadrature~: la fonction \(\tsigma_{t+\Delta
  t}\colon\tenseur{B}\) est évaluée en un certain nombre de points
choisis pour minimiser l'erreur commise. Ces points particuliers sont
les {\em points de \nom{Gauss}} de l'élément. La formule de quadrature
s'écrit alors~:
\begin{equation}
  \label{eq:bibliographie:quadrature_gauss}
  \forceintElem=\sum_{i=1}^{N^{G}} \paren{\tsigma_{t+\Delta
      t}\paren{\Delta\tepsilonto\paren{\vec{\eta}_{i}},\Delta\, t}\colon\tenseurq{B}\paren{\vec{\eta}_{i}}}w_{i}
\end{equation}
où \(N^{G}\) est le nombre de points de \nom{Gauss} de l'élément,
\(\vec{\eta}_{i}\) leurs coordonnées et \(w_{i}\) un poids associé au
\(i\)\textsuperscript{ième} point de \nom{Gauss}. Les conséquences de
cette formule sont importantes~: les contraintes et, par conséquent, les
variables d'états \(z_{i}\), n'ont besoin d'être connues qu'aux points
de \nom{Gauss}.

L'étape d'{\em intégration locale}, c'est à dire la détermination
de la contrainte en fin de pas de temps, est donc effectuée aux points
de \nom{Gauss}.

\paragraph{Notion de résidu} L'équation d'équilibre
discrétisée~\eqref{eq:bibliographie:equilibre_forces_interieures_forces_exterieures}
est écrite classiquement en introduisant le {\em résidu}
\(\residuEF\paren{\Delta\discret{\vec{u}}}\)~:
\begin{equation}
  \label{eq:bibliographie:equilibre_residu}
  \residuEF\paren{\Delta\discret{\vec{u}}}=\discret{\vec{O}}\quad\text{
    avec
  }\quad\residuEF\paren{\Delta\discret{\vec{u}}}=\forceintEF\paren{\Delta\discret{\vec{u}}}-\forceextEF
\end{equation}
La résolution de l'équilibre de la structure est équivalente à la
recherche d'un champ nodal annulant le résidu.

\subsection{Principe de la méthode de \nom{Newton-Raphson}}

La {\em méthode de \nom{Newton-Raphson}} propose de rechercher de
manière itérative la solution à l'équation
\eqref{eq:bibliographie:equilibre_residu}. Une estimation
\(\Delta\discret{\vec{u}}^{n+1}\) de cette solution est construite à
partir de l'estimation \(\Delta\discret{\vec{u}}^{n}\) à l'étape
précédente. La relation de récurrence entre
\(\Delta\discret{\vec{u}}^{n+1}\) et \(\Delta\discret{\vec{u}}^{n}\)
s'obtient en écrivant le développement limité à l'ordre \(1\) du résidu
\(\residuEF\paren{\Delta\discret{\vec{u}}^{n+1}}\) en supposant que les
estimations \(\Delta\discret{\vec{u}}^{n+1}\) et
\(\Delta\discret{\vec{u}}^{n}\) proches~:
\[
\residuEF\paren{\Delta\discret{\vec{u}}^{n+1}}\approx
\residuEF\paren{\Delta\discret{\vec{u}}^{n}}+\left.\deriv{\residuEF}{\Delta\discret{\vec{u}}}\right|_{\Delta\discret{\vec{u}}^{n}}.\paren{\Delta\discret{\vec{u}}^{n+1}-\Delta\discret{\vec{u}}^{n}}
\]
En écrivant que cette approximation est solution de l'équation
\eqref{eq:bibliographie:equilibre_residu} et en supposant
\(\left.\deriv{\residuEF}{\Delta\discret{\vec{u}}}\right|_{\Delta\discret{\vec{u}}^{n}}\)
inversible, nous obtenons la relation de récurrence
suivante~:
\begin{equation}
  \label{eq:newton_raphson_iteration}
  \Delta\discret{\vec{u}}^{n+1}=\Delta\discret{\vec{u}}^{n}-\paren{\left.\deriv{\residuEF}{\Delta\discret{\vec{u}}}\right|_{\Delta\discret{\vec{u}}^{n}}}^{-1}.\residuEF\paren{\Delta\discret{\vec{u}}^{n}}
\end{equation}
Nous sommes donc amenés à résoudre le problème
\eqref{eq:bibliographie:equilibre_residu} par une {\em succession de
  problèmes linéaires}.

\paragraph{Matrice de raideur} Le calcul de la dérivée
\(\deriv{\residuEF}{\Delta\discret{\vec{u}}}\) se fait en supposant les
forces extérieures constantes au cours de la résolution~:
\[
\left.\deriv{\residuEF}{\Delta\discret{\vec{u}}}\right|_{\Delta\discret{\vec{u}}^{n}}=\tenseurq{\mathbb{K}}\quad\text{
 avec }\quad
\tenseurq{\mathbb{K}}=\left.\deriv{\forceintEF}{\Delta\discret{\vec{u}}}\right|_{\Delta\discret{\vec{u}}^{n}}
\]
La matrice \(\tenseurq{\mathbb{K}}\) est appelée {\em matrice de
  raideur} de la structure.

\paragraph{Calcul de la matrice de raideur} La méthode
des éléments finis permet de calculer la matrice de raideur globale par
{\em assemblage} de {\em matrices de raideur élémentaires}
\(\tenseurq{\mathbb{K}}^{e}\). Pour chaque élément, l'expression de
cette matrice de raideur élémentaire est~:
\begin{equation}
  \label{eq:bibliographie:matrice_raideur_elementaire}
  \tenseurq{\mathbb{K}}^{e}=\sum_{i=1}^{N^{G}}
  \mbox{}^{t}\tenseurq{B}\paren{\vec{\eta}_{i}}\colon\deriv{\Delta\tsigma}{\Delta\tepsilonto}\paren{\vec{\eta}_{i}}\colon\tenseurq{B}\paren{\vec{\eta}_{i}}w_{i}
\end{equation}
où apparaît la {\em matrice tangente cohérente}
\(\deriv{\Delta\tsigma}{\Delta\tepsilonto}\).

\paragraph{Matrice de raideur tangente cohérente}
L'utilisation de la matrice tangente cohérente ou de la matrice sécante
est la solution la plus performante en théorie car elle conduit à une
convergence quadratique de la méthode de \nom{Newton-Raphson}. Le coût
de la réactualisation de la matrice de raideur du système à chaque
itération peut cependant être important.


\clearpage
\newpage
\include{glossary}

\clearpage
\newpage
\section{Principales opérations tensorielles}
\label{sec:oper-tens-dans}

Nous décrivons dans ce paragraphe les opérations tensorielles fournies
par \TFEL{} pour les tenseurs d'ordre \(2\) et d'ordre \(4\).

\subsection{Conventions de représentation des tenseurs}

Les tenseurs d'ordre \(2\) symétriques peuvent être représentés par
des vecteurs et \og~tenseurs d'ordre \(4\)~\fg{}, peuvent être
représentés par des matrices.

En \(3D\), les tenseurs d'ordre \(2\) symétriques ont \(6\)
composantes. Un tenseur \(\tenseur{a}\) est représenté par le vecteur
suivant~:
\[
\left(
\begin{array}{c}
  \mathcal{A}_{0} \\
  \mathcal{A}_{1} \\
  \mathcal{A}_{2} \\
  \mathcal{A}_{3} \\
  \mathcal{A}_{4} \\
  \mathcal{A}_{5} \\
\end{array}
\right)
=
\left(
\begin{array}{c}
  a_{xx} \\
  a_{yy} \\
  a_{zz} \\
  \sqrt{2}a_{xy} \\
  \sqrt{2}a_{xz} \\
  \sqrt{2}a_{yz} \\
\end{array}
\right)
\]

Le racine de \(2\) sur les termes extradiagonaux sont là pour assurer
que le produit doublement contracté de deux tenseurs d'ordre \(2\)
symétriques \(\tenseur{a}\) et \(\tenseur{b}\) est égal au produit
scalaire de leurs représentations vectorielles.

\subsection{Opérations sur les tenseurs d'ordre 2}

Pour les tenseurs d'ordre \(2\), \TFEL{} propose~:
\begin{itemize}
  \item la négation d'un tenseur~;
  \item la somme de deux tenseurs~;
  \item la différence de deux tenseurs~;
  \item la multiplication à droite ou à gauche d'un tenseur par
  un scalaire~;
  \item la division à droite d'un tenseur par un scalaire~;
  \item le produit tensoriel de deux tenseurs grâce à l'opérateur
  \texttt{\^}. Le résultat est un tenseur d'ordre \(4\) \og~tenseurs
  d'ordre \(4\)~\fg{} \(\tenseur{a}\otimes\tenseur{b}\) dont la
  représentation matricielle vérifie a pour élément
  \[
  \tenseurq{c}=\tenseur{a}\otimes\tenseur{b}
  \quad\Rightarrow\quad\tenseurq{C}_{ij}=\mathcal{A}_{i}\,\mathcal{B}_{j}
  \]
  \item le produit contracté de deux tenseurs grâce à l'opérateur
  \texttt{|}. Comme indiqué plus haut, ce produit contracté est égal au
  produit scalaire de leurs représentations tensorielles~:
  \[
  \tenseur{a}\,\colon\,\tenseur{b} = \sum\,\mathcal{A}_{i}\,\mathcal{B}_{i}
  \]
\end{itemize}

En \cpp{}, ces opérations respectent l'ordre de priorité usuel sauf
pour les deux derniers qui peuvent nécessiter l'emploi de parenthèses.

\subsection{Opérations sur les tenseurs d'ordre 4}

Pour les tenseurs d'ordre \(4\), \TFEL{} propose~:
\begin{itemize}
  \item la négation d'un tenseur~;
  \item la somme de deux tenseurs~;
  \item la différence de deux tenseurs~;
  \item la multiplication à droite ou à gauche d'un tenseur par
  un scalaire~;
  \item la division à droite d'un tenseur par un scalaire~;
  \item la multiplication de deux tenseurs~;
  \item l'application d'un tenseur d'ordre \(4\) à un tenseur
  d'ordre \(2\) par l'opérateur de multiplication \texttt{*}. La
  représentation vectorielle du résultat de cette opération est égal à
  la multiplication vectorielle de la représentation matricielle du
  tenseur d'ordre \(4\) par la représentation vectorielle du tenseur
  d'ordre \(2\)~:
  \[
  \tenseur{c}=\tenseurq{b}\,\colon\,\tenseur{a} \quad \Rightarrow \quad \mathcal{C}_{i} = \sum_{j} \mathcal{B}_{ij}A_{j}
  \]
  \item le produit contracté à gauche d'un tenseur d'ordre \(4\)
  par un tenseur d'ordre \(2\) grâce à l'opérateur \texttt{|}~:
  \[
  \tenseur{c}=\tenseur{a}\,\colon\,\tenseurq{B} \quad \Rightarrow \quad \mathcal{C}_{i} = \sum_{j} \mathcal{A}_{j}\mathcal{B}_{ji}
  \]
\end{itemize}

Ces opérations respectent l'ordre de priorité usuel sauf pour le
dernier qui peuvent nécessiter l'emploi de parenthèses.

\clearpage
\newpage
\section{Fonctions utiles}
\label{sec:mfront:fcts}

Cette section décrit quelques fonctions utiles à la création de lois de
comportement mécanique fournies par la librairie \TFEL{}.

\subsection{Manipulation des tenseurs d'ordre $2$ symétriques}

Le fichier d'entête
\headerb{TFEL/Math/stensor.hxx}{TFEL/\-Math/\-stensor.hxx}, inclut par
toutes les lois de comportement mécanique, définit la classe
\tfelb{stensor<N,T,Storage>}{stensor} représentant les tenseurs d'ordre
$2$ symétriques.

\paragraph{La méthode \texttt{computeEigenValues}} La méthode
\tfelb{stensor<N,T,Storage>::computeEigenValues}{compute\-Eigen\-Values}
calcule les valeurs propres du tenseur d'ordre $2$ symétriques. Ses
valeurs propres doivent être passées en référence.

\paragraph{La méthode \texttt{computeEigenVectors}} La méthode
\tfelb{stensor<N,T,Storage>::computeEigenVectors}{compute\-Eigen\-Vectors}
calcule les valeurs et vecteurs propres du tenseur d'ordre $2$
symétriques. Son premier argument doit être un objet de type
\tfel{tvector}\footnote{Les objets de type \tfel{tvector} sont définis
  dans le fichier d'entête
  \headerb{TFEL/Math/tvector.hxx}{TFEL/\-Math/\-tvector.hxx}.} et
contiendra les différentes valeurs propres. Son second argument doit
être un objet de type \tfel{tmatrix} contenant les vecteurs
propres\footnote{Les objets de type \tfel{tmatrix} sont définis dans le
  fichier d'entête
  \headerb{TFEL/Math/tmatrix.hxx}{TFEL/\-Math/\-tmatrix.hxx}.}. Cet
objet peut être passé directement à la méthode
\tfelb{stensor<N,T,Storage>::changeBasis}{change\-Basis} .

\paragraph{La méthode \texttt{changeBasis}} La méthode
\tfelb{stensor<N,T,Storage>::changeBasis}{change\-Basis} effectue le
changement de repère d'un tenseur d'ordre \(2\) symétrique. Ce
changement de repère est défini par un objet de type \tfel{tmatrix}.

\paragraph{La fonction \texttt{trace}} La fonction \tfel{trace} renvoie
la trace d'un tenseur d'ordre $2$ symétrique.

\paragraph{La fonction \texttt{sigmaeq}} La fonction \tfel{sigmaeq}
renvoie la contrainte équivalente au sens de \nom{Von Mises} d'un
tenseur d'ordre $2$ symétrique.

\paragraph{La fonction \texttt{abs}} La fonction \tfel{abs} renvoie la
somme des valeurs absolues des composantes d'un tenseur d'ordre $2$
symétrique, les termes extradiagonaux sont comptés une fois et sont
affectés d'un facteur \(\sqrt{2}\).

\subsection{Calcul des coefficents de \nom{Lamé}}

Le fichier d'entête
\headerb{TFEL/MaterialLaw/Lame.hxx}{TFEL/\-Material\-Law/\-Lame.hxx},
qu'il est nécessaire d'inclure au niveau du fichier \mfront{} en
utilisant la directive \mkey{Includes}, définit différentes fonctions
utilitaires~:
\begin{itemize}
  \item \tfel{compute\-Lambda} permet de calculer le premier
  coefficent de \nom{Lamé}. Elle prend en argument le module
  d'\nom{Young} et le coefficient de \nom{Poisson}~;
  \item \tfel{compute\-Mu} permet de calculer le second coefficent de
  \nom{Lamé}, c'est à dire le module de cisaillement. Elle prend en
  argument le module d'\nom{Young} et le coefficient de \nom{Poisson}
\end{itemize}

Ce fichier d'entête propose également une classe
\tfelb{computeElasticStiffness}{com\-pute\-Elastic\-Stiff\-ness},
paramétré par la dimension d'espace et le type numérique utilisé qui
permet de calculer la matrice d'élasticité à partir des coefficients de
\nom{Lamé}. Cette classe propose en méthode statique nommée
\tfelb{computeElasticStiffness::exe}{exe} qui prend en arguments~:
\begin{itemize}
\item une référence un objet de type \tfelb{st2tost2<N,T>}{st2tost2\-<N,T>}
qui représente la matrice d'élasticité.
\item le premier coefficient de \nom{Lamé}~;
\item le second  coefficient de \nom{Lamé}.
\end{itemize}

\subsection{Calcul du tenseur de \nom{Hill}}
\label{sec:calcul-du-tenseur}

La librairie {\tt TFEL} fournit la fonction \tfel{hillTensor}. Cette
fonction est déclarée dans le fichier d'entête
\headerb{TFEL/MaterialLaw/Hill.hxx}{TFEL/\-Material\-Law/\-Hill.hxx}
qu'il est nécessaire d'inclure au niveau du fichier \mfront{} en
utilisant la directive \mkey{Includes}.

Cette fonction est paramétrée par la dimension d'espace, notée
\texttt{N} dans le fichier \mfront{}, et par le type numérique utilisé,
noté T. Ce type est appelé {\tt real} dans le fichier \mfront{}. Elle
prend en argument les \(6\) coefficients \(F\), \(G\), \(H\), \(L\),
\(M\) et \(N\) présentés aux paragraphes \ref{sec:critere-de-nomhill}.

Elle retourne un objet de type \tfelb{st2tost2<N,T>}{st2tost2} qui
représente une forme linéaire sur les tenseurs d'ordre \(2\) symétriques
et qui se manipule comme une matrice\footnote{Les objets de type
  \tfelb{st2tost2<N,T>}{st2tost2} sont définis dans le fichier d'entête
  \headerb{TFEL/Math/st2tost2.hxx}{TFEL/\-Math/\-st2tost2.hxx}.}.

\paragraph{Utilisation dans le cas des tubes} Avec les conventions
décrites au paragraphe~\ref{sec:conv-de-rang}
et~\ref{sec:trait-des-probl} et les identifications faites aux
paragraphes~\ref{sec:critere-de-nomhill} et~\ref{sec:trait-des-probl},
cette fonction s'utilise ainsi~:
\begin{center}
  \begin{minipage}[htbp]{0.9\linewidth}
    \shorthandoff{:} \code{ \small \input{@abs_top_srcdir@/docs/mfront/mfront/hillTensor.tex} }
    \shorthandon{:}
  \end{minipage}
\end{center}

% \clearpage
% \newpage
% \section{Amélioration de la robustesse des algorithmes de
%   résolutions implicites} \label{sec:mfront:implicit:robustness}

% Les méthodes de résolution des systèmes implicites nécessitent de
% choisir une estimation de la solution \(Y_{0}\) pour initier la
% recherche. Par défaut, \mfront{} utilisera le vecteur nul comme
% estimation initiale, mais l'utilisateur peut en spécifier une autre par
% la directive \mkey{Predictor}.

% La convergence des algorithmes de résolutions implicites n'est garantie
% que si l'estimation initiale \(Y_{0}\) est suffisamment proche de
% la solution recherchée. Loin de la solution, l'algorithme peut
% éventuellement converger mais les premières itérations tendent
% généralement vers la solution assez lentement (phase de recherche)~: il
% est classique d'observer une convergence linéaire lors de la phase de
% recherche et une transition vers un ordre de convergence plus élevé
% proche de la solution.

% Pour les lois de comportement fortement non linéaires, il est difficile
% de proposer une estimation \(Y_{0}\) garantissant la convergence,
% il peut donc être intéressant d'améliorer la robustesse des algorithmes.

% \begin{table}
%   \centering
%   \begin{tabular}{|c|c|}
%     \hline Méthode & Expression de \(s_{n}\) \\
%     \hline \hline Algorithme de \textsc{Newton-Raphson} &
%     \(-J^{-1}\,.\,F\paren{Y_{n-1}}\) \\
%     \hline
%     \begin{minipage}{7cm}
%       Algorithme de \textsc{Newton-Raphson} et
%       dérivation numérique
%     \end{minipage}
%     & \(-\underset{\sim}{J}^{-1}\,.\,F\paren{Y_{n-1}}\) \\
%     \hline Premier algorithme de \textsc{Broyden} &
%     \begin{minipage}{5cm}
%       \(-\underset{\sim}{J}_{n-1}^{-1}\,.\,F\paren{Y_{n-1}}\)
%       où \(\underset{\sim}{J}_{n}\) est donné par
%       l'équation~\eqref{eq:Broyden}
%     \end{minipage}
%     \\
%     \hline Second algorithme de \textsc{Broyden} &
%     \begin{minipage}{5cm}
%       \(-\underset{\sim}{J}_{n-1}^{-1}\,.\,F\paren{Y_{n-1}}\)
%       où \(\underset{\sim}{J}_{n}^{-1}\) est donné par\(Y_{0}\)
%       l'équation~\eqref{eq:Broyden2}
%     \end{minipage}
%     \\
%     \hline
%   \end{tabular}
%   \caption{Expression de $s_{n}$ pour les
%     différents algorithmes de résolution de système d'équations non
%     linéaires} \label{tab:sn}
% \end{table}

% L'ensemble des algorithmes présentés jusqu'ici prennent la forme~:
% \[
% Y_{n}=Y_{n-1}+s^{N}_{n}
% \]
% L'expression du pas \(s^{N}_{n}\) pour les différentes algorithmes
% présentés est donnée au tableau~\ref{tab:sn}. D'après la discussion
% précédente, ce pas \(s^{N}_{n}\) est optimal proche de la solution
% finale.

% Il est intéressant de trouver une autre estimation de la solution, moins
% performante près de la solution mais plus robuste. Les méthodes
% précédentes donnant une approximation du jacobien, il peut être
% intéressant de substituer à la fonction \(F\) son application tangente
% et de choisir cette nouvelle estimation comme la solution du problème de
% minimisation suivant~:
% \[
% \min_{s^{G}\,\in\,\mathcal{R}^{N}}\,\norm{F\paren{Y_{n-1}}+J\paren{Y_{n-1}}\,.\,s^{G}}
% \]

% La solution de ce problème est donnée par\footnote{Pour le prouver,
%   nous pouvons écrire que~:
%   \[
%   \begin{array}{rcl}
%     \norm{F\paren{Y_{n-1}}+J\paren{Y_{n-1}}\,.\,s^{G}}^{2}
%     &=& {F}\paren{Y_{n-1}}\,.\,{F}\paren{Y_{n-1}}+2\,F\paren{Y_{n-1}}\,\mid\,J\paren{Y_{n-1}}\,.\,s^{G} \\
%     &&+\paren{J\paren{Y_{n-1}}\,.\,s^{G}}\,\mid\,\paren{J\paren{Y_{n-1}}\,.\,s^{G}}
%   \end{array}
%   \]
%   Le gradient de cette expression est égal à~:
%   \[
%   2\,\mbox{}^{t}J\paren{Y_{n-1}}F\paren{Y_{n-1}}+\mbox{}^{t}J\paren{Y_{n-1}}\,J\paren{Y_{n-1}}\,.s^{G}
%   \]
% }~:
% \[
% \mbox{}^{t}J\paren{Y_{n-1}}\,J\paren{Y_{n-1}}\,.s^{G}=-\,\mbox{}^{t}J\paren{Y_{n-1}}\,.\,F\paren{Y_{n-1}}
% \]
% Ce pas est intéressant loin de la solution.

% Les pas \(s^{N}_{n}\) et \(s_{n}^{G}\) sont combinés ainsi~:
% \[
% s = a_{n} \,s^{G}_{n} + \paren{1-a_{n}}\,s_{n}^{N}
% \]

\clearpage
\newpage
\section{Description de l'algorithme \nom{Runge-Kutta}-\castem{}}
\label{sec:algoRKcastem}
% $Y$ représente les vecteurs de variables internes et

G représente une fonction a priori non linéaire et que nous supposerons
a minima continûment dérivable. On cherche à résoudre le système
différentiel suivant:
\begin{equation}
  \label{eq:systeme_diff:RK}
  \dot{Y}=G\paren{Y,t}
\end{equation}

La résolution numérique de cette équation consiste à estimer, à partir
de la valeur \(\debutpas{Y}\) à un instant \(t\), la valeur
\(\finpas{Y}\) à un instant \(t+\Delta\, t\).

Voici une description de l'algorithme \nom{Runge-Kutta}-\castem{}~:
\begin{itemize}
\item première estimation d'une solution à \(\Delta\, t /2\)~:
\begin{subequations}
\begin{align}
  \dot{Y}_1 &= G\paren{\debutpas{Y},t} \\
  \demipas{Y_1} &= \debutpas{Y} + \Frac{\Delta\, t}{2} \dot{Y}_1
\end{align}
\end{subequations}
\item deuxième estimation d'une solution à \(\Delta\, t /2\)~:
\begin{subequations}
\begin{align}
  \dot{Y}_2 &= G\paren{Y_1,t} \\
  \demipas{Y_{12}} &= \debutpas{Y} + \Frac{\Delta\, t}{2} \paren{\Frac{\dot{Y}_1 + \dot{Y}_2}{2}}
\end{align}
\end{subequations}
\item première estimation d'une solution à \(\Delta\, t\)~:
\begin{subequations}
\begin{align}
  \dot{Y}_3 &= G\paren{Y_{12},t} \\
  \finpas{Y_{13}} &= \demipas{Y_{12}} + \Frac{\Delta\, t}{2} \dot{Y}_3
\end{align}
\end{subequations}
\item deuxième estimation d'une solution à \(\Delta\, t\)~:
\begin{subequations}
\begin{align}
  \dot{Y}_4 &= G\paren{Y_{13},t} \\
\label{eq:estimation_yf}
\finpas{Y_f} &= \demipas{Y_{12}} + \Frac{\Delta\,
  t}{2} \paren{\Frac{\dot{Y}_3 + \dot{Y}_4}{2}}
\end{align}
\end{subequations}
\item troisième estimation d'une solution à \(\Delta\, t\)~:
\begin{subequations}
\begin{align}
  \dot{Y}_5 &= G\paren{Y_f,t} \\
\label{eq:estimation_y5}
\finpas{Y_5} &= \debutpas{Y} + \Delta\, t \paren{\Frac{\dot{Y}_1 +
    4\dot{Y}_3 + \dot{Y}_5}{6}}
\end{align}
\end{subequations}
\end{itemize}

Dans le cas où le calcul d'une de ces estimations ne pourraient être
calculés, un sous-découpage du pas de temps est effectué comme le
montre l'algorithme décrivant le sous-découpage du pas temps. Ceci est
illustré en figure~\ref{fig:algo_temps}.

Une fois toutes ces estimations calculées, on compare les contraintes
obtenues lors des estimations \(\finpas{Y_f}\)
(\eqref{eq:estimation_yf}) et \(\finpas{Y_5}\)
(\eqref{eq:estimation_y5}), notées respectivement \(\finpas{\tsigma}\)
et \(\tsigma_5\).

Pour cela, on se donne une erreur dont l'amplitude est basée sur la
valeur des contraintes en début de pas \(\debutpas{\tsigma}\) ou sur la
valeur du module d'\nom{Young} \(E\)~:
\[
errabs = \left\{
\begin{aligned}
  \sqrt{\debutpas{\tsigma}:\debutpas{\tsigma}}\,\times\,10^{-5} &\textrm{ si }\sqrt{\debutpas{\tsigma}:\debutpas{\tsigma}}>E\,\times\,10^{-3} \\
  E\,\times\,10^{-8}                       &\textrm{ si }\sqrt{\debutpas{\tsigma}:\debutpas{\tsigma}}<E\,\times\,10^{-3}
\end{aligned}
\right.
\]

On définit deux variables représentant
l'écart des contraintes \(ra\) et \(sqra\) définies de la manière
suivante :
\begin{align}
  ra &=
  \Frac{\sqrt{\paren{\finpas{\tsigma}-\tsigma_5}:\paren{\finpas{\tsigma}-\tsigma_5}}}{errabs}\\
  sqra &= \sqrt{ra}\\
\end{align}

En testant ces deux paramètres avec des critères (fixés), on établit la
convergence ou la non-convergence du calcul ainsi que le nouvel
incrément de temps \(\Delta\, t\) qui sera utilisé. Dans le code de
calcul, cet incrément de temps \(\Delta\, t\) est noté \texttt{dt\_}.
L'algorithme déterminant le nouvel incrément de temps est présenté à la
figure~\ref{fig:algo_temps}. Dans cet algorithme, le temps est initial
est noté \texttt{t} tandis que le temps final est noté \texttt{dt}. La
variable \texttt{dtprec} est définie comme une valeur minimale du
sous-découpage \texttt{dt\_} en deçà de laquelle on renvoie une erreur.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\textwidth,height=22cm]{@top_srcdir@/docs/mfront/Images/Algo_temps.eps}
\caption{Description de l'algorithme pour déterminer le sous-découpage du pas de temps}
\label{fig:algo_temps}
\end{figure}

\clearpage
\newpage
\include{annexe-orthotropie}

\clearpage
\newpage
\section{Comparaisons numériques de différents
  algorithmes de résolution}
\label{QNR:ComparaisonNumeriques}

\subsection{Loi utilisée pour la comparaison}

Pour ces premières évaluations, nous avons utilisés un écoulement
viscoplastique de \nom{Norton}. Il s'agit d'une loi extrêmement simple
convergent en peu d'itérations. Le coût intrinsèque des différentes
méthodes est donc exacerbé. Cette loi présente deux variables internes
ont été retenues, la déformation élastique \(\tepsilonel\) et la
déformation viscoplastique cumulée \(p\).

\subsection{Résultats obtenus}

Nous comparons, dans le cas d'un essai en traction uniaxiale, l'
algorithme de \nom{Newton} à l'algorithme de \nom{Broyden}.

Pour l'algorithme de \nom{Newton}, nous avons introduit une variante
dans laquelle le jacobien est estimé numériquement (avec une différence
finie d'ordre \(1\) ou d'ordre \(2\)) et réactualisé à différentes
périodes (tableaux~\ref{tab:QNR:1} et~\ref{tab:QNR:2}).

Pour le premier algorithme de \nom{Broyden}, nous avons testé
différentes variantes suivantes (tableaux~\ref{tab:Broyden:1}
et~\ref{tab:Broyden:2})~:
\begin{itemize}
  \item que le jacobien initial était pris égal à l'identité ou
  évalué numériquement par une différence finie d'ordre \(2\)~;
  \item que des termes du jacobien étaient calculés
  explicitement.
\end{itemize}

Pour le second algorithme de \nom{Broyden}, nous avons distingué
deux variantes suivant que le jacobien inital était pris égal à
l'identité ou évalué numériquement par une différence finie d'ordre
\(2\).

\begin{figure}[htbp]
  \centering
  \includegraphics[height=0.8\linewidth,angle=-90]{@top_srcdir@/docs/mfront/Images/CompSeq.eps}
  \caption{Contraintes équivalentes obtenues respectivement par une
    résolution par un algorithme de \textsc{Newton-Raphson} ou par un
    algorithme de \textsc{Broyden}.}
  \label{fig:CompVMis}
\end{figure}

Quand ils convergent, les algorithmes donnent toujours la même
solution, ce qui est illustré en figure~\ref{fig:CompVMis}. Les coûts
des différentes méthodes, en termes de nombre cycles
processeurs\footnote{Il s'agit d'une mesure fiable du coût d'une
  fonction, indépendante de la charge du système. Nous avons utilisé
  l'outil \valgrind{} pour avoir accès à cette mesure.}, sont fournis
dans les tableaux~\ref{tab:NR}, \ref{tab:QNR:1},
\ref{tab:QNR:2}, \ref{tab:Broyden:1}, \ref{tab:Broyden:2}
et~\ref{tab:Broyden:3}.

Le test effectué étant extrêmement frustre, il serait mal avisé d'être
conclusif.

\begin{table}
  \centering
  \begin{tabular}[htbp]{|c|c|c|}
    \hline
    Variante & Nombre de cycles &
    \begin{minipage}{4cm}
      \begin{center}
        Ratio par rapport à \\
        l'algorithme de \textsc{Newton}
      \end{center}
    \end{minipage} \\
    \hline
    \hline
    \(J\) exact & \(4\,707\,335\)  & 1\\
    \hline
    \begin{minipage}[p]{5cm}
      \begin{center}
        \(J\) égal à l'identité
      \end{center}
    \end{minipage}
    & pas de convergence  & \\
    \hline
    \begin{minipage}[p]{5cm}
      \begin{center}
        \(\deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel}\) et
        \(\deriv{f_{p}}{\Delta\,p}\) exacts
      \end{center}
    \end{minipage} &
    pas de convergence & \\
    \hline
    \begin{minipage}[p]{5cm}
      \begin{center}
        \(\deriv{f_{p}}{\Delta\,\tepsilonel}\) et
        \(\deriv{f_{p}}{\Delta\,p}\) exacts, et
        \(\deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel}\) égal à
        l'identité
      \end{center}
    \end{minipage} &
    pas de convergence & \\
    \hline
    \begin{minipage}[p]{5cm}
      \begin{center}
        \(\deriv{f_{\tepsilonel}}{\Delta\,p}\) et
        \(\deriv{f_{p}}{\Delta\,p}\) exacts, et
        \(\deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel}\) égal à
        l'identité
      \end{center}
    \end{minipage} &
    \(152\,721\,841\) & \(32.44\)\\
    \hline
  \end{tabular}
  \caption{Temps calculs obtenus avec l'algorithme de
    \textsc{Newton-Raphson}.}
  \label{tab:NR}
\end{table}

Le tableau~\ref{tab:NR} montre à quel point l'algorithme de \nom{Newton}
est sensible à la qualité de la matrice jacobienne.

\begin{table}
  \centering
  \begin{tabular}[htbp]{|c|c|c|}
    \hline
    Variante & Nombre de cycles &
    \begin{minipage}{4cm}
      \begin{center}
        Ratio par rapport à \\
        l'algorithme de \textsc{Newton}
      \end{center}
    \end{minipage} \\
    \hline
    \hline
    \(\underset{\sim}{J}\) réactualisé à chaque itération          & \(9\,900\,064\) & \(2,10\) \\
    \hline
    \(\underset{\sim}{J}\) réactualisé toutes les \(2\) itérations & \(8\,674\,464\) & \(1,84\) \\
    \hline
    \(\underset{\sim}{J}\) réactualisé toutes les \(3\) itérations & \(9\,159\,968\) & \(1,94\) \\
    \hline
  \end{tabular}
  \caption{Temps calculs obtenus avec l'algorithme de
    quasi-\textsc{Newton} et un calcul numérique du jacobien par une
    différence finie (approximation d'ordre $1$).}
  \label{tab:QNR:1}
\end{table}

\begin{table}
  \centering
  \begin{tabular}[htbp]{|c|c|c|}
    \hline
    Variante & Nombre de cycles &
    \begin{minipage}{4cm}
      \begin{center}
        Ratio par rapport à \\
        l'algorithme de \textsc{Newton}
      \end{center}
    \end{minipage} \\
    \hline
    \hline
    \(\underset{\sim}{J}\) réactualisé à chaque itération & \(14\,916\,316\)  & \(3,16\)        \\
    \hline
    \(\underset{\sim}{J}\) réactualisé toutes les \(2\) itérations & \(11\,773\,880\) & \(2,5\) \\
    \hline
    \(\underset{\sim}{J}\) réactualisé toutes les \(3\) itérations & \(13\,330\,216\) & \(2,83\) \\
    \hline
  \end{tabular}
  \caption{Temps calculs obtenus avec l'algorithme de
    quasi-\textsc{Newton} et un calcul numérique du jacobien par une
    différence finie centrée (approximation d'ordre $2$).}
  \label{tab:QNR:2}
\end{table}

Les tableaux~\ref{tab:QNR:1} et~\ref{tab:QNR:2} montrent que les
évalutions numériques des jacobiens sont assez coûteuses. On peut
également noté que si l'on gagne du temps à n'évaluer la matrice que
toutes les deux itérations, on en perd à ne l'évaluer toutes les trois
itérations car la convergence est dégradée. Il n'y a pas de gain à
utiliser une différence finie d'ordre \(2\).

Il est intéressant de noter que les algorithmes de \nom{Broyden}
convergent toujours.

Le cas où tous les termes de la jacobienne sont fournis permet de
comparer le coût intrinsèque plus élevé du premier algorithme de
\nom{Broyden} (tableau~\ref{tab:Broyden:2}).

Si aucun terme du jacobien n'est fourni, le premier algorithme de
\nom{Broyden} peut être relativement efficace, à condition d'estimer
numériquement la jacobienne initiale (tableau~\ref{tab:Broyden:2}).

De plus, il est intéressant de s'intéresser au terme
\(\deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel}\) de la matrice
jacobienne. Il s'agit du terme le plus complexe (et le plus coûteux) à
calculer. Nous pouvons noter que si l'on ne calcule pas ce terme dans
l'algorithme de \nom{Broyden}, mais que l'on calcule les autres,
l'algorithme converge et est légèrement meilleur que l'algorithme de
\nom{Newton}.

\begin{table}
  \centering
  \begin{tabular}[htbp]{|c|c|c|}
    \hline
    Variante & Nombre de cycles &
    \begin{minipage}{5cm}
      \begin{center}
        Ratio par rapport à \\
        l'algorithme de \textsc{Newton}
      \end{center}
    \end{minipage} \\
    \hline
    \hline
    Défaut    & \(20\,197\,423\) & \(4,29\)\\
    \hline
    \(J\) exact & \(6\,120\,862\) & \(1,3\)\\
    \hline
    \(\deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel}\) exact &
    \(36\,821\,766\) & \(7,82\)\\
    \hline
    \(\deriv{f_{p}}{\Delta\,p}\) exact &
    \(33\,826\,368\) & \(7,186\)\\
    \hline
    \(\deriv{f_{p}}{\Delta\,\tepsilonel}\) exact &
    \(19\,577\,698\) & \(4,15\)\\
    \hline
    \(\deriv{f_{\tepsilonel}}{\Delta\,p}\) exact &
    \(12\,132\,956\) & \(2,58\)\\
    \hline
    \(\deriv{f_{\tepsilonel}}{\Delta\,p}\) et \(\deriv{f_{p}}{\Delta\,\tepsilonel}\) exacts &
    \(4\,686\,228\) & \(0,995\) \\
    \hline
  \end{tabular}
  \caption{Temps calculs obtenus avec le premier algorithme de
    \textsc{Broyden} et jacobien initial égal à l'identité.}
  \label{tab:Broyden:1}
\end{table}

\begin{table}
  \centering
  \begin{tabular}[htbp]{|c|c|c|}
    \hline
    Variante & Nombre de cycles &
    \begin{minipage}{5cm}
      \begin{center}
        Ratio par rapport à \\
        l'algorithme de \textsc{Newton}
      \end{center}
    \end{minipage} \\
    \hline
    \hline
    Défaut & \(9\,535\,347\) & \(2,02\) \\
    \hline
    \(\deriv{f_{\tepsilonel}}{\Delta\,p}\) et \(\deriv{f_{p}}{\Delta\,\tepsilonel}\) exacts &
    \(8\,178\,547\) & \(1,73\) \\
    \hline
  \end{tabular}
  \caption{Temps calculs obtenus avec le premier algorithme de
    \textsc{Broyden} et jacobien initial approximé numériquement.}
  \label{tab:Broyden:2}
\end{table}

\begin{table}
  \centering
  \begin{tabular}[htbp]{|c|c|c|}
    \hline
    Variante & Nombre de cycles &
    \begin{minipage}{5cm}
      \begin{center}
        Ratio par rapport à \\
        l'algorithme de \textsc{Newton}
      \end{center}
    \end{minipage} \\
    \hline
    \hline
    Défaut    & \(19\,530\,077\) & \(4,15\)\\
    \hline
  \end{tabular}
  \caption{Temps calculs obtenus avec le second algorithme de
    \textsc{Broyden} et jacobien initial égal à l'identité.}
  \label{tab:Broyden:3}
\end{table}

\clearpage
\newpage
\include{InternalNames}

%\clearpage
%\newpage
%\printindex{env}{Index des variables d'environnement}

\clearpage
\newpage
\printindex{bheaders}{Index des fichiers d'entête fournis par la
  librairie \TFEL{}}

\clearpage
\newpage
\printindex{btfel}{Index des classes, des méthodes et des fonctions
  fournies par la librairie \TFEL{}}

\clearpage
\newpage
\printindex{bmkeys}{Index des directives}

\end{document}