This page describes a general framework to implement finite strain single crystal behaviors. The reader is strongly encouraged to read the following page on how to describe of the slip systems, orientation tensors and of the interaction matrices of a single crystal.

The first section is dedicated to the description of the thermodynamic bases of the framework presented, based on the work of Xhu Han (see [1]).

The second part describes the implicit scheme used to integrate this behaviour and the computation of the consistent tangent operator. This part is based on the Code_Aster documentation ([2]).

The third part details the MFront implementation of this behaviour.

The last part introduces the FiniteStrainSingleCrystal brick which hides all the details described in the previous sections.

A generic framework for finite strain single crystal plasticity

Thermodynamical basis

Les notions principales liées au formalisme en grandes déformations sont présentées. Cette partie s’appuie sur les références et .

Le vecteur \(\vec{X}\) désigne les coordonnées d’un point dans la configuration initiale \(\mathit{C}_0\) et le vecteur \(\vec{x}\) désigne les coordonnées dans la configuration finale \(\mathit{C}\). Le gradient de la transformation s’écrit : \begin{equation} = ]

Pour la plasticité cristalline, la décomposition multiplicative du gradient de la transformation \({\underset{\tilde{}}{\mathbf{F}}}\) est adoptée selon :

\[ {\underset{\tilde{}}{\mathbf{F}}} = {\underset{\tilde{}}{\mathbf{F}}}_{e}\,\cdot\, {\underset{\tilde{}}{\mathbf{F}}}_{p} \]

avec la partie plastique \({\underset{\tilde{}}{\mathbf{F}}}_{p}\) pour l’écoulement plastique induit par les glissements entre différents systèmes de glissement, et la partie élastique \({\underset{\tilde{}}{\mathbf{F}}}_{e}\) pour la déformation élastique pure et la rotation de corps rigide du réseau cristallin. De plus, la déformation élastique de Green-Lagrange est définie par :

\[ \underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}} ={{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}} \left( {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{T}}} {\underset{\tilde{}}{\mathbf{F}}}_{e} - \underline{I} \right) \]

Le gradient du vecteur vitesse \({\underset{\tilde{}}{\mathbf{L}}}\) est défini par :

\[ {\underset{\tilde{}}{\mathbf{L}}} = \vec{\nabla} \vec{v} = \dot{{\underset{\tilde{}}{\mathbf{F}}}} {{\underset{\tilde{}}{\mathbf{F}}}^{\mathrm{-1}}} \]

où \(\vec{v}\) est le vecteur vitesse. Compte tenu de la décomposition multiplicative du gradient de transformation, le gradient du vecteur vitesse \({\underset{\tilde{}}{\mathbf{L}}}\) s’exprime :

\[ {\underset{\tilde{}}{\mathbf{L}}} = \dot{{\underset{\tilde{}}{\mathbf{F}}}} {{\underset{\tilde{}}{\mathbf{F}}}^{\mathrm{-1}}} = \dot{{\underset{\tilde{}}{\mathbf{F}}}_{e}} {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{-1}}} + {\underset{\tilde{}}{\mathbf{F}}}_{e} \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{-1}}} = {\underset{\tilde{}}{\mathbf{L}}}^e + {\underset{\tilde{}}{\mathbf{L}}}^p \]

avec \({\underset{\tilde{}}{\mathbf{L}}}^e=\dot{{\underset{\tilde{}}{\mathbf{F}}}_{e}} {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{-1}}}\) et \({\underset{\tilde{}}{\mathbf{L}}}^p={\underset{\tilde{}}{\mathbf{F}}}_{e} \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{-1}}}\).

Pour la plasticité cristalline, \(\dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}}\) peut être déterminé sur les 12 systèmes de glissement des matériaux CFC de type \(\{111\}\left<110\right>\) par :

\[ \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} = \sum \limits_{s=1}^{12} \dot{\gamma}^s {\underset{\tilde{}}{\mathbf{N}}}^s \qquad{(1)}\]

avec le tenseur d’orientation \({\underset{\tilde{}}{\mathbf{N}}}^s\) pour le système de glissement \(s\) et le taux de glissement \(\dot{\gamma}^s\) pour le système de glissement \(s\).

Dans la suite, la formulation thermodynamique est présentée pour le monocristal. Le deuxième principe de la thermodynamique dans sa forme locale, qui est connu sous le nom d’Inégalité de Clausius-Duhem s’écrit :

\[ -\rho (\dot e - T \dot s) + \underline{\sigma}\,\colon\, \underline{D} - {{\displaystyle \frac{\displaystyle \vec{q}}{\displaystyle T}}} \cdot \vec{\nabla} T \geqslant 0 \qquad{(2)}\]

avec \(\rho\) la masse volumique du monocristal dans la configuration actuelle, \(e\) la densité d’énergie interne, \(T\) la température, \(s\) la densité d’entropie, \(\underline{D}\) le tenseur vitesse de déformation \(\underline{D} = {{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}}{\left({\underset{\tilde{}}{\mathbf{L}}} + {{\underset{\tilde{}}{\mathbf{L}}}^{\mathrm{T}}}\right)}\) et \(\vec{q}\) le flux de chaleur. Pour le monocristal, la puissance des efforts intérieurs est:

\[ \begin{aligned} {{\displaystyle \frac{\displaystyle 1}{\displaystyle \rho}}} \underline{\sigma}\,\colon\, \underline{D} &=& {{\displaystyle \frac{\displaystyle 1}{\displaystyle \rho}}}\underline{\sigma}\,\colon\, {\underset{\tilde{}}{\mathbf{L}}} = {{\displaystyle \frac{\displaystyle 1}{\displaystyle \rho}}}\underline{\sigma}\,\colon\, \dot{{\underset{\tilde{}}{\mathbf{F}}}} {{\underset{\tilde{}}{\mathbf{F}}}^{\mathrm{-1}}} \\ &=& {{\displaystyle \frac{\displaystyle 1}{\displaystyle \rho}}}\underline{\sigma}\,\colon\, (\dot{{\underset{\tilde{}}{\mathbf{F}}}_{e}} {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{-1}}}+{\underset{\tilde{}}{\mathbf{F}}}_{e} \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{-1}}}) \\ &=& {{\displaystyle \frac{\displaystyle 1}{\displaystyle \rho}}}\underline{\sigma}\,\colon\, (\dot{{\underset{\tilde{}}{\mathbf{F}}}_{e}} {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{-1}}}) + {{\displaystyle \frac{\displaystyle 1}{\displaystyle \rho}}}\underline{\sigma}\,\colon\, ({\underset{\tilde{}}{\mathbf{F}}}_{e} \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{-1}}}) \end{aligned} \qquad{(3)}\]

Dans l’équation (3), comme \(\underline{\sigma}\) est un tenseur symétrique, on peut donc remplacer \(\underline{D}\) par \({\underset{\tilde{}}{\mathbf{L}}}\). Dans l’équation (¿eq:power_internal_forces_1?), tenant compte de la relation entre le second tenseur de Piola-Kirchhoff \(\underline{S}^e\) et le tenseur de Cauchy \(\underline{\sigma}\), on a

\[ \begin{aligned} \frac{1}{\rho} \underline{\sigma}\,\colon\, \underline{D} &=& \frac{1}{\rho_i} ({\underset{\tilde{}}{\mathbf{F}}}_{e} \underline{S}^e {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{T}}})\,\colon\, (\dot{{\underset{\tilde{}}{\mathbf{F}}}_{e}} {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{-1}}}) + \frac {1}{\rho_i} ({\underset{\tilde{}}{\mathbf{F}}}_{e} \underline{S}^e {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{T}}})\,\colon\, ({\underset{\tilde{}}{\mathbf{F}}}_{e} \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{-1}}}) \\ &=& \frac{1}{\rho_i} \underline{S}^e\,\colon\, ({{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{T}}} \dot{{\underset{\tilde{}}{\mathbf{F}}}_{e}}) + \frac{1}{\rho_i} ({{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{T}}} {\underset{\tilde{}}{\mathbf{F}}}_{e} \underline{S}^e)\,\colon\, (\dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}}) \end{aligned} \]

Comme le second tenseur de Piola-Kirchhoff est symétrique, on a :

\[ \begin{aligned} \frac{1}{\rho} \underline{\sigma}\,\colon\, \underline{D} &=& \frac{1}{\rho_i} \underline{S}^e\,\colon\, {{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}}\left( {{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{T}}} \dot{{\underset{\tilde{}}{\mathbf{F}}}_{e}} + {\dot{{\underset{\tilde{}}{\mathbf{F}}}_{e}}^{\mathrm{T}}} {\underset{\tilde{}}{\mathbf{F}}}_{e}\right)+ \frac{1}{\rho_i} ({{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{T}}} {\underset{\tilde{}}{\mathbf{F}}}_{e} \underline{S}^e)\,\colon\, (\dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}}) \\ &=& \frac{1}{\rho_i} \underline{S}^e\,\colon\, \underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}}+ \frac{1}{\rho_i} {\underset{\tilde{}}{\mathbf{M}}}\,\colon\, (\dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}}) \end{aligned} \qquad{(4)}\]

Pour le monocristal sans écrouissage, la densité d’énergie interne \(e\) est une fonction de la déformation élastique de Green-Lagrange \(\underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}}\) et de la densité d’entropie \(s\) :

\[ e=e \left( \underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}},s \right) \]

la densité d’énergie libre \(\psi\) est une fonction de \(\underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}}\) et de la température \(T\) :

\[ \psi=\psi \left( \underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}},T \right) \]

On sait que la relation entre la densité d’énergie libre \(\psi\) et la densité d’énergie interne \(e\) s’écrit :

\[ \psi=e-Ts \qquad{(5)}\]

Compte tenu des équations (4) et (5), l’inégalité de Clausisu-Duhem devient :

\[ \rho \left( \frac{\underline{S}^{e}}{\rho_i} - \frac{\partial \psi}{\partial \underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}}} \right)\,\colon\, \underline{\dot{\epsilon}}^{\mathrm{el}}_{\mathrm{GL}} - \rho \left( s + \frac{\partial \psi}{\partial T}\right) \dot T +\frac{\rho}{\rho_i} {\underset{\tilde{}}{\mathbf{M}}}\,\colon\, \left( \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} \right) - \frac{\vec{q}}{T} \cdot \vec{\nabla} T \geqslant 0 \qquad{(6)}\]

où on a la dissipation thermique \(D^{t\!h}\) :

\[ D^{t\!h} = - \frac{\vec{q}}{T} \cdot \vec{\nabla} T \]

et la dissipation intrinsèque \(D^i\):

\[ D^i = \rho \left( \frac{\underline{S}^{e}}{\rho_i} - \frac{\partial \psi}{\partial \underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}}} \right)\,\colon\, \underline{\dot{\epsilon}}^{\mathrm{el}}_{\mathrm{GL}} - \rho \left( s + \frac{\partial \psi}{\partial T}\right) \dot T +\frac{\rho}{\rho_i} {\underset{\tilde{}}{\mathbf{M}}}\,\colon\, \left( \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} \right) \]

Considérons un processus réversible tel que \(\dot T = 0\), \(\dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}}=0\) et \(\vec{\nabla} T = 0\). Pour le processus réversible, l’inégalité \(\ref{eq:Clausius-Duhem_crystal}\) réduit à l’égalité, puisqu’il n’y a pas de dissipation. Par conséquent, on obtient

\[ \underline{S}^e = \rho_i \frac{\partial \psi}{\partial \underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}}} \]

De même, considérons un autre processus réversible tel que \(\underline{\dot{\epsilon}}^{\mathrm{el}}_{\mathrm{GL}}=0\), \(\dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}}=0\) et \(\vec{\nabla} T = 0\), on a

\[ s=-\frac{\partial \psi}{\partial T} \]

Enfin, la dissipation intrinsèque se réduit à :

\[ D^i = \frac{\rho}{\rho_i} {\underset{\tilde{}}{\mathbf{M}}}\,\colon\, \left( \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} \right) \qquad{(7)}\]

Dans l’équation (7), \(\rho / \rho_i\) vient du fait que l’inégalité \(\ref{eq:Clausius-Duhem}\) est écrite pour une volume \(dV\) dans la configuration actuelle et que \({\underset{\tilde{}}{\mathbf{M}}} : \left( \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} \right)\) est pour une volume \(dV_i\) dans la configuration intermédiaire. En raison de la conservation de masse, on a

\[ dV_i = \frac{\rho}{\rho_i} dV \]

Dans le cadre de matériaux standards généralisés, un potentiel convexe \(\Omega= \Omega\left( {\underset{\tilde{}}{\mathbf{M}}} \right)\) existe et permet de définir l’évolution de \(\dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}}\) :

\[ \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} = \frac{\partial \Omega}{\partial {\underset{\tilde{}}{\mathbf{M}}}} \]

Selon les équations (4) et (7), le second tenseur de Piola-Kirchhoff \(\underline{S}^e\) est utilisé pour la loi d’élasticité et le tenseur de Mandel \({\underset{\tilde{}}{\mathbf{M}}}\) pour la loi d’écoulement. Vu que la dissipation décrite par \({\underset{\tilde{}}{\mathbf{M}}}\,\colon\, \left( \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} \right)\) est liée au glissement des dislocations dans le monocristal, on a

\[ {\underset{\tilde{}}{\mathbf{M}}}\,\colon\, \left( \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} \right) = \sum_{s=1}^{12} \dot \gamma^s \tau^s \qquad{(8)}\]

où \(\tau^s\) est la cission résolue du système de glissement \(s\). Compte tenu l’équation (1), on a

\[ {\underset{\tilde{}}{\mathbf{M}}}\,\colon\, \left( \dot{{\underset{\tilde{}}{\mathbf{F}}}_{p}} {{\underset{\tilde{}}{\mathbf{F}}}_{p}^{\mathrm{-1}}} \right) = {\underset{\tilde{}}{\mathbf{M}}}\,\colon\, \left(\sum_{s=1}^{12} \dot \gamma^s {\underset{\tilde{}}{\mathbf{N}}}^s \right) = \sum_{s=1}^{12} \dot \gamma^s {\underset{\tilde{}}{\mathbf{M}}}\,\colon\, {\underset{\tilde{}}{\mathbf{N}}}^s \qquad{(9)}\]

Si on compare l’équation (8) et l’équation (9), la cission résolue \(\tau^s\) du système de glissement \(s\) s’écrit donc :

\[ \tau^s = {\underset{\tilde{}}{\mathbf{M}}}\,\colon\, {\underset{\tilde{}}{\mathbf{N}}}^s \]

The FiniteStrainSingleCrystal brick

The previous framework allows the implementation of a wide range of finite strain single crystals that differ from their flow rules and the constitutive equations used to describe their internal state variables’ evolution.

The FiniteStrainSingleCrystal brick hides all the details related to this framework and simplifies the declaration of the slip systems and the interaction matrix.

Variables automatically defined

The plastic slip \(g\) are automatically defined.

Implementation of the implicit system

Le second tenseur de Piola-Kirchhoff est défini par la loi élastique de Saint-Venant Kirchhoff: \[ \underline{S}=\underline{\mathbf{D}}\,.\,{\left({\left.\underline{\epsilon}^{\mathrm{el}}\right|_{t}}+\Delta\,\underline{\epsilon}^{\mathrm{el}}\right)} \]

Nous calculons (une approximation) de l’inverse de l’incrément de la partie plastique du gradient de la transformation: \[ {\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1} = {\underset{\tilde{}}{\mathbf{I}}}-\sum_{i=1}^{12}\Delta\,g_{i}\,{\underset{\tilde{}}{\mathbf{\mu}}}_{i} \]

Nous calculons la partie élastique de gradient de la transformation \({\underset{\tilde{}}{\mathbf{F}}}_{e}\) en fin de pas de temps en fin de pas de temps: \[ {\underset{\tilde{}}{\mathbf{F}}}_{e}={\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}\,\star\,{\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1} \] où \(\star\) représente le produit des représentations matricielles de \({\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}\) et \({\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}\).

L’équation associée à l’incrément de déformation élastique impose qu’elle soit égale au tenseur de Green-Lagrange \(\underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}}\) associé à la partie élastique \({\underset{\tilde{}}{\mathbf{F}}}_{e}\) du gradient de la transformation: \[ {\left.\underline{\epsilon}^{\mathrm{el}}\right|_{t}}+\Delta\,{\left.\underline{\epsilon}^{\mathrm{el}}\right|_{t}}- \underline{\epsilon}^{\mathrm{el}}_{\mathrm{GL}}{\left({\underset{\tilde{}}{\mathbf{F}}}_{e}\right)}=0 \]

Le calcul des termes \({{\displaystyle \frac{\displaystyle \partial f_{\underline{\epsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,g_{i}}}}\) est assez technique et devra se faire en plusieurs étapes: \[ {{\displaystyle \frac{\displaystyle \partial \underline{\varepsilon_{\mathrm{GL}}}}{\displaystyle \partial \Delta\,g_{i}}}}= {{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}}\,{{\displaystyle \frac{\displaystyle \partial \underline{C}_{e}}{\displaystyle \partial \Delta\,g_{i}}}}= {{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}}\,{{\displaystyle \frac{\displaystyle \partial \underline{C}_{e}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}}\,.\, {{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}{\displaystyle \partial \Delta\,g_{i}}}} \]

La classe t2tost2 fournit une méthode statique permettant de calculer la dérivée du tenseur de Cauchy élastique \(\underline{C}_{e}\) par rapport au gradient de la transformation élastique \(F_{e}\).

Il faut alors calculer le terme \({{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}{\displaystyle \partial \Delta\,g_{i}}}}\). Le tenseur \({\underset{\tilde{}}{\mathbf{F}}}_{e}\) est le produit \(\star\) des représentation matricielle de deux tenseurs: \({\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}\) et le tenseur \({\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}\).

Cette démarche est maintenant appliquée au calcul de la dérivée du tenseur \({\underset{\tilde{}}{\mathbf{F}}}_{e}\): \[ {{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}{\displaystyle \partial {\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}}}} = \partial^{r}_{\star}{\left({\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}\right)} \]

En négligeant la correction plastique, la dérivée \({{\displaystyle \frac{\displaystyle \partial {\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}}{\displaystyle \partial \Delta\,g_{i}}}}\) s’écrit simplement: \[ {{\displaystyle \frac{\displaystyle \partial {\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}}{\displaystyle \partial \Delta\,g_{i}}}} = -{\underset{\tilde{}}{\mathbf{\mu}}}_{i} \]

Nous obtenons finalement l’expression la dérivée \({{\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \Delta\,g_{i}}}}\): \[ {{\displaystyle \frac{\displaystyle \partial f_{\underline{\epsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,g_{i}}}}= -{{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}}\,{{\displaystyle \frac{\displaystyle \partial \underline{C}_{e}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}}\,.\, \partial^{r}_{\star}{\left({\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}\right)}\,.\, {\underset{\tilde{}}{\mathbf{\mu}}}_{i} \]

Implementation of the consistent tangent operator

Le point de départ du calcul est la relation: \[ {{\displaystyle \frac{\displaystyle \partial \underline{\tau}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}={{\displaystyle \frac{\displaystyle \partial \underline{\tau}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}}{{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}} \]

Nous calculons maintenant chacun des termes du produit.

La contrainte de Kirchhoff \(\underline{\tau}\) se calcule par le transport dans la configuration finale du second tenseur de Piola-Kirchhoff calculé dans la configuration intermédiaire. \[ \underline{\tau} = {\underset{\tilde{}}{\mathbf{F}}}_{e}\,\star\,\underline{S}\,\star\,{{\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{T}}} \] où \(\star\) représente le produit matriciel.

TFEL propose la méthode computePushForwardDerivative permettant de relier la dérivée \({{\displaystyle \frac{\displaystyle \partial \underline{\tau}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}}\) à la dérivée du second tenseur de Piola-Kirchhoff \({{\displaystyle \frac{\displaystyle \partial \underline{S}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}}\). Cette dérivée se calcule par dérivation en chaînes: \[ {{\displaystyle \frac{\displaystyle \partial \underline{S}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}}={{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}}\,D\,{{\displaystyle \frac{\displaystyle \partial \underline{C}_{e}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}} \]

La dérivée \({{\displaystyle \frac{\displaystyle \partial \underline{C}_{e}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}}\) se calcule par la méthode statique dCdF de la classe t2tost2.

Il nous fait maintenant calculer le terme \({{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\).

Nous savons que: \[ {\underset{\tilde{}}{\mathbf{F}}}_{e}={\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}\,\star\,{\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1} \] Par dérivation, nous obtenons: \[ {{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}= \partial^{l}_{\star}{\left({\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}\right)}\,{{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}+ \partial^{r}_{\star}{\left({\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}\right)}\,{{\displaystyle \frac{\displaystyle \partial {\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}} \]

Le calcul du terme \({{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\) est aisé: \[ {{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}= \partial^{l}_{\star}{\left({\left.{\underset{\tilde{}}{\mathbf{F}}}_e\right|_{t}}\right)} \]

Pour le calcul de la dérivée \({{\displaystyle \frac{\displaystyle \partial {\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\), nous pouvons nous baser sur l’approximation faite lors de l’intégration. Nous avons alors directement: \[ {{\displaystyle \frac{\displaystyle \partial {\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}= -\sum_{i=1}^{12}{\underset{\tilde{}}{\mathbf{\mu}}}_{i}\,\otimes\,{{\displaystyle \frac{\displaystyle \partial \Delta\,g_{i}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}} \]

Le calcul de \({{\displaystyle \frac{\displaystyle \partial \Delta\,g_{i}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\) représente la seule vraie difficulté. Son calcul peut se faire en généralisant les techniques utilisées en petites déformations.

Pour toute variation de \(\delta\,\Delta\,{\underset{\tilde{}}{\mathbf{F}}}\) de \(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}\), la variation du système implicite \(\delta\,F\) est nulle: \[ \delta\,F{\left(\Delta\,Y{\left(\delta\,\Delta\,{\underset{\tilde{}}{\mathbf{F}}}\right)},\delta\,\Delta\,{\underset{\tilde{}}{\mathbf{F}}}\right)} =J\,\delta\,\Delta\,Y+ \left.{{\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\right|_{{\left(\Delta\,\underline{\epsilon}^{\mathrm{el}},\Delta\,g_{i}\right)}}\,\delta\,\Delta\,{\underset{\tilde{}}{\mathbf{F}}}=0 \]

Nous avons explicitement spécifié, en utilisant la notation classique de la thermodynamique, que le calcul du terme \(\left.{{\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\right|_{{\left(\Delta\,\underline{\epsilon}^{\mathrm{el}},\Delta\,g_{i}\right)}}\) se faisait en supposant que \(\Delta\,\underline{\epsilon}^{\mathrm{el}}\) et \(\Delta\,g_{i}\) étaient des variables indépendantes.

Il est donc nécessaire de calculer la dérivée \(\left.{{\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\right|_{{\left(\Delta\,\underline{\epsilon}^{\mathrm{el}},\Delta\,g_{i}\right)}}\) dont la seule composante non nulle est \(\left.{{\displaystyle \frac{\displaystyle \partial f_{\underline{\epsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\right|_{{\left(\Delta\,\underline{\epsilon}^{\mathrm{el}},\Delta\,g_{i}\right)}}\): \[ \left.{{\displaystyle \frac{\displaystyle \partial f_{\underline{\epsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\right|_{{\left(\Delta\,\underline{\epsilon}^{\mathrm{el}},\Delta\,g_{i}\right)}}= -{{\displaystyle \frac{\displaystyle \partial \underline{C}_{e}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}}\,.\, {{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}= -{{\displaystyle \frac{\displaystyle \partial \underline{C}_{e}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}}\,.\, {{\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}{\left({\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}\,\star\,{\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}\right)} \]

\[ {{\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}{\left({\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}\,\star\,{\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}\right)}= \partial^{l}_{\star}{\left({\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}\right)}\,.\, {{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}^{\mathrm{tr}}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}= \partial^{l}_{\star}{\left({\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}\right)}\,.\, \partial^{l}_{\star}{\left({\left.{\underset{\tilde{}}{\mathbf{F}}}_{e}\right|_{t}}\right)} \]

Finalement, \({{\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\) s’écrit: \[ {{\displaystyle \frac{\displaystyle \partial f_{\underline{\epsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}={{\displaystyle \frac{\displaystyle \partial \underline{C}_{e}}{\displaystyle \partial {\underset{\tilde{}}{\mathbf{F}}}_{e}}}}\,.\, \partial^{l}_{\star}{\left({\left(\Delta\,{\underset{\tilde{}}{\mathbf{F}}}_{p}\right)}^{-1}\right)}\,.\, \partial^{l}_{\star}{\left({\left.{\underset{\tilde{}}{\mathbf{F}}}_{e}\right|_{t}}\right)} \]

La dérivée \({{\displaystyle \frac{\displaystyle \partial \Delta\,g_{i}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}\) est alors: \[ {{\displaystyle \frac{\displaystyle \partial \Delta\,g_{i}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}}=-J^{-1}_{\Delta\,g_{i}}|{{\displaystyle \frac{\displaystyle \partial f_{\underline{\epsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,{\underset{\tilde{}}{\mathbf{F}}}}}} \]

Code blocks automatically declared

The ComputeFinalStress code block is automatically declared.

The TangentOperator code block is automatically declared.

Variables automatically computed

The residual \(f_{\underline{\epsilon}^{\mathrm{el}}}\) associated with the elastic strain \(\underline{\epsilon}^{\mathrm{el}}\) is automatically implemented.

The Mandel stress \({\underset{\tilde{}}{\mathbf{M}}}\) is automatically computed before the integrator code block. Its value is stored in the variables M. If the algorithm requires an analytic jacobian, the derivative \({{\displaystyle \frac{\displaystyle \partial {\underset{\tilde{}}{\mathbf{M}}}}{\displaystyle \partial \Delta\,\underline{\epsilon}^{\mathrm{el}}}}}\) is also computed and stored in the variable dM_ddeel respectively.

References