r5.03.11 Comportements élastoviscoplastiques mono et polycristallins

Résumé :

Le but de ce document est de décrire l’intégration des comportements mono et polycristallins.

On traite ici de l’intégration de ces lois de comportement associées à des systèmes de glissement correspondant aux familles cristallines habituelles ou spécifiées par l’utilisateur. Cette intégration peut se faire de façon explicite (méthode de Runge-Kutta avec contrôle de la précision et redécoupage local du pas de temps) ou de façon implicite (méthode de Newton avec redécoupage local du pas de temps).

Ces comportements peuvent être employés pour le calcul de microstructures (maillage d’un agrégat, avec représentation géométrique de chaque grain physique monocristallin) ou pour le calcul de type POLYCRISTAL, milieu «homogénéisé» possédant en tout point matériel (ou point d’intégration ou de calcul) plusieurs phases simultanées, dans des proportions variables.

On décrit aussi les types de calculs possibles avec ces comportements ainsi que la procédure de génération de maillages.

Formulation des comportements mono et polycristallins

Principe des lois de comportement du monocristal

Les comportements relatifs aux système de glissement d’un monocristal sont (dans l’ensemble des comportements envisagés) de type élasto-visco-plastique. Il peuvent être construits sur des bases phénoménologiques ou sur des bases physiques (dynamique des dislocations). Pour chacune des directions de glissement, le comportement est mono dimensionnel . Il peut se décomposer en trois types d’équations:

• Relation d’écoulement: \(\dot{{\gamma}_{s}}\) représente le glissement plastique du système \(s\)

\(\dot{{\gamma}_{s}}=g({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\) , avec \(\dot{{p}_{s}}=∣\dot{{\gamma}_{s}}∣\) et:

• cas élastoplastique, un critèredu type: \(F({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\le 0\) et \(F\cdot \dot{{p}_{s}}=0\)

• cas élasto-viscoplastique, \(\dot{{p}_{s}}=f({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\)

• \(\dot{{\alpha}_{s}}\) représente l’écrouissage cinématique dans le cas d’une loi phénoménologique et la densité de dislocations pour un comportement à base physique.

Son évolution est la forme: \(\dot{{\alpha}_{s}}=h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\)

• pour les lois phénoménologiques, l’écrouissage isotropeest défini par une fonction: \(R({p}_{s})\) .

Ces relations deviennent, après discrétisation en temps:

• Relation d’écoulement:

\(\Delta {\gamma}_{s}=g({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\) , avec, \(\Delta {p}_{s}=∣\Delta {\gamma}_{s}∣\)

pour un comportement élastoplastique un critèredu type: \(F({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\le 0\) et \(F\cdot \Delta {p}_{s}=0\)

et pour un comportement élasto-viscoplastique, \(\Delta {p}_{s}=f({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\Delta t\)

• Évolution de \(\dot{{\alpha}_{s}}\) : \(\Delta {\alpha}_{s}=h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\)

• Évolution de l’écrouissage isotrope: \(R({p}_{s})\)

Les quantités \(({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\) sont évaluées à l’instant courant pour une discrétisation implicite et à l’instant précédent pour une discrétisation explicite.

Les lois monocristallines disponibles sont décrites ci-après. Les noms de ces relations correspondent à leur appellation dans la commande DEFI_MATERIAU [U4.43.01].

Lois de comportement du monocristal disponibles

Lois phénoménologiques

Elles sont construites à partir d’une loi d’écoulement, une loi d’écrouissage isotrope et une loi d’écrouissage cinématique. La combinaison des lois MONO_VISC1, MONO_ISOT1, MONO_CINE1 correspond à la loi cristalline connue sous le nom de Meric-Cailletaud 1 .

Loi d’écoulement MONO_VISC1

\(\begin{array}{}\Delta {\gamma}_{s}=g({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})=\Delta {p}_{s}\frac{{\tau}_{s}-c{\alpha}_{s}}{∣{\tau}_{s}-c{\alpha}_{s}∣}\\ \Delta {p}_{s}=\Delta t{\langle \frac{∣{\tau}_{s}-c{\alpha}_{s}∣-{R}_{s}({p}_{s})}{k}\rangle }^{n}\end{array}\)

Les paramètres sont: \((c,k,n)\) .

Loi d’écoulement MONO_VISC2

\(\begin{array}{c}\Delta {\gamma}_{s}=g({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})=\Delta {p}_{s}\frac{{\tau}_{s}-c{\alpha}_{s}-a{\gamma}_{s}}{∣{\tau}_{s}-c{\alpha}_{s}-a{\gamma}_{s}∣}\\ \Delta {p}_{s}=\Delta t{\langle \frac{∣{\tau}_{s}-c{\alpha}_{s}-a{\gamma}_{s}∣-{R}_{s}({p}_{s})+\frac{d}{\mathrm{2c}}{(c{\alpha}_{s})}^{2}}{k}\rangle }^{n}\end{array}\)

Les paramètres sont : \(c,k,n,a,d\) .

Loi d’écrouissage cinématique MONO_CINE1

\(\Delta {\alpha}_{s}=h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})=\Delta {\gamma}_{s}-d{\alpha}_{s}\Delta {p}_{s}\)

Le paramètre est : \(d\) .

Loi d’écrouissage cinématique MONO_CINE2

\(\Delta {\alpha}_{s}=h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})=\Delta {\gamma}_{s}-d{\alpha}_{s}\Delta {p}_{s}-{(\frac{∣c{\alpha}_{s}∣}{M})}^{m}\frac{{\alpha}_{s}}{∣{\alpha}_{s}∣}\)

Les paramètres étant alors : \(d,M,m,c\) .

Loi d’écrouissage isotrope MONO_ISOT1

\({R}_{s}(p)={R}_{0}+Q(\sum_{r=1}^{\mathrm{ns}}{h}_{\mathrm{sr}}(1-{e}^{{\mathrm{bp}}_{r}}))\) ; les paramètres sont: \({R}_{0,}Q,b,{h}_{\mathrm{sr}}\) .

Loi d’écrouissage isotrope MONO_ISOT2

\({R}_{s}(p)={R}_{0}+{Q}_{1}\sum_{r}{h}_{\mathrm{sr}}(1-{e}^{-{b}_{1}{p}_{r}})+{Q}_{2}(1-{e}^{-{b}_{2}{p}_{s}})\) ; les paramètres sont: \({R}_{0,}{Q}_{1,}{b}_{1,}{b}_{2,}{h}_{\mathrm{sr}},{Q}_{2}\) .

Matrice d’interaction

\({h}_{\mathrm{sr}}\) une matrice carrée symétrique d’interaction entre systèmes de glissement (voir annexe 4).

Dans le cas où un seul un coefficient \(H\) est fourni, la matrice \({h}_{\mathrm{sr}}\) prend une forme très simple :

\({h}_{\mathrm{rs}}=1\mathrm{si}r=s\) et \({h}_{\mathrm{rs}}=H\mathrm{si}r\ne s\) , avec par défaut \(H=0\) (matrice identité).

Loi issue de la Dynamique des Dislocations : MONO_DD_KR

Le modèle MONO_DD_KR (Kocks-Rauch) est applicable aux aciers de type \(\mathit{CC}\) (structure cristalline Cubique Centré). Sa variable interne principale est la densité de dislocations sur chaque système de glissement. Il offre par ailleurs l’avantage d’être valide pour une large gamme de températures [bib7]. L’écoulement se formule de la façon suivante:

Si \(∣{\tau}^{s}∣>{\tau}_{0}\)

\({\tau}_{\mathrm{eff}}^{s}=\langle ∣{\tau}^{s}∣-{\tau}_{0}-{\tau}_{\mu}^{s}\rangle\) \(\langle \rangle\) partie positive

\({\tau}_{\mu}^{s}=\frac{{(\mu )}^{2}{\sum}_{u}{a}^{\mathrm{su}}{\alpha}^{u}}{∣{\tau}^{s}∣-{\tau}_{0}}\) avec \({\alpha}^{u}={\rho}^{u}{b}^{2}\)

sinon \({\tau}_{\mathrm{eff}}^{s}=0\) . Dans ces expressions, \(\mu\) désigne le module de cisaillement du matériau, \(b\) est la constante de Burgers, \({\rho}^{u}\) est la densité de dislocations, \({a}^{\mathrm{su}}\) est la matrice d’interaction entre les 24 systèmes de glissement (sa forme particulière, dépendant de 5 coefficients, est fournie en annexe). Les termes diagonaux représentent l’auto-écrouissage d’un système de glissement, les autres termes représentent l’écrouissage latent.

Si \(∣{\tau}^{s}∣>{\tau}_{0}\) et \({\tau}_{\mathrm{eff}}^{s}>0\)

\(\Delta {\gamma}_{s}=g({\tau}_{s},{\gamma}_{s})=\dot{{\gamma}_{0}}\exp(\frac{-\Delta G({\tau}_{\mathrm{eff}}^{s})}{{k}_{B}T})\cdot \frac{{\tau}_{s}}{∣{\tau}_{s}∣}\cdot \Delta t\)

avec \(\Delta G({\tau}_{\mathit{eff}}^{s})=\Delta {G}_{0}{(1-{(\frac{\langle {\tau}_{\mathit{eff}}^{s}\rangle }{{\tau}_{R}})}^{p})}^{q}\)

sinon

\(\Delta {\gamma}_{s}=0\)

\(\Delta {\alpha}_{s}=0\)

\(\Delta G({\tau}_{\mathrm{eff}}^{s})\) est une activation d’enthalpie correspondant à une activation libre dans un domaine athermique. \({\gamma}_{0,}{\tau}_{R},p,q\) sont des données matériaux, \({k}_{B}\) la constante de Boltzman et \(T\) la température.

L’évolution de la densité de dislocation est proportionnelle à la racine carrée de la somme des densités de dislocations de tous les autres systèmes de glissement :

\(\dot{{\rho}_{s}}=\frac{∣\dot{{\gamma}^{s}}∣}{b}(\frac{1}{{\Lambda}^{s}}-{g}_{c}(T){\rho}^{s})\) avec \(\frac{1}{{\Lambda}^{s}}=\frac{1}{d}+\frac{\sqrt{\sum_{u\ne s}{\rho}^{u}}}{K(T)}\)

et \({g}_{c}(T)={g}_{\mathrm{c0}}\exp\left[-\frac{{E}_{\mathrm{gc}}}{{k}_{B}T}\right]\) où \({E}_{\mathrm{gc}},{k}_{B},b,d,{g}_{\mathrm{c0}},K(T)\) sont des paramètres matériaux,

d’où

\(\Delta {\rho}^{s}=∣\Delta {\gamma}^{s}∣(\frac{1}{\mathrm{bd}}+\frac{\sqrt{\sum_{u\ne s}{\alpha}^{u}}}{{b}^{2}K}-\frac{{g}_{c}{\alpha}^{s}}{{b}^{3}})\) et par conséquent:

\(\Delta {\alpha}^{s}=∣\Delta {\gamma}^{s}∣(\frac{b}{d}+\frac{\sqrt{\sum_{u\ne s}{\alpha}^{u}}}{K}-\frac{{g}_{c}{\alpha}^{s}}{b})=∣g({\tau}_{s},{\alpha}_{s})∣\cdot h({\alpha}_{s})\)

Les paramètres sont: \(\dot{{\gamma}_{0,}}{k}_{B},T,{\tau}_{R},\Delta {G}_{0,}p,q,{\tau}_{0,}\mu ,b,d,K,{g}_{c}\)

\({k}_{B}\) est la constante de Boltzman.

Remarque: La résolution de Newton peut amener à des valeurs négatives pour \(\Delta {\alpha}^{u}\) , il faut donc que les termes sous la racine soient protégés. Pour ce faire, nous choisissons la modification suivante: \(\sqrt{\sum_{u\ne s}\langle {\alpha}^{u}\rangle }\) où \(\langle .\rangle\) représente la parti positive : \(\langle x\rangle =\lbrace \begin{array}{}x\text{si}x>0\\ 0\text{sinon}\end{array}\)

Remarquons aussi que l’expression de \(\Delta {\gamma}_{s}\) , qui est fonction de \(\Delta G({\tau}_{\mathrm{eff}}^{s})\) , conduit à une discontinuité au voisinage de \({\tau}_{\mathrm{eff}}^{s}=0\) : pour \({\tau}_{\mathrm{eff}}^{s}<0\) , \(\Delta {\gamma}_{s}=0\) , et \(\underset{{\tau}_{\mathrm{eff}}^{s}\to 0\text{*}}{\lim}\Delta {\gamma}_{s}=\dot{{\gamma}_{0}}\exp(\frac{-\Delta {G}_{0}}{{k}_{B}T})\cdot \frac{{\tau}_{s}}{∣{\tau}_{s}∣}\cdot \Delta t\)

Lois issue de la Dynamique des Dislocations : MONO_DD_CFC / MONO_DD_CFC_IRRA /MONO_DD_FAT

La formulation des lois DD_CFC [8], [9] et DD_FAT [18] est construite à partir de calculs de dynamique des dislocations. Elle s’applique aux matériaux à structure cristalline Cubique à Faces Centrées (\(\mathit{CFC}\) ) tels que les aciers austénitiques. A priori le modèle DD_CFC n’est pas compatible avec un changement de trajet de chargement, en particulier les paramètres ne sont pas adaptés à une sollicitation cyclique (lorsqu’une dislocation “revient sur ses pas” sa cinétique est différente du fait d’une interaction différente avec les obstacles).

Les équations ci-dessous sont écrites pour chaque système de glissement \(s\) .

La direction de l’écoulement est donnée par :

\(\dot{{\gamma}_{s}}=\dot{{p}_{s}}\frac{{\tau}_{s}}{∣{\tau}_{s}∣}=\dot{{p}_{s}}\text{signe}({\tau}_{s})\) et l’intensité de l’écoulement par :

\(\dot{{p}_{s}}=\dot{{\gamma}_{0}}({(\frac{∣{\tau}_{s}∣}{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}})}^{n}-1)\text{si}∣{\tau}_{s}∣\ge {\tau}_{f}+{\tau}_{s}^{\mathrm{forest}},\text{sinon}\dot{{p}_{s}}=0\)

avec \(n\) grand (pour représenter un modèle quasi élasto-plastique : pour le CFC, la dépendance en temps est en effet négligeable).

Par rapport au modèle initial [8], le seuil a été introduit ainsi que le terme \(–\dot{{\gamma}_{0}}\) pour assurer la continuité lorsque le seuil est franchi.

En pratique, le paramètre \(\dot{{\gamma}_{0}}\) est petit et le terme \({(\frac{∣{\tau}_{s}∣}{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}})}^{n}\) devient très vite largement supérieur à \(1\) , à cause des grandes valeurs de l’exposant \(n\) .

\({\tau}_{s}^{\mathrm{forest}}\) est une fonction d’écrouissage dont les coefficients sont fournis par le mot-clé MONO_DD_* (*=CFC,FAT ou CC), et qui est décrite au § suivant.

Dans DD_CFC, l’évolution de la densité de dislocation est donnée par :

\(\dot{{\rho}_{s}}=\frac{\dot{{p}_{s}}}{b}(A\frac{\sum_{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}^{\mathrm{eff}}}{\rho}_{j}}{\sum_{j=1,12}\sqrt{{a}_{\mathrm{sj}}^{\mathrm{eff}}{\rho}_{j}}}+B\sum_{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}^{\mathrm{eff}}{\rho}_{j}}-y{\rho}_{s})\)

• \(j\in \mathrm{forest}(s)\) sont les systèmes de normales distinctes \({n}_{j}\ne {n}_{s}\)

• \(j\in \mathrm{copla}(s)\) sont les systèmes de normales identiques \({n}_{j}={n}_{s}\)

La valeur initiale de la densité de dislocation est : \({\rho}_{s}(t=0)={\rho}_{s}^{0}\) (champ fourni par l’utilisateur. Bien vérifier s’il s’agit d’une densité globale (pour l’ensemble des 12 systèmes de glissement) ou bien d’une densité par système).

Comme \({a}_{ij}^{\mathrm{eff}}=C{(\rho )}^{2}.{a}_{ij}\) , \({\tau}_{s}^{\mathrm{forest}}(\rho )=\mu bC(\rho )\sqrt{\sum_{j=1,12}{a}_{\mathrm{sj}}{\rho}_{j}}\) et \(\dot{{\rho}_{s}}=\frac{\dot{{p}_{s}}}{b}(A\frac{\sum_{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}{\rho}_{j}}{\sum_{j=1,12}\sqrt{{a}_{\mathrm{sj}}{\rho}_{j}}}+BC(\rho )\sum_{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}{\rho}_{j}}-y{\rho}_{s})\)

Pour éviter tout problème de conditionnement de la matrice jacobienne, on choisit d’adimensionner le système en remplaçant la densité de dislocation par : \({\omega}_{s}={b}^{2}\ast {\rho}_{s}\) et résoudre en \({\omega}_{s}\) . Cela donne :

\(\dot{{\omega}_{s}}=\dot{{p}_{s}}{h}_{s}(\omega )\text{avec}{h}_{s}(\omega )=(A\frac{\sum_{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}{\omega}_{j}}{\sum_{j=1,12}\sqrt{{a}_{\mathrm{sj}}{\omega}_{j}}}+BC(\omega )\sum_{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}{\omega}_{j}}-\frac{y}{b}{\omega}_{s})\)

Remarque :physiquement, la densité de dislocation doit rester positive. Numériquement, pour que tous les termes en \(\sqrt{\rho}\) soient licites, il faut tester la positivité de \(\rho\) et soit redécouper le pas de temps si cette valeur devient négative (par «hésitation» de l’algorithme au cours des itérations) soit introduire une partie positive dans les équations d’évolution. Cette dernière solution est choisie pour l’intégration explicite comme pour l’intégration implicite.

L’ensemble des équations différentielles régissant l’écoulement viscoplastique du modèle DD_CFC s’écrit donc :

\(\dot{{\gamma}_{s}}=\dot{{p}_{s}}\frac{{\tau}_{s}}{∣{\tau}_{s}∣}\) (12 équations)

\(\dot{{p}_{s}}=\dot{{\gamma}_{0}}({(\frac{∣{\tau}_{s}∣}{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}})}^{n}-1)\text{si}∣{\tau}_{s}∣\ge {\tau}_{f}+{\tau}_{s}^{\mathrm{forest}},\text{sinon}\dot{{p}_{s}}=0\) (12 équations)

\(\dot{{\omega}_{s}}=\dot{{p}_{s}}{h}_{s}(\langle \omega \rangle )\) (12 équations) avec : \(\langle {\omega}_{i}\rangle =\lbrace \begin{array}{c}{\omega}_{i}\text{si}{\omega}_{i}>0\\ 0\text{sinon}\end{array}\)

\(C(\omega )=0.2+0.8\frac{\ln(\alpha \sqrt{\sum_{i=1,12}\langle {\omega}_{i}\rangle })}{\ln(\alpha b\sqrt{{\rho}_{\mathrm{ref}}})}\)

\({\tau}_{s}^{\mathrm{forest}}(\omega )=\mu C(\omega )\sqrt{\sum_{j=1,12}{a}_{\mathrm{sj}}\langle {\omega}_{j}\rangle }\)

\({h}_{s}(\omega )=(A\frac{\sum_{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}\langle {\omega}_{j}\rangle }{\sum_{j=1,12}\sqrt{{a}_{\mathrm{sj}}\langle {\omega}_{j}\rangle }}+BC(\omega )\sum_{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}\langle {\omega}_{j}\rangle }-\frac{y}{b}\langle {\omega}_{s}\rangle )\)

avec \({\omega}_{s}={b}^{2}\ast {\rho}_{s}\)

Dans le cas MONO_DD_CFC_IRRA, l’évolution des variables liées à l’irradiation est donnée par :

\({\dot{\rho}}_{s}^{\mathit{loops}}=-\xi \left(\sum_{j\in \mathit{copla}(s)}∣{\dot{\gamma}}_{j}∣\right)\left({\rho}_{s}^{\mathit{loops}}-{\rho}^{\mathit{sat}}\right)\) \({\dot{\varphi}}_{s}^{\mathit{voids}}=-\zeta \left(\sum_{j\in \mathit{copla}(s)}∣{\dot{\gamma}}_{j}∣\right)\left({\varphi}_{s}^{\mathit{voids}}-{\varphi}^{\mathit{sat}}\right)\)

l’expression de \({\tau}_{s}^{\mathit{forest}}\) devient : \({\tau}_{s}^{\mathit{forest}}=\mu \sqrt{\sum_{j=1,12}{a}_{\mathit{sj}}^{\mathit{eff}}{\omega}_{j}+{\alpha}^{\mathit{loops}}{\varphi}^{\mathit{loops}}{\rho}_{s}^{\mathit{loops}}{b}^{2}+{\alpha}^{\mathit{voids}}{\varphi}_{s}^{\mathit{voids}}{\rho}^{\mathit{voids}}{b}^{2}}\)

ce qui, après multiplication par \({b}^{\mathrm{2 }}\) : \({\omega}_{s}^{\mathit{loops}}={b}^{2}\ast {\rho}_{s}^{\mathit{loops}}\) donne:

\({\dot{\omega}}_{s}^{\mathit{loops}}=-\xi \left(\sum_{j\in \mathit{copla}(s)}∣{\dot{\gamma}}_{j}∣\right)\left({\omega}_{s}^{\mathit{loops}}-{\omega}^{\mathit{sat}}\right)\) et \({\tau}_{s}^{\mathit{forest}}=\mu \sqrt{\sum_{j=1,12}{a}_{\mathit{sj}}^{\mathit{eff}}{\omega}_{j}+{\alpha}^{\mathit{loops}}{\varphi}^{\mathit{loops}}{\omega}_{s}^{\mathit{loops}}+{\alpha}^{\mathit{voids}}{\varphi}_{s}^{\mathit{voids}}{\omega}^{\mathit{voids}}}\)

on peut remarquer au passage que les équations d’évolution de \({\omega}_{s}^{\mathit{loop}}\) et \({\varphi}_{s}^{\mathit{voids}}\) s’intègrent analytiquement:

\({\omega}_{s}^{\mathit{loops}}(t+\Delta t)={\omega}^{\mathit{sat}}+\left({\omega}_{s}^{\mathit{loops}}(t)-{\omega}^{\mathit{sat}}\right)\exp\left(-\xi \sum_{j\in \mathit{copla}(s)}∣{\Delta \gamma }_{j}∣\right)\)

\({\varphi}_{s}^{\mathit{voids}}(t+\Delta t)={\varphi}^{\mathit{sat}}+\left({\varphi}_{s}^{\mathit{voids}}(t)-{\varphi}^{\mathit{sat}}\right)\exp\left(-\zeta \sum_{j\in \mathit{copla}(s)}∣{\Delta \gamma }_{j}∣\right)\)

Dans le cas DD_FAT, l’évolution de la densité de dislocation est donnée par :

\(\dot{{\rho}_{s}}=\frac{\dot{{p}_{s}}}{b}(\frac{1}{d}+\frac{\sqrt{\sum_{u\mathrm{\ne }s}{\rho}_{u}}}{K}-{g}_{\mathit{c0}}{\rho}_{s})\)

En adimensionnant comme ci-dessus, c’est à dire en posant :

\({\omega}_{s}={b}^{2}\times {\rho}_{s}\)

on obtient :

\(\dot{{\omega}_{s}}=\dot{{p}_{s}}{h}_{s}(\omega )\) avec \({h}_{s}(\omega )=(\frac{b}{d}+\frac{\sqrt{\sum_{u\mathrm{\ne }s}{\omega}_{u}}}{K}-{g}_{\mathit{c0}}\frac{{\omega}_{s}}{b})\)

D’où en introduisant les parties positives :

\(\dot{{\gamma}_{s}}=\dot{{p}_{s}}\frac{{\tau}_{s}}{∣{\tau}_{s}∣}\) (12 équations)

\(\dot{{p}_{s}}=\dot{{\gamma}_{0}}({(\frac{∣{\tau}_{s}∣}{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}})}^{n}-1)\text{si}∣{\tau}_{s}∣\ge {\tau}_{f}+{\tau}_{s}^{\mathrm{forest}},\text{sinon}\dot{{p}_{s}}=0\) (12 équations)

\(\dot{{\omega}_{s}}=\dot{{p}_{s}}{h}_{s}(\langle \omega \rangle )\) (12 équations) avec : \(\langle {\omega}_{i}\rangle =\lbrace \begin{array}{c}{\omega}_{i}\text{si}{\omega}_{i}>0\\ 0\text{sinon}\end{array}\)

\(C(\omega )=1.\)

\({\tau}_{s}^{\mathrm{forest}}(\omega )=\mu \sqrt{\sum_{j=1,12}{a}_{\mathrm{sj}}\langle {\omega}_{j}\rangle }\)

\({h}_{s}(\omega )=(\frac{b}{d}+\frac{\sqrt{\sum_{u\ne s}\langle {\omega}_{u}\rangle }}{K}-{g}_{\mathrm{c0}}\frac{\langle {\omega}_{s}\rangle }{b})\)

avec \({\omega}_{s}={b}^{2}\ast {\rho}_{s}\)

La loi d’écrouissage des modèles DD_CFC est : \({\tau}_{s}^{\mathrm{forest}}=\mu b\sqrt{\sum_{j=1,{n}_{s}}{a}_{\mathrm{sj}}^{\mathrm{eff}}{\rho}_{j}}\)

et \({a}_{ij}^{\mathrm{eff}}=C{(\rho )}^{2}.{a}_{ij}\text{avec}C(\rho )=0.2+0.8\frac{\ln(\alpha b\sqrt{\sum_{i=1,12}{\rho}_{i}})}{\ln(\alpha b\sqrt{{\rho}_{\mathrm{ref}}})}\)

où \({a}_{ij}\) représente la matrice d’interaction entre systèmes de glissement, dont la forme particulière est fournie en annexe 4.

La loi d’écrouissage du modèle DD_FAT est : \({\tau}_{s}^{\mathrm{forest}}(\rho )=\mu b\sqrt{\sum_{j=1,12}{a}_{\mathrm{sj}}{\rho}_{j}}\) , ce qui correspond dans le modèle précédent à \(C(\rho )=1.\)

Lois issues de la Dynamique des Dislocations : MONO_DD_CC / MONO_DD_CC_IRRA

La formulation des lois DD_CC (avec ou sans irradiation) est construite à partir de calculs de dynamique des dislocations [ 17 ]. Elle s’applique aux matériaux à structure cristalline Cubique Centrée (\(\mathit{CC}\) )

Les équations ci-dessous sont écrites pour chaque système de glissement \(s\) .

L’écoulement visco-plastique est obtenu par moyenne harmonique de deux régimes (basse température, et haute température).

On choisit ici de faire apparaître les inconnues \(\rho\) (représentant l’ensemble des densités de dislocation \({\rho}^{s}\text{et}{\rho}_{\mathit{irr}}^{s}\) ) et \({\tau}_{\mathrm{eff}}^{s}\) (qui est une inconnue supplémentaire a priori. On verra que l’algorithme de résolution proposé permet de tout exprimer en fonction de \(\rho\) ).

\({\rho}_{\mathit{irr}}^{s}\) désigne la densité de dislocations induite par l’irradiation (cas de MONO_DD_CC_IRRA ), et n’est pas utilisée pour le comportement MONO_DD_CC afin de limiter le nombre de variables internes.

\({\tau}^{s}\) étant la cission résolue sur le système \(s\) :

• le glissement plastique est obtenue par : \(\frac{1}{{\dot{\gamma}}^{s}}=\frac{1}{{\dot{\gamma}}_{\mathrm{nuc}}^{s}}+\frac{1}{{\dot{\gamma}}_{\mathrm{prob}}^{s}}\)

• \({\dot{\gamma}}_{\mathrm{nuc}}^{S}={\rho}_{m}^{s}bH\cdot {l}^{s}(\rho ,{\tau}_{\mathrm{eff}}^{s})\cdot \exp(\frac{-\Delta G({\tau}_{\mathrm{eff}}^{s})}{{k}_{B}T})\mathrm{sgn}({\tau}_{s})\)

• \(\Delta G({\tau}_{\mathrm{eff}}^{s})=\Delta {G}_{0}(1-{(\frac{\langle {\tau}_{\mathrm{eff}}^{s}\rangle }{{\tau}_{0}})}^{0.5})\) (avec : \({\tau}_{\mathrm{eff}}^{s}<{\tau}_{0}\) et \(0<\Delta G<\Delta {G}_{0}\) )

• \({\tau}_{\mathrm{eff}}^{s}({\tau}^{s},\rho )=∣{\tau}^{s}∣-{\tau}_{c}\) \({\tau}_{c}={\tau}_{F}+\sqrt{({\tau}_{\mathit{LR}}^{s}{(\rho )}^{2}+{\tau}_{\text{LT}}^{s}{(\rho ,{\tau}_{\mathit{eff}}^{s})}^{2})}\)

• \({\tau}_{\mathrm{LR}}^{s}(\rho )=\mu b\sqrt{{a}_{\mathrm{AT}}^{s}(\rho ){\rho}^{s}}\)

• \({\tau}_{\text{LT}}^{s}(\rho ,{\tau}_{\mathrm{eff}}^{s})=\max(0;{\alpha}_{\mathrm{AT}}^{s}(\rho )\mu b(\frac{1}{{\lambda}^{s}(\rho ,{\tau}_{\mathrm{eff}}^{s})}-\frac{1}{2{\alpha}_{\mathrm{AT}}^{s}(\rho ){R}^{s}(\rho ,{\tau}_{\mathrm{eff}}^{s})+{l}_{c}(T)}))\)

• \({\rho}_{\mathrm{tot}}^{s}(\rho )=\sum_{j\ne s}{\rho}^{j}+{\rho}_{\mathrm{irr}}^{s}\)

• \({\alpha}_{\mathit{AT}}^{s}(\rho )=\sqrt{\sum_{j\ne s}{a}_{\mathit{AT}}^{\mathit{sj}}\frac{{\rho}^{j}}{{\rho}_{\mathit{tot}}^{s}}+{a}_{\mathit{irr}}\frac{{\rho}_{\mathit{irr}}^{s}}{{\rho}_{\mathit{tot}}^{s}}}\)

• \(\frac{1}{{\lambda}^{s}(\rho ,{\tau}_{\mathrm{eff}}^{s})+D}=\min(\sqrt{{\rho}_{\mathrm{tot}}^{s}(\rho )};(D+2{R}^{s}(\rho ,{\tau}_{\mathrm{eff}}^{s})){\rho}_{\mathrm{tot}}^{s}(\rho ))\)

• \({R}^{s}(\rho ,{\tau}_{\mathrm{eff}}^{s})=\frac{\mu b}{2{\tau}_{0}{(1-\frac{\Delta G({\tau}_{\mathrm{eff}}^{s})}{\Delta {G}_{0}})}^{2}}\)

• \({l}^{s}(\rho ,{\tau}_{\mathrm{eff}}^{s})=\max({\lambda}^{s}(\rho ,{\tau}_{\mathrm{eff}}^{s})-2{\alpha}_{\mathrm{AT}}^{s}(\rho ){R}^{s}(\rho ,{\tau}_{\mathrm{eff}}^{s});{l}_{c})\)

• \({\dot{\gamma}}_{\mathrm{prob}}^{s}=\dot{{\gamma}_{0}}{(\frac{∣{\tau}^{s}∣}{{\tau}_{c}})}^{n}\mathrm{sgn}({\tau}^{s})\)

Les paramètres sont : \({\rho}_{m}^{s},b,H,\Delta {G}_{0,}{\tau}_{0,}{\tau}_{F},\mu ,{a}_{\mathit{irr}},{l}_{c}(T),D,\dot{{\gamma}_{0,}}n\) et la matrice d’interaction \({\alpha}_{\mathrm{AT}}^{\mathrm{sr}}\)

Les lois d’évolution des densités de dislocations (\({\rho}^{s}`et :math:`{\rho}_{\mathrm{irr}}^{s}\)) sont :

• \(\dot{{\rho}^{s}}=\frac{∣\dot{{\gamma}^{s}}∣}{b}\left(\frac{1}{{d}_{\mathit{lath}}}+\frac{{c}_{\mathit{eff}}\sqrt{{a}_{\mathit{AT}}^{\mathit{ss}}{\rho}^{s}}}{{K}_{\mathit{self}}}+\frac{{\alpha}_{\mathit{AT}}^{s}{c}_{\mathit{eff}}{\lambda}^{s}{\rho}^{s}}{{K}_{f}}-{y}^{s}{\rho}^{s}\right)\)

• \(\dot{{\rho}_{\mathrm{irr}}^{s}}=-\xi {\rho}_{\mathrm{irr}}^{s}∣\dot{{\gamma}^{s}}∣\)

avec:

\(\frac{1}{{y}^{s}}=\frac{1}{{y}_{\mathrm{AT}}^{s}}+\frac{2\pi {\tau}_{\mathrm{eff}}^{s}(\rho )}{\mu b}\) et \({c}_{\mathit{eff}}^{s}=\left(1-\frac{{\tau}_{\mathit{eff}}^{s}}{{\tau}_{0}}\right)\) ,

avec les paramètres \({d}_{\mathit{lath}}{K}_{\mathit{self}}{K}_{f}\xi {y}_{\mathit{AT}}^{s}\)

Remarques :

• Physiquement, la densité de dislocation doit rester positive. Numériquement, pour que tous les termes en \(\sqrt{\rho}\) soient licites, il faut tester la positivité de \(\rho\) et soit redécouper le pas de temps si cette valeur devient négative (par «hésitation» de l’algorithme au cours des itérations) soit introduire une partie positive dans les équations d’évolution. Cette dernière solution est choisie pour l’intégration explicite comme pour l’intégration implicite.

• Dans le cas où l’élasticité n’est pas isotrope, la valeur de \(\mu\) utilisée ci-dessus est indépendante des coefficient d’élasticité, car obtenue par homogénéisation.

Systèmes de glissement et comportement global du monocristal

Un monocristal est composé de une ou plusieurs familles de systèmes de glissement, (cubique, octaédrique, basal, prismatique,… ), chaque famille comprenant un certain nombre de systèmes (12pour la famille octaédrique par exemple).

A chaque famille de système de glissement sont associés une loi d’écoulement, un type d’écrouissage cinématique et isotrope, et des valeurs des paramètres pour ces lois. En d’autres termes, on ne prévoit pas de faire varier les relations de comportement ou les coefficients au sein d’une même famille de systèmes de glissement. Par contre, d’une famille à l’autre, les lois de comportement peuvent changer, ainsi que les valeurs des paramètres.

Un système de glissement est déterminé par un tenseur d’orientation (ou tenseur de Schmid) construit à partir des définitions cristallographiques de:

• la direction de glissement du système \(s\) (définie par le vecteur unitaire \({m}_{s}\) )

• et de la normale au plan de glissement (définie par le vecteur unitaire \({n}_{s}\) ).

Petites déformations, configuration initiale

En petites déformations, le tenseur d’orientation est défini par :

\({\mu}_{ij}^{s}=\frac{1}{2}({m}_{i}\otimes {n}_{j}+{n}_{j}\otimes {m}_{i})\)

Du point de vue du comportement au point matériel, ce tenseur intervient pour le calcul de la cission résolue.

\({\tau}^{s}=\Sigma :{\mu}^{s}\)

et celui de la vitesse de déformation viscoplastique globale \({E}^{\mathrm{vp}}\) , définie à partir de la connaissance des vitesses de glissement \(\dot{{\gamma}^{s}}\) pour tous les systèmes de glissement:

\({\dot{E}}_{ij}^{\mathrm{vp}}=\sum_{s\in g}\dot{{\gamma}^{s}}{\mu}_{ij}^{s}\)

De plus, le monocristal peut être orienté par rapport aux axes globaux de définition des coordonnées. Cette orientation est définie pour chaque maille ou groupe de mailles (typiquement pour chaque grain) par la donnée de 3 angles nautiques ou 3 angles d’Euler.

Les composantes du tenseur d’orientation \({\mu}^{s}\) , définies dans le repère lié au monocristal, sont alors exprimées dans le repère global en utilisant ces angles nautiques. En petites déformations, elles sont fixes dans la configuration initiale.

De plus, au niveau d’un point de Gauss, on postule une relation d’élasticité sur les tenseurs globaux:

• déformation totale macroscopique \(E\)

• déformation viscoplastique macroscopique \({E}^{\mathrm{vp}}\)

• contrainte macroscopique : \(\Sigma\)

\(\Sigma =\Lambda (E-{E}^{\mathrm{th}}-{E}^{\mathrm{vp}})\)

où \(\Lambda\) représente l’opérateur d’élasticité (isotrope ou orthotrope)

Petites déformations, rotation du réseau cristallin

Une approximation de l’effet de la déformation sur le comportement du monocristal consiste à prendre en compte la rotation du réseau cristallin. Cette rotation est intéressante à deux titres :

• sa prise en compte effective dans les calculs est cohérente avec des déformations non infiniment petites (dans les calculs courants, il n’est pas rare d’atteindre 10%). Elle est un préliminaire à un modèle complet de grandes déformations sur le monocristal.

• C’est une information intéressante pour le post-traitement, car elle peut être confrontée à des résultats expérimentaux.

On décrit ici (suivant 9 , 10 , 11 ) la méthode de prise en compte de la rotation de réseau cristallin dans le cadre général des petites déformations. Les principales hypothèses physiques sont :

• connaissant le tenseur vorticité \(\stackrel{ˆ}{\dot{\omega}}\) ou taux de rotation, qui est la partie antisymétrique du tenseur \(L=\dot{F}{F}^{\text{-1}}\) ), on postule une décomposition additive en une partie dite rotation plastique \(\stackrel{ˆ}{\dot{{\omega}^{\mathrm{p}}}}\) (partie antisymétrique de la rotation plastique) et une partie «élastique» \(\stackrel{ˆ}{\dot{{\omega}^{e}}}\) , qui permettra de calculer la rotation du réseau cristallin ;

• on admet pour le moment que l’on peut calculer le tenseur taux de rotation, tout en effectuant l’intégration du comportement avec un tenseur des déformations linéarisées. Cette limitation tombe naturellement avec le modèle complet en grandes déformations, mais on suppose ici que l’on peut raisonner de façon découplée. On a donc, sous forme incrémentale :

\(\Delta \stackrel{ˆ}{\omega}=\frac{1}{2}(\Delta F{F}^{\text{-1}}-{(\Delta F{F}^{\text{-1}})}^{T})\) et \(\Delta \stackrel{ˆ}{{\omega}^{e}}=\Delta \stackrel{ˆ}{\omega}-\Delta \stackrel{ˆ}{{\omega}^{p}}\) qui sont des tenseurs antisymétriques.

Le calcul de \(\Delta \stackrel{ˆ}{{\omega}^{p}}\) est effectué en fin d’intégration : il est stocké et utilisé à l’itération suivante (un calcul purement implicite impliquerait de résoudre le système d’équations décrit au paragraphe 3 en ajoutant toutes les équations sur l’orientation du réseau, ce qui alourdirait beaucoup le calcul).

Une fois l’incrément de rotation \(\Delta \stackrel{ˆ}{{\omega}^{e}}\) connu, on extrait le vecteur «axial» [voir par exemple R5.03.40] : \(\Delta {\omega}^{e}\) tel que \(\Delta {\omega}^{e}\wedge U=\Delta \stackrel{ˆ}{{\omega}^{e}}\cdot U\) et on calcule l’angle de rotation et le vecteur qui correspond à l’axe de rotation :

\(\Delta \theta =\parallel \Delta {\omega}^{e}\parallel\) et \({n}_{\omega}=\frac{\Delta {\omega}^{e}}{\Delta \theta }\) .

La matrice d’incrément de rotation est alors obtenue par la formule d’Euler-Rodrigues :

\(\Delta Q={I}_{d}+\sin(\Delta \theta )\stackrel{ˆ}{{n}_{\omega}}+(1-\cos(\Delta \theta ))\stackrel{ˆ}{{n}_{\omega}}\cdot \stackrel{ˆ}{{n}_{\omega}}\) où \(\stackrel{ˆ}{{n}_{\omega}}\) est la matrice antisymétrique issue de \({\mathrm{n}}_{\omega}\)

Ce qui permet de mettre à jour la matrice de rotation du réseau : \({Q}_{+}=\Delta Q\cdot {Q}_{-}\)

Elle permet d’actualiser le tenseur d’orientation par : \(\begin{array}{}{n}^{s}={Q}_{+}{n}_{0}^{s}\\ {m}^{s}={Q}_{+}{m}_{0}^{s}\end{array}\) .

Le comportement est alors intégré avec cette définition du réseau cristallin.

Il reste alors à actualiser l’incrément de rotation plastique : \(\Delta \stackrel{ˆ}{{\omega}^{p}}=\frac{1}{2}{\sum}_{s}\Delta {\gamma}^{s}({m}^{s}\otimes {n}^{s}-{n}^{s}\otimes {m}^{s})\)

Les variables internes supplémentaires sont : \(Q,\Delta {\omega}^{p},\Delta {\omega}^{e},\theta\)

En pratique la rotation de réseau est prise en compte via le mot-clé ROTA_RESEAU de la commande DEFI_COMPOR [U4.43.06]. Elle n’est effective que pour l’intégration implicite du comportement MONOCRISTAL.

Grandes déformations

Propriétés attendues : le formalisme choisi doit être objectif, et doit conserver l’incompressibilité plastique ( \(\det({F}^{p})=1\) ). On utilise la décomposition multiplicative \(F={F}^{e}{F}^{p}\) avec la configuration intermédiaire relâchée isocline de Mandel.

La loi d”(hyper-)élasticité choisie utilise le second tenseur de Piola-Kirchhoff \(S={J}^{e}{F}^{e-1}\sigma {F}^{e-T}\) et la déformation de Green-Lagrange associée au gradient de déformation élastique :

\({E}_{\mathrm{GL}}^{e}=\frac{1}{2}({{F}^{e}}^{T}{F}^{e}–{I}_{d})\) \(S=\Lambda :{E}_{\mathrm{GL}}^{e}\) .

Sur chaque système de glissement \(s\) , la cission résolue est écrite avec le tenseur de Mandel \(M\) défini par \(M={J}_{e}{{F}^{e}}^{T}\sigma {{F}^{e}}^{\text{-T}}\) de la façon suivante : \({\tau}_{s}=M:{m}_{s}\otimes {n}_{s}\) (cf. 7 , 14 ).

En effet, la puissance des efforts intérieurs peut s’écrire : \(\frac{1}{\rho}\sigma :D=\frac{1}{\rho}\sigma :L=\frac{1}{\rho}\sigma :\dot{F}{F}^{\text{-1}}\)

et \(\dot{F}{F}^{\text{-1}}=(\dot{{F}^{\text{e}}}{F}^{p}+{F}^{\text{e}}\dot{{F}^{p}}){({F}^{p})}^{\text{-1}}{({F}^{e})}^{\text{-1}}={L}^{\text{e}}+{F}^{\text{e}}{L}^{p}{({F}^{e})}^{\text{-1}}\)

donc \(\frac{1}{\rho}\sigma :D=\frac{1}{\rho}\left[\sigma :{L}^{\text{e}}+\sigma :({F}^{\text{e}}{L}^{p}{({F}^{\text{e}})}^{\text{-1}})\right]=\frac{1}{{\rho}_{0}}S:\dot{{E}_{\mathrm{GL}}^{\text{e}}}+\frac{1}{{\rho}_{0}}M:{L}^{p}\)

La puissance dissipée plastiquement s’écrit : \(M:{L}^{p}\) et aussi : \(\sum_{s}\dot{{\gamma}_{s}}{\tau}_{s}\) , donc \(M:{L}^{p}=M:\sum_{s}\dot{{\gamma}_{s}}{\mu}_{s}=\sum_{s}\dot{{\gamma}_{s}}M:{\mu}_{s}\) , d’où \({\tau}_{s}=M:{\mu}_{s}\)

Pour chaque loi de comportement cristalline, les équations d’écoulement et d’évolution des variables internes s’écrivent classiquement sur chaque système de glissement :

\(\dot{{\gamma}_{s}}=g({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\) ou \(\dot{{\alpha}_{s}}=h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\)

Par contre la relation reliant la déformation plastique aux glissements des systèmes actifs s’écrit maintenant : \({\dot{F}}^{p}=(\sum_{s}\dot{{\gamma}_{s}}{m}_{s}\otimes {n}_{s}){F}^{p}\)

Comportement du polycristal homogénéisé

Dans le cas de polycristal homogénéisé, il faut définir chaque phase mono-cristalline par son orientation, sa proportion (fraction volumique) et le comportement associé. Une orientation globale du polycristal peut être définie. Il faut de plus définir une règle de localisation.

Le comportement monocristallin est construit comme précédemment à partir du comportement elasto‑visco-plastique précédent et de la donnée de familles de systèmes de glissement.

En plus du comportement monocristallin décrit précédemment, on ajoute une échelle de modélisation, qui représente l’assemblage des phases.

Au niveau d’un point de Gauss, on postule une relation d’élasticité sur les tenseurs globaux:

• déformation totale macroscopique \(E\)

• déformation viscoplastique macroscopique \({E}^{\mathrm{VP}}\)

• contrainte macroscopique : \(\Sigma\)

\(\Sigma =\Lambda (E-{E}^{\mathrm{th}}-{E}^{\mathrm{VP}})\)

\(\Lambda\) représente l’opérateur d’élasticité (isotrope ou orthotrope)

• De plus, connaissant l’ensemble des variables internes relatives aux systèmes de glissement de chaque phase, les paramètres de comportement de chaque phase, les orientations et fractions volumiques de chaque phase, et le type de méthode de localisation,

• pour chaque phase monocristalline (ou «grain»), définie par une orientation et une proportion \({f}_{g}\) , une relation de localisation des contraintes, de la forme générale (à exprimer dans le repère local de chaque phase)

\({\sigma}_{g}=L(\Sigma ,{E}^{\mathrm{vp}},{\varepsilon}_{g}^{\mathrm{vp}},{\beta}_{g})\)

et pour chaque système de glissementde chaque phase, des relations de comportement relatives à chaque système de glissement, similaires au cas du monocristal:

• Relation d’écoulement:

\(\dot{{\gamma}_{s}}=g({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\) , avec, \(\dot{{p}_{s}}=∣\dot{{\gamma}_{s}}∣\) et \(\dot{{p}_{s}}=f({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\)

• Évolution de l’écrouissage cinématique ou de la densité de dislocations: \(\dot{{\alpha}_{s}}=h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\)

• Évolution de l’écrouissage isotropedéfini par une fonction: \(R({p}_{s})\)

• Déformations viscoplastiques de la phase: \({\dot{\varepsilon}}_{ij}^{v{p}_{g}}=\sum_{s\in g}\dot{{\gamma}^{s}}{\mu}_{ij}^{s}\)

Il reste alors les équations d’homogénéisation: \({\dot{E}}^{\mathrm{vp}}=\sum_{g}{f}_{g}{\dot{\varepsilon}}_{g}^{\mathrm{vp}}\)

Deux relations de localisation de type \({\sigma}_{g}=L(\Sigma ,{E}^{\mathrm{vp}},{\varepsilon}_{g}^{\mathrm{vp}},{\beta}_{g})\) sont disponibles:

• La relation de Berveiller-Zaoui [bib5] établie sur la notion d’autocohérence. Cette relation est validée sous certaines conditions, à savoir: isotropie du matériau (des travaux récent mettent en oeuvre des extensions à l’élasticité cubique), comportement élastique homogène et chargement monotone:

\({\sigma}_{ij}^{g}={\Sigma}_{ij}+\alpha {\mu}^{\mathit{loca}}\left({E}_{ij}^{\mathit{vp}}-{\epsilon}_{ij}^{v{p}_{g}}\right)\) \(\frac{1}{\alpha}=1+\frac{3}{2}{\mu}^{\mathit{loca}}\frac{∥{E}_{ij}^{\mathit{vp}}∥}{{J}_{2}({\Sigma}_{ij})}\)

• La deuxième relation, inspirée de la précédente, et développée plus particulièrement pour des chargements cycliques [bib4] permet de donner une bonne description pour schématiser les interactions entre les grains:

• \({\sigma}_{ij}^{g}={\Sigma}_{ij}+{\mu}^{\mathit{loca}}\left({B}_{ij}-{\beta}_{ij}^{g}\right)\) \({B}_{ij}=\sum_{g}{f}_{g}{\beta}_{ij}^{g}\)

• \({\dot{\beta}}_{ij}^{g}={\dot{\varepsilon}}_{ij}^{v{p}_{g}}-D({\beta}_{ij}^{g}-\delta {\varepsilon}_{ij}^{v{p}_{g}})∥{\dot{\varepsilon}}_{ij}^{v{p}_{g}}∥\)

où \({\mu}^{\mathit{loca}}\) , \(D\) et \(\delta\) sont des paramètres caractéristiques du matériau et de la température.

Intégration locale et mise en œuvre numérique

Système d’équations à résoudre

Comportement de type MONOCRISTAL

Le comportement local du monocristal est défini, à un instant donné de la discrétisation en temps et au niveau d’un point d’intégration d’un élément fini, par la donnée:

• du tenseur de contraintes macroscopiques à l’instant précédent \(\Sigma ({t}_{i-1})={\Sigma}^{-}\) ,

• des variables internes à l’instant précédent, pour chaque système de glissement: \({\alpha}_{s}({t}_{i-1}),{\gamma}_{s}({t}_{i-1}),{p}_{s}({t}_{i-1})\) ,

• et du tenseur d’accroissement de déformation totale fournie par l’itération n de l’algorithme global de résolution \(\Delta E=\Delta {E}_{i}^{n}\) (avec les notations de [R5.03.01]).

L’intégration consiste à trouver:

• le tenseur de contraintes macroscopique \(\Sigma =\Sigma ({t}_{i})\)

• et les variables internes \({\alpha}_{s}={\alpha}_{s}({t}_{i})\) , \({\gamma}_{s}={\gamma}_{s}({t}_{i})\) , \({p}_{s}={p}_{s}({t}_{i})\)

vérifiant les équations de comportement dans chaque système de glissement (qui sont des relations mono dimensionnelles), et les relations de passage entre les tenseurs macroscopique et l’ensemble des directions de glissement. Notation: on écrit les équations sous forme discrétiséede façon:

• explicite, si les quantités notées \({A}^{\text{+/ -}}\) sont évaluées à l’instant \({t}_{i-1}\) : \({A}^{\text{+/ -}}={A}^{-}=A({t}_{i-1})\)

• implicite, si elles sont évaluées à l’instant \({t}_{i}\) : \({A}^{\text{+/ -}}={A}^{+}=A({t}_{i})\)

Les équations à intégrer peuvent se mettre sous la forme générale suivante:

Etant donnés, en un point de Gauss, les tenseurs:

\(\Delta E\) : variation de déformation à l’instant \({t}_{i}\) ,

\(E({t}_{i-1})={E}^{-}\) : déformation à l’instant \({t}_{i-1}\) ,

\(\Sigma ({t}_{i-1})={\Sigma}^{-}\) : contrainte macroscopique à l’instant \({t}_{i-1}\) ,

\({\alpha}_{s}({t}_{i-1}),{\gamma}_{s}({t}_{i-1}),{p}_{s}({t}_{i-1})\) : variables internespour chaque système de glissement à \({t}_{i-1}\) ,

Il faut trouver:

\(\Sigma =\Sigma ({t}_{i})\) : contrainte macroscopique à l’instant \({t}_{i}\) , dans le repère correspondant à l’orientation globale

\({\alpha}_{s}={\alpha}_{s}({t}_{i})\)

\({\gamma}_{s}={\gamma}_{s}({t}_{i})\)

\({p}_{s}={p}_{s}({t}_{i})\)

vérifiant:

\({D}^{-1}\Sigma =({D}_{-}^{-1}){\Sigma}^{-}+(\Delta E-\Delta {E}^{\text{th}}-\Delta {E}^{\text{vp}})\) , où \(\Lambda\) peut dépendre de la température, et peut correspondre à une élasticité orthotrope.

\(\Delta {E}^{\text{vp}}=\sum_{s}{\mu}_{s}\Delta {\gamma}_{s}\)

pour chaque système de glissement (de l’ensemble des familles de systèmes):

\({n}_{s}\) relations d’écoulement: soit en viscoplasticité \(\Delta {\gamma}_{s}=g({\tau}_{s}^{\text{+/ -}},{\alpha}_{s}^{\text{+/ -}},{\gamma}_{s}^{\text{+/ -}},{p}_{s}^{\text{+/ -}})\)

avec \(\Delta {p}_{s}=f({\tau}_{s}^{\text{+/ -}},{\alpha}_{s}^{\text{+/ -}},{\gamma}_{s}^{\text{+/ -}},{p}_{s}^{\text{+/ -}})\)

soit en plasticité \(F({\tau}_{s}^{\text{+/ -}},{\alpha}_{s}^{\text{+/ -}},{\gamma}_{s}^{\text{+/ -}},{p}_{s}^{\text{+/ -}})\le 0\) , \(F\Delta {p}_{s}=0\) ,

avec :math:`Delta {p}_{s}=mid Delta {gamma}_{s}mid `

\({n}_{s}\) équations d’évolutions de l’écrouissage cinématique: \(\Delta {\alpha}_{s}=h({\tau}_{s}^{\text{+/ -}},{\alpha}_{s}^{\text{+/ -}},{\gamma}_{s}^{\text{+/ -}},{p}_{s}^{\text{+/ -}})\)

\({n}_{s}\) équations d’évolution de l’écrouissage isotrope: \({R}_{s}({p}_{s}^{\text{+/ -}})\)

Ceci est résolu soit de façon explicite ( Runge_Kutta ), soit implicite ( Newton ).

Comportement de type POLYCRISTAL

Les relations de comportement discrétisées sont:

Etant donnés (en un point de Gauss) les tenseursglobaux:

• accroissement de déformation totale \(\Delta E\) ,

• déformation totale à l’instant précédent , \(E({t}_{i-1})={E}^{-}\)

• contrainte à l’instant précédent: \(\Sigma (t{}_{i-1}\text{})={\Sigma}^{-},\)

• l’ensemble des variables internes \({\alpha}_{s}^{-},{\gamma}_{s}^{-},{p}_{s}^{-}\) relatives aux systèmes de glissement de chaque phase,

• les paramètres de comportement de chaque phase, les orientations et fractions volumiques de chaque phase, et le type de méthode de localisation.

Il faut trouver \(\Sigma =\Sigma ({t}_{i}),{\alpha}_{s}=\alpha {}_{s}\text{}({t}_{i}),{\gamma}_{s}={\gamma}_{s}({t}_{i}),{p}_{s}={p}_{s}({t}_{i})\) vérifiant:

• au niveau du point de Gauss: \(\Sigma =\Lambda ({\Lambda}_{-}^{-1}){\Sigma}^{-}+\Lambda (\Delta E-\Delta {E}^{\text{th}}-\Delta {E}^{\text{vp}})\) , dans le repère global,

• pour chaque phase (ou «grain»), définie par une orientation et une proportion \({f}_{g}\) , une relation de localisation des contraintes, de la forme générale(à exprimer dans le repère local de chaque phase)

\({\sigma}_{g}=L(\Sigma ,{E}^{\text{vp}},{\varepsilon}_{g}^{\text{vp}},{\beta}_{g})\)

• et pour chaque système de glissementde chaque phase:

• \(\Delta {\varepsilon}_{{g}^{\text{vp}}}=\sum_{s}{m}_{s}\Delta {\gamma}_{s}\)

• \({n}_{s}\) équations: \({\tau}^{s}={\Sigma}_{ij}:{\mu}_{ij}^{s}\)

• \({n}_{s}\) relations d’écoulement: \(\Delta {\gamma}_{s}=g({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\) avec \(\Delta {p}_{s}=\mid \Delta {\gamma}_{s}\mid\)

• \({n}_{s}\) évolutions de l’écrouissage: \(\Delta {\alpha}_{s}=h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\)

• \(F({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\le 0,F\cdot \Delta {p}_{s}=0\) , (en plasticité indépendante du temps)

Il reste alors les équations d’homogénéisation: \(\Delta {E}^{\text{vp}}=\sum_{g}{f}_{g}\Delta {\varepsilon}_{g}^{\text{vp}}\)

Les comportements viscoplastiques relatifs à chaque système de glissement sont identiques au cas de la microstructure.

Dans la version actuelle de Code_Aster , ces relations de comportement sont intégrés uniquement de façon explicite.

Résolution implicite du comportement MONOCRISTAL

Résolution implicite des lois phénoménologiques

Il faut résoudre un système de la forme générale suivante:

\(R(Y)=\lbrace \begin{array}{}s(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ e(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ {n}_{s}\left\lbrace \begin{array}{}a(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ g(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ p(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\end{array}\right\rbrace \end{array}\) \(=\lbrace \begin{array}{}{\Lambda}^{\text{-1}}\Sigma -({\Lambda}_{-}^{\text{- 1}}){\Sigma}^{-}-(\Delta E-\Delta {E}^{\text{th}}-\Delta {E}^{\text{vp}})\\ \Delta {E}^{\text{vp}}-{\sum}_{s}{\mu}_{s}\Delta {\gamma}_{s}\\ {n}_{s}\left\lbrace \begin{array}{}\Delta {\alpha}_{s}-h({\tau}_{s}^{+},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ \Delta {\gamma}_{s}-g({\tau}_{s}^{+},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ {\mathrm{\Delta p}}_{s}-f({\tau}_{s}^{+},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\end{array}\right\rbrace \end{array}\) =0

De manière plus contractée, on pose:

\(R(\Delta y)=0=\left[\begin{array}{c}s(\Delta y)\\ e(\Delta y)\\ a(\Delta y)\\ g(\Delta y)\\ p(\Delta y)\end{array}\right]\) avec \(\Delta y=(\begin{array}{c}\Delta \Sigma \\ \Delta {E}^{\text{vp}}\\ \Delta {\alpha}_{s}\\ \Delta {\gamma}_{s}\\ \Delta {p}_{s}\end{array})\)

Pour résoudre ce système de \(6+6+3{n}_{s}\) équations non linéaires (en 3D), on utilise une méthode de Newton: on construit une suite de vecteurs solution de la façon suivante:

\({Y}_{k+1}={Y}_{k}-(\frac{\text{dR}}{{\text{dY}}_{k}}{)}^{-1}R({Y}_{k})\)

Il faut donc définir les valeurs initiales \({Y}_{0}\) , et calculer la matrice jacobienne du système: \(\frac{\text{dR}}{{\text{dY}}_{k}}\) (celle-ci est détaillée en annexe pour les comportements viscoplastiques décrits précédemment). Elle a la forme suivante:

\(J=\left[\begin{array}{ccccc}\frac{\partial s}{\partial \text{ΔΣ}}& \frac{\partial s}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial s}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial s}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial s}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial e}{\partial \text{ΔΣ}}& \frac{\partial e}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial e}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial e}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial e}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial a}{\partial \text{ΔΣ}}& \frac{\partial a}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial a}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial a}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial a}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial g}{\partial \text{ΔΣ}}& \frac{\partial g}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial g}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial g}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial g}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial p}{\partial \text{ΔΣ}}& \frac{\partial p}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial p}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial p}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial p}{\partial {\mathrm{\Delta p}}_{s}}\end{array}\right]\)

On peut toutefois simplifier ce système d’équations en réduisant sa taille, en vue d’améliorer les performances de l’intégration.

Dans l’expression des 6 composantes du tenseur des contraintes, \(\Delta {E}^{\text{vp}}\) peut être exprimée en fonction de \(\sum_{s}{\mu}_{s}\Delta {\gamma}_{s}\) donc une équation peut être éliminée du système global à résoudre.

Comme \(\Delta p=\mid \Delta \gamma \mid\) , l’équation en \(\Delta p\) peut être éliminée.

Dans le cas de MONO_CINE1 , \(\Delta \alpha\) s’exprime directement en fonction de \(\Delta \gamma\) . La taille du système à résoudre se réduit donc à \(6+{n}_{s}\) équations en \(\Delta \gamma\) , avec \({n}_{s}\) le nombre de systèmes de glissement.

Dans le cas de MONO_CINE2 , on calcule numériquement \(\Delta \alpha\) en fonction de \(\Delta \gamma\) , par une méthode de sécante.

L’équation se rapportant aux composantes des contraintes se réduit à:

\({R}_{1}^{i}({\Sigma}_{i},{\tau}_{s},\Delta {\gamma}_{s},\Delta {\alpha}_{s},{R}_{s}({p}_{r}))={\Lambda}^{-1}{\Sigma}_{i}-({\Lambda}_{-}^{-1}{\Sigma}_{{i}_{-}})-\Delta \varepsilon +\Delta {\varepsilon}^{\text{th}}+\sum_{s}{\mu}_{s}{g}_{s}({\tau}_{s},\Delta {\gamma}_{s},\Delta {\alpha}_{s},{R}_{s}({p}_{r}))=0\)

L’équation se rapportant au paramètre d’écrouissage \(\Delta \gamma\) s’écrit sous la forme:

\({R}_{2}^{s}({\Sigma}_{i},{\tau}_{s},\Delta {\gamma}_{s},\Delta {\alpha}_{s},{R}_{s}({p}_{r}))=\Delta {\gamma}_{s}-{g}_{s}({\tau}_{s},\Delta {\gamma}_{s},\Delta {\alpha}_{s},{R}_{s}({p}_{r}))=0\)

avec \(i=1à6\) , \(s=1\) à \({n}_{s}\) et \(r=1\) à \({n}_{s}\)

Dans le cas de MONO_VISC1 :

\({g}_{s}=\Delta t{\langle \frac{\mid {\tau}_{s}-c{\alpha}_{s}\mid -{R}_{s}({p}_{r})}{k}\rangle }^{n}\frac{{\tau}_{s}-{\mathrm{c\alpha }}_{s}}{\mid {\tau}_{s}-c{\alpha}_{s}\mid }\) avec \({\tau}_{s}=\Sigma :{\mu}_{s}\)

Dans le cas de MONO_VISC2 :

\({g}_{s}=\Delta t{\langle \frac{\mid {\tau}_{s}-c{\alpha}_{s}-a{\gamma}_{s}\mid -{R}_{s}({p}_{r})+\frac{d}{\mathrm{2c}}{(c{\alpha}_{s})}^{2}}{k}\rangle }^{n}\frac{{\tau}_{s}-c{\alpha}_{s}-a{\gamma}_{s}}{\mid {\tau}_{s}-c{\alpha}_{s}-a{\gamma}_{s}\mid }\)

La matrice jacobienne du système réduit s’écrit sous la forme:

\(J=\left[\begin{array}{cc}\frac{\partial {R}_{1}^{i}}{\partial {\Sigma}_{j}}& \frac{\partial {R}_{1}^{i}}{\partial \Delta {\gamma}_{s}}\\ \frac{\partial {R}_{2}^{s}}{\partial {\Sigma}_{i}}& \frac{\partial {R}_{2}^{s}}{\partial {\mathrm{\Delta \gamma }}_{r}}\end{array}\right]\) de dimension \((6+{n}_{s})\times (6+{n}_{s})\)

Les différents termes de la matrice jacobienne sont:

\(\frac{\partial {R}_{1}^{i}}{\partial {\Sigma}_{j}}={({\Lambda}^{-1})}_{ij}+\sum_{s}{({\mu}_{s})}_{i}\frac{\partial {g}_{s}}{\partial {\Sigma}_{j}}={({\Lambda}^{-1})}_{ij}+\sum_{s}{({\mu}_{s})}_{i}\frac{\partial {g}_{s}}{\partial {\tau}_{s}}\frac{\partial {\tau}_{s}}{\partial {\Sigma}_{j}}={({\Lambda}^{-1})}_{ij}+\sum_{s}{({\mu}_{s})}_{i}{({\mu}_{s})}_{j}\frac{\partial {g}_{s}}{\partial {\tau}_{s}}\)

\(\frac{\partial {R}_{2}^{s}}{\partial {\Sigma}_{i}}=\frac{\partial (-{g}_{s}({\tau}_{s},{\mathrm{\Delta \alpha }}_{s},{\mathrm{\Delta \gamma }}_{s},{R}_{s}({p}_{r})))}{\partial {\Sigma}_{i}}=\frac{\partial (-{g}_{s})}{\partial {\tau}_{s}}.\frac{\partial {\tau}_{s}}{\partial {\Sigma}_{i}}=-{\mu}_{s}^{i}.\frac{\partial {g}_{s}}{\partial {\tau}_{s}}\)

\(\frac{\partial {R}_{1}^{i}}{\partial \Delta {\gamma}_{s}}=\sum_{s}{\mu}_{s}^{i}.\frac{\partial {g}_{s}}{\partial \Delta {\gamma}_{s}}\)

\(\frac{\partial {R}_{2}^{s}}{\partial \Delta {\gamma}_{\mathrm{rs}}}={\delta}_{\mathrm{sr}}-\frac{\partial {g}_{s}}{\partial \Delta {\gamma}_{r}}\) avec \(\frac{\partial {g}_{s}}{\partial \Delta {\gamma}_{r}}=\frac{\partial {g}_{s}}{\partial \Delta {\alpha}_{s}}.\frac{\partial \Delta {\alpha}_{s}}{\partial \Delta {\gamma}_{s}}{\delta}_{\mathit{sr}}+\frac{\partial {g}_{s}}{\partial {R}_{s}}.\frac{\partial {R}_{s}}{\partial \Delta {\gamma}_{r}}\)

Dans les expressions intervenant ci-dessus , les termes \(\frac{\partial {g}_{s}}{\partial {\tau}_{s}},\frac{\partial {g}_{s}}{\partial \Delta {\alpha}_{s}},\frac{\partial \Delta {\alpha}_{s}}{\partial \Delta {\gamma}_{s}},\frac{\partial {g}_{s}}{\partial {R}_{s}},\frac{\partial {R}_{s}}{\partial \Delta {p}_{r}}\) dépendent les lois d’écoulement et d’écrouissage:

• Dans le cas de MONO_VISC1 :

\(\frac{\partial {g}_{s}}{\partial {\tau}_{s}}=\frac{n\Delta t}{{k}^{n}}{\langle \mid {\tau}_{s}-c{\alpha}_{s}\mid -{R}_{s}\left({p}_{r}\right)\rangle }^{n-1}={H}_{s}\left({\tau}_{s},\Delta {\alpha}_{s},\Delta {\gamma}_{s},{R}_{s}\left({p}_{r}\right)\right)\)

\(\frac{\partial {g}_{s}}{\partial \Delta {\alpha}_{s}}=-c.{H}_{s}({\tau}_{s},\Delta {\gamma}_{s},{R}_{s})\)

\(\frac{\partial {g}_{s}}{\partial {R}_{s}}=-{H}_{s}\left({\tau}_{s},\Delta {\gamma}_{s},{R}_{s}\right)\text{sgn}s\) avec \(\text{sgn}s=\frac{{\tau}_{s}-c{\alpha}_{s}}{\mid {\tau}_{s}-c{\alpha}_{s}\mid }\)

• Dans le cas de MONO_VISC2 :

\(\frac{\partial {g}_{s}}{\partial {\tau}_{s}}={H}_{s}\left({\tau}_{s},\Delta {\alpha}_{s},\Delta {\gamma}_{s},{R}_{s}\left({p}_{r}\right)\right)=\frac{n}{{k}^{n}}{\langle \mid {\tau}_{s}-c{\alpha}_{s}-a{\gamma}_{s}\mid -{R}_{s}\left({p}_{r}\right)+\frac{d}{\mathrm{2c}}{\left(c{\alpha}_{s}\right)}^{2}\rangle }^{n-1}\)

\(\frac{\partial {g}_{s}}{\partial \Delta {\alpha}_{s}}={H}_{s}({\tau}_{s},\Delta {\gamma}_{s},{R}_{s})(-c\text{sgn}s+d{\alpha}_{s}\Delta {\alpha}_{s})\)

\(\frac{\partial {g}_{s}}{\partial {R}_{s}}=-{H}_{s}\left({\tau}_{s},\Delta {\gamma}_{s},{R}_{s}\right)\text{sgn}s\) avec \(\text{sgn}s=\frac{{\tau}_{s}-c{\alpha}_{s}-a{\gamma}_{s}}{\mid {\tau}_{s}-c{\alpha}_{s}-a{\gamma}_{s}\mid }\)

• Dans le cas de MONO_CINE1:

\(\frac{\partial \Delta {\alpha}_{s}}{\partial \Delta {\gamma}_{r}}=\frac{1-d{\alpha}_{s}^{-}\text{sgn}s}{{\left(1+d\mid \Delta {\gamma}_{s}\mid \right)}^{2}}\)

• Dans le cas de MONO_CINE2:

\(\frac{\partial {\alpha}_{s}}{\partial \Delta {\gamma}_{r}}\) ; Cette dérivation nécessite la résolution d’une équation non-linéaire par la méthode de Newton

• Dans le cas de MONO_ISOT1:

\(\frac{\partial {R}_{s}}{\partial \Delta {\gamma}_{r}}=\frac{\partial \left[Q\sum_{k}{h}_{\text{sk}}(1-{e}^{-{\text{bp}}_{k}})\right]}{\partial \Delta {\gamma}_{r}}=\frac{\partial}{\partial {p}_{r}}\left[Q\sum_{k}{h}_{\text{sk}}(1-{e}^{-{\text{bp}}_{k}})\right]\frac{\partial {p}_{r}}{\partial \Delta {\gamma}_{r}}={\text{bQh}}_{\text{sr}}{e}^{-{\text{bp}}_{r}}.\text{sgn}r\)

avec \(\text{sgn}r=\frac{\partial {p}_{r}}{\partial \Delta {\gamma}_{r}}=\frac{\partial \mid \Delta {\gamma}_{r}\mid }{\partial \Delta {\gamma}_{r}}=\frac{\Delta {\gamma}_{r}}{\mid \Delta {\gamma}_{r}\mid }\)

• Dans le cas de MONO_ISOT2:

\(\begin{array}{}\frac{\partial {R}_{s}}{\partial \Delta {\gamma}_{r}}=\frac{\partial \left[{Q}_{1}\sum_{k}{h}_{\text{sk}}(1-{e}^{-{b}_{1}{p}_{k}})+{Q}_{2}(1-{e}^{-{b}_{2}{p}_{r}})\right]}{\partial \Delta {\gamma}_{r}}=\\ \frac{\partial}{\partial {p}_{r}}\left[{Q}_{1}\sum_{k}{h}_{\text{sk}}(1-{e}^{-{b}_{1}{p}_{k}})\right]\frac{\partial {p}_{r}}{\partial \Delta {\gamma}_{r}}+\frac{\partial}{\partial {p}_{r}}\left[{Q}_{2}(1-{e}^{-{b}_{2}{p}_{r}})\right]\frac{\partial {p}_{r}}{\partial \Delta {\gamma}_{r}}=\\ ({b}_{1}{Q}_{1}{h}_{\text{sr}}{e}^{-{b}_{1}{p}_{r}}+{b}_{2}{Q}_{2}{e}^{-{b}_{2}{p}_{r}}{\delta}_{\text{sr}})\text{sgn}r\end{array}\)

avec \(\text{sgn}r=\frac{\partial {p}_{r}}{\partial \Delta {\gamma}_{r}}=\frac{\partial \mid \Delta {\gamma}_{r}\mid }{\partial \Delta {\gamma}_{r}}=\frac{\Delta {\gamma}_{r}}{\mid \Delta {\gamma}_{r}\mid }\)

Résolution implicite – Intégration du modèle MONO_DD_KR

Le système à résoudre se réduit aux équations relatives aux composantes des contraintes

\({R}_{1}^{i}({\Sigma}_{i},{\tau}_{s},{\alpha}_{s})={({\Lambda}^{-1}\Sigma )}_{i}-{({\Lambda}_{-}^{\text{-1}}{\Sigma}_{-})}_{i}-\Delta {\varepsilon}_{i}+\Delta {\varepsilon}_{i}^{\text{th}}+\sum_{s}{\mu}_{i}^{s}{g}_{s}({\tau}_{s},{\alpha}_{s})=0\)

et celles relatives au paramètre d’écrouissage:

\({R}_{2}^{s}({\sum}_{i},{\tau}_{s},{\alpha}_{s})=\Delta {\alpha}_{s}-\mid {g}_{s}({\tau}_{s},{\alpha}_{s})\mid .h({\alpha}_{s})=0\)

avec \(i\) = 1 à 6 , \(s\) = 1 à \(\mathrm{ns}\) et \(r\) =1 à \(\mathrm{ns}\)

Ceci permet d’obtenir un système d’équations formulé de façon entièrement implicite, où \({R}_{1}^{i}\) et \({R}_{2}^{s}\) ne dépendent que de \({\Sigma}_{i},{\alpha}_{s}\) (car \({\tau}_{s}\) se calcule explicitement à partir de \({\Sigma}_{i}\)

\({g}_{s}({\tau}_{s},{\alpha}_{s})={\dot{\gamma}}_{0}\exp(\frac{-\Delta G({\tau}_{\text{eff}}^{s})}{{k}_{B}T}).\frac{{\tau}_{s}}{\mid {\tau}_{s}\mid }.\Delta t\)

\(\Delta G({\tau}_{\text{eff}}^{s})=\Delta {G}_{0}{(1-{(\frac{\langle {\tau}_{\text{eff}}^{s}\rangle }{{\tau}_{R}})}^{p})}^{q}\)

\(h({\alpha}_{s})=(\frac{b}{d}+\frac{\sqrt{\sum_{u\ne s}{\alpha}^{u}}}{K}-\frac{{g}_{c}{\alpha}^{s}}{b})\)

\({\tau}_{\text{eff}}^{s}=\langle \mid {\tau}^{s}\mid -{\tau}_{0}-{\tau}_{\mu}^{s}\rangle\) \(\langle \rangle\) : partie positive

\({\tau}_{\mu}^{s}=\frac{{(\mu )}^{2}\sum_{u}{a}^{\text{su}}{\alpha}^{u}}{\mid {\tau}^{s}\mid -{\tau}_{0}}\) avec \({\alpha}^{u}={\rho}^{u}{b}^{2}\)

Remarque: pour chaque système de glissement \(s\) où les conditions \({\tau}_{\text{eff}}^{s}>0\mid {\tau}^{s}\mid >{\tau}_{0}\) ne sont pas vérifiées, on obtient immédiatement: \({g}_{s}({\tau}_{s},{\alpha}_{s})=0,\Delta {\alpha}_{s}=0\)

La matrice jacobienne du système s’écrit sous la forme:

\(J=\left[\begin{array}{cc}\frac{\partial {R}_{1}^{i}}{\partial {\Sigma}_{j}}& \frac{\partial {R}_{1}^{i}}{\partial \Delta {\alpha}_{s}}\\ \frac{\partial {R}_{2}^{s}}{\partial {\Sigma}_{i}}& \frac{\partial {R}_{2}^{s}}{\partial \Delta {\alpha}_{r}}\end{array}\right]\) de dimension \((6+{n}_{s})\times (6+{n}_{s})\)

• \(\frac{\partial {R}_{1}^{i}}{\partial {\Sigma}_{j}}={({\Lambda}^{-1})}_{ij}+\sum_{s}{({\mu}_{s})}_{i}\frac{\partial {g}_{s}}{\partial {\Sigma}_{j}}={({\Lambda}^{-1})}_{ij}+\sum_{s}{({\mu}_{s})}_{i}\frac{\partial {g}_{s}}{\partial {\tau}_{s}}\frac{\partial {\tau}_{s}}{\partial {\Sigma}_{j}}={({\Lambda}^{-1})}_{ij}+\sum_{s}{({\mu}_{s})}_{i}{({\mu}_{s})}_{j}\frac{\partial {g}_{s}}{\partial {\tau}_{s}}\)

=> Il faut calculer \(\frac{\partial {g}^{s}}{\partial {\tau}_{s}}=\text{sgn}({\tau}_{s}).\frac{\partial \mid {g}^{s}\mid }{\partial {\tau}_{s}}\)

• \(\frac{\partial {R}_{2}^{s}}{\partial {\Sigma}_{i}}=\frac{\partial (-\mid {g}_{s}\mid .{h}_{s})}{\partial {\tau}_{s}}.\frac{\partial {\tau}_{s}}{\partial {\Sigma}_{i}}=-{\mu}_{s}^{i}.(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau}_{s}}.{h}_{s}+\mid {g}_{s}\mid .\frac{\partial {h}_{s}}{\partial {\tau}_{s}})\)

=> Le terme \(\frac{\partial {h}_{s}}{\partial {\tau}_{s}}\) étant nul, il reste à calculer \(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau}_{s}}\) .

• \(\frac{\partial {R}_{1}^{i}}{\partial \Delta {\alpha}_{s}}=\sum_{s}{\mu}_{s}^{i}\frac{{\tau}_{s}}{\mid {\tau}_{s}\mid }.\frac{\partial \mid {g}_{s}\mid }{\partial {\alpha}_{s}}\)

• \(\frac{\partial {R}_{2}^{s}}{\partial \Delta {\alpha}_{r}}={\delta}_{\text{sr}}-\frac{\partial \mid {g}_{s}\mid }{\partial \Delta {\gamma}_{r}}={\delta}_{\text{sr}}-(\frac{\partial \mid {g}_{s}\mid }{\partial \Delta {\alpha}_{r}}{h}_{s}+\mid {g}_{s}\mid \frac{\partial {h}_{s}}{\partial \Delta {\alpha}_{r}})\)

=> Il faut aussi calculer \(\frac{\partial \mid {g}_{s}\mid }{\partial {\alpha}_{r}}\) et \(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau}_{s}}\) ;

\({h}_{s}=\frac{\Delta {\alpha}_{s}}{\mid \Delta {\gamma}_{s}\mid }\) et \(\frac{\partial {\tau}_{\mu}^{s}}{\partial {\tau}_{s}}=\frac{-{\tau}_{\mu}^{s}}{\mid {\tau}^{s}\mid -{\tau}_{0}}.\text{sgn}s\)

Calcul de \(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau}_{r}}\) :

\(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau}_{s}}=-\frac{\Delta t.{\dot{\gamma}}_{0}}{{k}_{B}T}\exp(-\frac{\Delta G({\tau}_{\text{eff}}^{s})}{{k}_{B}T}).\frac{\Delta G({\tau}_{\text{eff}}^{s})}{\Delta {\tau}_{s}}=\frac{-\mid {g}_{s}({\tau}_{s},{\alpha}_{s})\mid }{{k}_{B}T}.\frac{\partial \Delta G({\tau}_{\text{eff}}^{s})}{\partial {\tau}_{s}}\)

où

\(\frac{\partial \Delta G({\tau}_{\text{eff}}^{s})}{\partial {\tau}_{s}}=-\Delta {G}_{0}{(1-{(\frac{{\tau}_{\text{eff}}^{s}}{{\tau}_{R}})}^{p})}^{q-1}.\frac{q.p}{{\tau}_{R}}.{(\frac{{\tau}_{\text{eff}}^{s}}{{\tau}_{R}})}^{p-1}.\frac{\partial {\tau}_{\text{eff}}^{s}}{\partial {\tau}_{s}}\)

et

\(\frac{\partial {\tau}_{\text{eff}}^{s}}{\partial {\tau}_{s}}=\frac{{\tau}_{s}}{\mid {\tau}_{s}\mid }-\frac{\partial {\tau}_{\mu}^{s}}{\partial {\tau}_{s}}=\text{sgn}({\tau}_{s})+{\mu}^{2}(\sum{a}^{\text{su}}{\alpha}^{u}).\frac{\text{sgn}({\tau}_{s})}{{(\mid {\tau}_{s}\mid -{\tau}_{0})}^{2}}=\text{sgn}({\tau}_{s})+{\tau}_{\mu}^{s}.\frac{\text{sgn}({\tau}_{s})}{\mid {\tau}_{s}\mid -{\tau}_{0}}\)

Calcul de \(\frac{\partial \mid {g}_{s}\mid }{\partial {\alpha}_{r}}\) :

\(\frac{\partial \mid {g}_{s}\mid }{\partial {\alpha}_{r}}=-\frac{\Delta t.{\dot{\gamma}}_{0}}{{k}_{B}T}\exp(-\frac{\Delta G({\gamma}_{\text{eff}}^{s})}{{k}_{B}T}).\frac{\Delta G({\tau}_{\text{eff}}^{s})}{\Delta {\alpha}_{r}}=\frac{-\mid {g}_{s}({\tau}_{s},{\alpha}_{s})\mid }{{k}_{B}T}.\frac{\partial \Delta G({\tau}_{\text{eff}}^{s})}{\partial {\alpha}_{r}}\)

où

\(\frac{\partial \Delta G({\tau}_{\text{eff}}^{s})}{\partial {\alpha}_{r}}=-\frac{\text{qp}}{{\tau}_{R}}\Delta {G}_{0}{(1-{(\frac{{\tau}_{\text{eff}}}{{\tau}_{R}})}^{p})}^{q-1}.{(\frac{{\tau}_{\text{eff}}}{{\tau}_{R}})}^{p-1}.\frac{\partial {\tau}_{\text{eff}}^{s}}{\partial {\alpha}_{r}}\)

et

\(\frac{\partial {\tau}_{\text{eff}}^{s}}{\partial {\alpha}_{r}}=-\frac{\partial {\tau}_{\mu}^{s}}{\partial {\alpha}_{r}}=-\frac{{\mu}^{2}}{\mid {\tau}_{s}\mid -{\tau}_{0}}\frac{\partial (\sum_{u}{a}^{\text{su}}{\alpha}^{u})}{\partial {\alpha}_{r}}=-\frac{{\mu}^{2}}{\mid {\tau}_{s}\mid -{\tau}_{0}}{a}^{\text{sr}}\)

Calcul de \(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau}_{s}}\) :

\(\frac{\partial {h}^{s}}{\partial {\alpha}_{r}}=\frac{1}{K}\frac{\partial}{\partial {\alpha}_{r}}\sqrt{\sum_{u\ne s}{\alpha}^{u}}-\frac{{g}_{c}}{b}{\delta}_{\text{sr}}=\frac{1}{\mathrm{2K}}\frac{1-{\delta}_{\text{rs}}}{\sqrt{\sum_{u\ne s}{\alpha}^{u}}}-\frac{{g}_{c}}{b}{\delta}_{\text{sr}}\)

Résolution implicite – modèles MONO_DD_CFC(IRRA) / MONO_DD_FAT

Expression du résidu Connaissant la solution à l’instant \({t}_{i-1}\) , notée \({\sigma}^{-}=\sigma ({t}_{i-1}),{\rho}_{s}^{-}=\rho {}_{s}\text{}({t}_{i-1}),{\gamma}_{s}^{-}={\gamma}_{s}({t}_{i-1}),{p}_{s}^{-}={p}_{s}({t}_{i-1})\) ,

il s’agit de trouver, à l’instant \({t}_{i}\) , les quantités \(\sigma =\sigma ({t}_{i}),{\rho}_{s}=\rho {}_{s}\text{}({t}_{i}),{\gamma}_{s}={\gamma}_{s}({t}_{i}),{p}_{s}={p}_{s}({t}_{i})\) vérifiant:

\({({\Lambda}^{-1}\sigma )}_{i}-{({({\Lambda}^{-})}^{\text{-1}}{\sigma}^{-})}_{i}=\Delta {\varepsilon}_{i}-\Delta {\varepsilon}_{i}^{\text{th}}-\sum_{s}({({\mu}_{s})}_{i}\Delta {\gamma}_{s}(\sigma ,\omega ))=0\)

La cission résolue par système de glissement est : \({\tau}_{s}=\sigma :{\mu}_{s}=({\sigma}^{-}+\Delta \sigma ):{\mu}_{s}\)

La loi d’écoulement après discrétisation implicite devient :

\({\mathrm{\Delta \gamma }}_{s}={\dot{\gamma}}_{0}\mathrm{\Delta t}({(\frac{\mid {\tau}_{s}\mid }{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}})}^{n}-1).\frac{{\tau}_{s}}{\mid {\tau}_{s}\mid }=\Delta {p}_{s}\frac{{\tau}_{s}}{\mid {\tau}_{s}\mid }\)

\(\Delta {p}_{s}(\sigma ,\omega )=\dot{{\gamma}_{0}}\Delta t({(\frac{\mid {\tau}_{s}\mid }{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}(\omega )})}^{n}-1)\) si \(∣{\tau}_{s}∣\ge {\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}\) ; sinon \(\Delta {p}_{s}=0\)

La loi d’écrouissage est : \({\tau}_{s}^{\mathrm{forest}}(\omega )={\tau}_{s}^{\mathrm{forest}}(\langle {\omega}^{-}+\Delta \omega \rangle )=\mu C(\omega )\sqrt{{\sum}_{j=1,12}\langle {a}_{\mathrm{sj}}{\omega}_{j}\rangle }\)

En présence d’irradiation, la loi d’écrouissage devient:

\({\tau}_{s}^{\mathit{forest}}=\mu \sqrt{\sum_{j=1,12}C{(\omega )}^{2}{a}_{\mathit{sj}}{\omega}_{j}+{\alpha}^{\mathit{loops}}{\varphi}^{\mathit{loops}}{\omega}_{s}^{\mathit{loops}}+{\alpha}^{\mathit{voids}}{\varphi}_{s}^{\mathit{voids}}{\omega}^{\mathit{voids}}}\)

avec \(C(\omega )=0.2+0.8\frac{\ln\left(\alpha \sqrt{\sum_{i=1,12}\langle {\omega}_{i}\rangle }\right)}{\ln\left(\alpha b\sqrt{{\rho}_{\mathit{ref}}}\right)}\) ou, dans le cas DD_FAT, \(C(\omega )=1.\)

L’évolution de la densité de dislocationest donnée (toujours sous forme implicite) par :

\(\Delta {\omega}_{s}=\Delta {p}_{s}{h}_{s}(\omega )=\Delta {p}_{s}{h}_{s}({\omega}^{-}+\Delta \omega )\) avec :

\({h}_{s}(\omega )=(A\frac{\sum_{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}\langle {\omega}_{j}\rangle }{\sum_{j=1,12}\sqrt{{a}_{\mathrm{sj}}\langle {\omega}_{j}\rangle }}+BC(\omega )\sum_{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}\langle {\omega}_{j}\rangle }-\frac{y}{b}{\omega}_{s})\)

ou, dans le cas DD_FAT, \({h}_{s}(\omega )=(\frac{b}{d}+\frac{\sqrt{\sum_{u\ne s}\langle {\omega}_{u}\rangle }}{K}-{g}_{\mathrm{c0}}\frac{\langle {\omega}_{s}\rangle }{b})\)

L’évolution des densités de dislocations liées à l’irradiation est de la forme suivante:

\({\omega}_{s}^{\mathit{loops}}(t+\Delta t)={\omega}^{\mathit{sat}}+\left({\omega}_{s}^{\mathit{loops}}(t)-{\omega}^{\mathit{sat}}\right)\exp\left(-\xi \sum_{j\in \mathit{copla}(s)}∣{\Delta \gamma }_{j}∣\right)\)

\({\varphi}_{s}^{\mathit{voids}}(t+\Delta t)={\varphi}^{\mathit{sat}}+\left({\varphi}_{s}^{\mathit{voids}}(t)-{\varphi}^{\mathit{sat}}\right)\exp\left(-\zeta \sum_{j\in \mathit{copla}(s)}∣{\Delta \gamma }_{j}∣\right)\)

En utilisant la notation vectorielle de Kelvin [chapitre 5 de 15 ) pour les tenseurs symétriques (\(\sigma\) , \({\mu}_{s}\) ), et en prenant en compte les éventuelles variations des coefficients de l’élasticité avec la température, on aboutit au système à résoudre, comportant 18 équations ; \(i=1,6\) , \(s=1,12\) , les 18 inconnues étant les 6 composantes de \(\sigma\) et les 12 valeurs \(\omega\) .

Système à résoudre en \(\sigma\) (6 composantes) et \(\Delta \omega\) (12 composantes) :

\({R}_{i}^{(1)}(\sigma ,\omega )={({\Lambda}^{-1}\sigma )}_{i}-{({({\Lambda}^{-})}^{\text{-1}}{\sigma}^{-})}_{i}-\Delta {\varepsilon}_{i}+\Delta {\varepsilon}_{i}^{\text{th}}+\sum_{s}({({\mu}_{s})}_{i}\Delta {p}_{s}(\sigma ,\omega ).\frac{{\tau}_{s}}{\mid {\tau}_{s}\mid })=0\)

\({R}_{s}^{(2)}(\sigma ,\omega )=\Delta {\omega}_{s}-\Delta {p}_{s}{h}_{s}(\omega )=0\)

avec : \({\tau}_{s}=\sigma :{\mu}_{s}=({\sigma}^{-}+\Delta \sigma ):{\mu}_{s}\)

\(\Delta {p}_{s}(\sigma ,\omega )=\dot{{\gamma}_{0}}\Delta t({(\frac{\mid {\tau}_{s}\mid }{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}(\omega )})}^{n}-1)\) si \(∣{\tau}_{s}∣\ge {\tau}_{0}+{\tau}_{s}^{f}\) ; sinon \(\Delta {p}_{s}=0\)

et si \(\Delta {p}_{s}>0\) :

\({\tau}_{s}^{\mathrm{forest}}(\omega )={\tau}_{s}^{\mathrm{forest}}({\omega}^{-}+\Delta \omega )=\mu C(\omega )\sqrt{{\sum}_{j=1,12}{a}_{\mathrm{sj}}\langle {\omega}_{j}\rangle }\)

ou \({\tau}_{s}^{\mathit{forest}}=\mu \sqrt{\sum_{j=1,12}C{(\omega )}^{2}{a}_{\mathit{sj}}{\omega}_{j}+{\alpha}^{\mathit{loops}}{\varphi}^{\mathit{loops}}{\omega}_{s}^{\mathit{loops}}+{\alpha}^{\mathit{voids}}{\varphi}_{s}^{\mathit{voids}}{\omega}^{\mathit{voids}}}\)

\(C(\omega )=0.2+0.8\frac{\ln(\alpha \sqrt{\sum_{i=1,12}\langle {\omega}_{i}\rangle })}{\ln(\alpha b\sqrt{{\rho}_{\mathrm{ref}}})}\)

\({h}_{s}(\omega )=(A\frac{\sum_{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}\langle {\omega}_{j}\rangle }{\sum_{j=1,12}\sqrt{{a}_{\mathrm{sj}}\langle {\omega}_{j}\rangle }}+BC(\omega )\sum_{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}\langle {\omega}_{j}\rangle }-\frac{y}{b}\langle {\omega}_{s}\rangle )\)

Dans le cas DD_FAT, ce système devient :

\({R}_{i}^{(1)}(\sigma ,\omega )={({\Lambda}^{-1}\sigma )}_{i}-{({({\Lambda}^{-})}^{\text{-1}}{\sigma}^{-})}_{i}-\Delta {\varepsilon}_{i}+\Delta {\varepsilon}_{i}^{\text{th}}+\sum_{s}({({\mu}_{s})}_{i}\Delta {p}_{s}(\sigma ,\omega ).\frac{{\tau}_{s}}{\mid {\tau}_{s}\mid })=0\)

\({R}_{s}^{(2)}(\sigma ,\omega )=\Delta {\omega}_{s}-\Delta {p}_{s}{h}_{s}(\omega )=0\)

avec : \({\tau}_{s}=\sigma :{\mu}_{s}=({\sigma}^{-}+\Delta \sigma ):{\mu}_{s}\)

\(\Delta {p}_{s}(\sigma ,\omega )=\dot{{\gamma}_{0}}\Delta t({(\frac{\mid {\tau}_{s}\mid }{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}(\omega )})}^{n}-1)\) si \(∣{\tau}_{s}∣\ge {\tau}_{0}+{\tau}_{s}^{f}\) ; sinon \(\Delta {p}_{s}=0\)

et si \(\Delta {p}_{s}>0\) :

\({\tau}_{s}^{\mathrm{forest}}(\omega )={\tau}_{s}^{\mathrm{forest}}({\omega}^{-}+\Delta \omega )=\mu \sqrt{{\sum}_{j=1,12}{a}_{\mathrm{sj}}\langle {\omega}_{j}\rangle }\)

\({h}_{s}(\omega )=(\frac{b}{d}+\frac{\sqrt{\sum_{u\ne s}\langle {\omega}_{u}\rangle }}{K}-{g}_{\mathrm{c0}}\frac{\langle {\omega}_{s}\rangle }{b})\)

La matrice jacobienne du système peut être alors calculée pour l’intégration par la méthode de Newton (R5.03.14).

Matrice jacobienne

Ce paragraphe ne concerne que les lois DD_CFC et DD_CC : pour la loi DD_FAT, la matrice jacobienne n’a pas été programmée et on utilise la matrice calculée par perturbation. La matrice jacobienne du système s’écrit sous la forme:

\(J=\left[\begin{array}{cc}{(\frac{\partial {R}_{i}^{(1)}}{\partial {\sigma}_{j}})}_{i=1,6;j=1,6}& {(\frac{\partial {R}_{i}^{(1)}}{\partial \Delta {\omega}_{s}})}_{i=1,6;s=1,12}\\ {(\frac{\partial {R}_{s}^{(2)}}{\partial {\sigma}_{j}})}_{s=1,12;j=1,6}& {(\frac{\partial {R}_{s}^{(2)}}{\partial \Delta {\omega}_{r}})}_{s=1,12;r=1,12}\end{array}\right]\) de dimension \(18\times 18\)

• \(\begin{array}{c}\frac{\partial {R}_{i}^{(1)}}{\partial {\sigma}_{j}}={\left({\Lambda}^{-1}\right)}_{ij}+\sum_{s}{({\mu}_{s})}_{i}\frac{\partial \left(\Delta {p}_{s}\frac{{\tau}_{s}}{∣{\tau}_{s}∣}\right)}{\partial {\sigma}_{j}}={\left({\Lambda}^{-1}\right)}_{ij}+\sum_{s}{({\mu}_{s})}_{i}\frac{{\tau}_{s}}{∣{\tau}_{s}∣}\frac{\partial \Delta {p}_{s}}{\partial {\tau}_{s}}\frac{\partial {\tau}_{s}}{\partial {\sigma}_{j}}=\\ {\left({\Lambda}^{-1}\right)}_{ij}+\sum_{s}{({\mu}_{s})}_{i}{({\mu}_{s})}_{j}\frac{{\tau}_{s}}{∣{\tau}_{s}∣}\frac{\partial \Delta {p}_{s}}{\partial {\tau}_{s}}\end{array}\)

avec \(\frac{\partial \Delta {p}_{s}}{\partial {\tau}_{s}}=n\frac{\left(\Delta {p}_{s}+\dot{{\gamma}_{0}}\Delta t\right)}{{\tau}_{s}}\)

en effet : \(\frac{\partial \Delta {p}_{s}}{\partial {\tau}_{s}}=\frac{n\dot{{\gamma}_{0}}\Delta t}{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}(\omega )}{(\frac{∣{\tau}_{s}∣}{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}(\omega )})}^{n-1}\cdot \frac{{\tau}_{s}}{∣{\tau}_{s}∣}=n\frac{(\Delta {p}_{s}+\dot{{\gamma}_{0}}\Delta t)}{∣{\tau}_{s}∣}\cdot \frac{{\tau}_{s}}{∣{\tau}_{s}∣}\)

car \({\tau}_{s}^{\mathrm{forest}}(\omega )\) ne dépend pas de \({\tau}_{s}\) et \(\Delta {p}_{s}+\dot{{\gamma}_{0}}\Delta t=\dot{{\gamma}_{0}}\Delta t{(\frac{\mid {\tau}_{s}\mid }{{\tau}_{f}+{\tau}_{s}^{\mathrm{forest}}(\omega )})}^{n}\)

• \(\frac{\partial {R}_{i}^{(1)}}{\partial \Delta {\omega}_{s}}=\sum_{u=1,12}{({\mu}_{u})}_{i}\frac{{\tau}_{u}}{∣{\tau}_{u}∣}\frac{\partial \Delta {p}_{u}}{\partial \Delta {\omega}_{s}}\)

avec \(\frac{\partial \Delta {p}_{u}}{\partial \Delta {\omega}_{s}}=-n\dot{{\gamma}_{0}}\Delta t∣{\tau}_{u}∣{(\frac{∣{\tau}_{u}∣}{{\tau}_{f}+{\tau}_{u}^{\mathrm{forest}}(\omega )})}^{n-1}\frac{\frac{\partial {\tau}_{u}^{\mathrm{forest}}}{\partial \Delta {\omega}_{s}}}{{({\tau}_{f}+{\tau}_{u}^{\mathrm{forest}}(\omega ))}^{2}}\) ou encore

\(\frac{\partial \Delta {p}_{u}}{\partial \Delta {\omega}_{s}}=-n\dot{{\gamma}_{0}}\Delta t{(\frac{∣{\tau}_{u}∣}{{\tau}_{f}+{\tau}_{u}^{\mathrm{forest}}(\omega )})}^{n}\frac{\frac{\partial {\tau}_{u}^{\mathrm{forest}}}{\partial \Delta {\omega}_{s}}}{({\tau}_{f}+{\tau}_{u}^{\mathrm{forest}}(\omega ))}\)

\(\frac{\partial \Delta {p}_{u}}{\partial \Delta {\omega}_{s}}=-n\frac{(\Delta {p}_{u}+\dot{{\gamma}_{0}}\Delta t)}{({\tau}_{f}+{\tau}_{u}^{\mathrm{forest}}(\omega ))}\frac{\partial {\tau}_{u}^{\mathrm{forest}}}{\partial \Delta {\omega}_{s}}\)

Il reste à exprimer : \(\frac{\partial {\tau}_{u}^{\mathit{forest}}}{\partial \Delta {\omega}_{s}}\) . Sans irradiation, on obtient:

\(\frac{\partial {\tau}_{u}^{\mathrm{forest}}}{\partial \Delta {\omega}_{s}}=\mu \frac{\partial C(\omega )}{\partial \Delta {\omega}_{s}}\cdot \sqrt{\sum_{j=1,12}{a}_{\mathrm{uj}}{\omega}_{j}}+\mu C(\omega )\frac{{a}_{\mathrm{us}}}{2\sqrt{\sum_{j=1,12}{a}_{\mathrm{uj}}{\omega}_{j}}}\frac{\langle {\omega}_{s}\rangle }{{\omega}_{s}}\)

et avec irradiation: \(\frac{\partial ({\tau}_{u}^{\mathit{forest}})}{\partial (\Delta {\omega}_{s})}=\frac{{\mu}^{2}C(\omega )}{2{\tau}_{u}^{\mathit{forest}}}\left(2\frac{\partial C(\omega )}{\Delta {\omega}_{s}}\sum_{j=1,12}{a}_{\mathit{uj}}{\omega}_{j}+C(\omega ){a}_{\mathit{us}}\right)\)

avec \(\frac{\partial C(\omega )}{\partial \Delta {\omega}_{s}}=\frac{0.8}{2\ln(\alpha b\sqrt{{\rho}_{\mathrm{ref}}})}\frac{1}{\sum_{k=1,12}{\omega}_{k}}\frac{\langle {\omega}_{s}\rangle }{{\omega}_{s}}\)

ce qui fournit la dérivée : \(\frac{\partial {R}_{i}^{(1)}}{\partial \Delta {\omega}_{s}}=-\sum_{u=1,12}{({\mu}_{u})}_{i}\frac{{\tau}_{u}}{∣{\tau}_{u}∣}n\frac{(\Delta {p}_{u}+\dot{{\gamma}_{0}}\Delta t)}{({\tau}_{f}+{\tau}_{u}^{\mathrm{forest}}(\omega ))}\frac{\partial {\tau}_{u}^{\mathrm{forest}}}{\partial \Delta {\omega}_{s}}\)

• \(\frac{\partial {R}_{s}^{2}}{\partial {\sigma}_{i}}=-{({\mu}_{s})}_{i}\frac{\partial \Delta {p}_{s}}{\partial {\tau}_{s}}{h}_{s}(\omega )\)

avec \(\frac{\partial \Delta {p}_{s}}{\partial {\tau}_{s}}=n\frac{\left(\Delta {p}_{s}+\dot{{\gamma}_{0}}\Delta t\right)}{{\tau}_{s}}\)

• \(\frac{\partial {R}_{s}^{(2)}}{\partial {\omega}_{r}}={\delta}_{\mathrm{sr}}-\Delta {p}_{s}\frac{\partial {h}_{s}(\omega )}{\partial {\omega}_{r}}-{h}_{s}(\omega )\frac{\partial \Delta {p}_{s}}{\partial {\omega}_{r}}\)

• On a déjà calculé \(\frac{\partial \Delta {p}_{s}}{\partial \Delta {\omega}_{r}}\) . Il reste à dériver : \(\frac{\partial {h}_{s}(\omega )}{\partial \Delta {\omega}_{r}}\) :

  • \(\frac{\partial}{\partial \Delta {\omega}_{r}}(A\frac{\sum_{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}{\omega}_{j}}{\sum_{j=1,12}\sqrt{{a}_{\mathrm{sj}}{\omega}_{j}}}+BC(\omega )\sum_{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}{\omega}_{j}}-\frac{y}{b}{\omega}_{s})=A\frac{\partial {T}_{1}(\omega )}{\partial \Delta {\omega}_{r}}+B\frac{\partial {T}_{2}(\omega )}{\partial \Delta {\omega}_{r}}-\frac{y}{b}{\delta}_{\mathrm{sr}}\frac{\langle {\omega}_{r}\rangle }{{\omega}_{r}}\)

\(\frac{\partial}{\partial \Delta {\omega}_{r}}({T}_{1}(\omega ))=\frac{\sqrt{{a}_{\mathrm{sr}}}}{\sum_{j=1,12}\sqrt{{a}_{\mathrm{sj}}{\omega}_{j}}}\cdot {I}_{s}^{\mathrm{forest}}(r)\frac{\langle {\omega}_{r}\rangle }{{\omega}_{r}}-\frac{{a}_{\mathrm{sr}}\sum_{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}{\omega}_{j}}{2\sqrt{{a}_{\mathrm{sr}}{\omega}_{r}}{(\sum_{j=1,12}\sqrt{{a}_{\mathrm{sj}}{\omega}_{j}})}^{2}}\frac{\langle {\omega}_{r}\rangle }{{\omega}_{r}}\)

\(\frac{\partial}{\partial \Delta {\omega}_{r}}{T}_{2}(\omega )=\frac{\partial C}{\partial \Delta {\omega}_{r}}\sum_{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}{\omega}_{j}}+C(\omega )\frac{{a}_{\mathrm{sr}}}{2\sqrt{{a}_{\mathrm{sr}}{\omega}_{r}}}{I}_{s}^{\mathrm{copla}}(r)\frac{\langle {\omega}_{r}\rangle }{{\omega}_{r}}\)

\({I}_{s}^{\mathrm{copla}}(r)=1\) si \(r\in \mathrm{copla}(s)\) , \(=0\) sinon et \({I}_{s}^{\mathrm{forest}}(r)=1\) si \(r\in \mathrm{forest}(s)\) , \(0\) sinon

Résolution implicite – modèles MONO_DD_CC / MONO_DD_CC_IRRA

On choisit d’adimensionner les inconnues en posant \(\omega =\rho {b}^{2}\) . L’algorithme proposé (cf. 17 ) consiste à éliminer la dépendance de \({R}^{s}\) à \({\tau}_{\mathrm{eff}}^{s}\) en approchant \(\Delta G\) par :

1. \(\Delta {G}_{\mathit{app}}=min\left(\Delta {G}_{0};kT\ln\left[\frac{{\omega}_{m}^{s}H}{\sqrt{{\omega}_{\mathit{tot}}^{s}}{\dot{\epsilon}}_{\mathit{eq}}}\right]\right)\) avec \({\omega}_{\mathit{tot}}^{s}(\omega )={\omega}_{F}^{s}(\omega )=\sum_{j\ne s}{\omega}^{j}+{\omega}_{\mathit{irr}}^{s}\)

dans cette expression, \(\dot{{\varepsilon}_{\mathrm{eq}}}\) est un paramètre fixé.

2. \({R}^{s}(\omega )=\frac{\mu b}{2{\tau}_{0}{\left(1-\frac{\Delta {G}_{\mathit{app}}}{\Delta {G}_{0}}\right)}^{2}}\) en vérifiant que \(\Delta {G}_{\mathit{app}}\le \Delta {G}_{0}\) (sinon, \({R}^{s}\) est choisi grand)

3. \(\frac{1}{{\lambda}^{s}(\omega )+D}=min\left(\frac{\sqrt{{\omega}_{\mathit{tot}}^{s}(\omega )}}{b};(D+2{R}^{s}(\omega ))\frac{{\omega}_{\mathit{tot}}^{s}(\omega )}{{b}^{2}}\right)\)

4. \({l}^{s}(\omega )=max\left({\lambda}^{s}(\omega )-2{\alpha}_{\mathit{AT}}^{s}(\omega ){R}^{s}(\omega );{l}_{c}\right)\)

5. \({\alpha}_{\mathit{AT}}^{s}(\omega )=\sqrt{\sum_{j\ne s}{h}^{\mathit{sj}}\frac{{\omega}^{j}}{{\omega}_{\mathit{tot}}^{s}}+{a}_{\mathit{irr}}\frac{{\omega}_{\mathit{irr}}^{s}}{{\omega}_{\mathit{tot}}^{s}}}\) avec \({h}_{\mathrm{sj}}\) le terme courant de la matrice d’interaction

6. \({\tau}_{\text{LT}}^{s}(\omega )=\max\left(0;{\alpha}_{\mathit{AT}}^{s}(\omega )\mu b\left(\frac{1}{{\lambda}^{s}(\omega )}-\frac{1}{2{\alpha}_{\mathit{AT}}^{s}(\omega ){R}^{s}(\omega )+{l}_{c}(T)}\right)\right)\)

7. \({\tau}_{\mathit{LR}}^{s}(\omega )=\mu \sqrt{{h}^{\mathit{ss}}{\omega}^{s}}\)

8. \({\tau}_{\mathit{eff}}^{s}(\omega ,{\tau}^{s})=∣{\tau}^{s}∣-{\tau}_{c}(\omega )\) et \({\tau}_{c}^{s}={\tau}_{F}+\text{}\sqrt{({\tau}_{\mathit{LR}}^{s}{(\omega )}^{2}+{\tau}_{\text{LT}}^{s}{(\omega )}^{2})}\)

9. \({\Delta \gamma }_{\mathit{nuc}}^{S}(\omega ,{\tau}^{s})=\frac{\Delta t{\omega}_{m}^{s}H}{b}\cdot {l}^{s}(\omega )\exp\left(\frac{-\Delta G(\omega ,{\tau}^{s})}{{k}_{B}T}\right)\mathit{sgn}({\tau}_{s})\) avec \(\Delta G(\omega ,{\tau}^{s})=\Delta {G}_{0}\left(1-\sqrt{\frac{\text{<}{\tau}_{\mathit{eff}}^{s}(\omega ,{\tau}^{s})\text{>}}{{\tau}_{0}}}\right)\) si \({\tau}_{\mathit{eff}}^{s}\le {\tau}_{0}\) , sinon \(\Delta G=\Delta {G}_{0}\)

10. \({\Delta \gamma }_{\mathit{prob}}^{s}(\omega ,{\tau}^{s})=\Delta t\dot{{\gamma}_{0}}{\left(\frac{∣{\tau}^{s}∣}{{\tau}_{c}^{s}(\omega )}\right)}^{n}\mathit{sgn}({\tau}^{s})\)

11. \(\frac{1}{{\Delta \gamma }^{s}(\omega ,{\tau}^{s})}=\frac{1}{{\Delta \gamma }_{\mathit{nuc}}^{s}(\omega ,{\tau}^{s})}+\frac{1}{{\Delta \gamma }_{\mathit{prob}}^{s}(\omega ,{\tau}^{s})}\)

12. \(\Delta {\omega}^{s}(\omega ,{\tau}^{s})=∣\Delta {\gamma}^{s}(\omega ,{\tau}^{s})∣\left(\frac{b}{{d}_{\mathit{lath}}}+\frac{{c}_{\mathit{eff}}^{s}(\omega ,{\tau}^{s})\sqrt{{h}^{\mathit{ss}}{\omega}^{s}}}{{K}_{\mathit{self}}}+\frac{{\alpha}_{\mathit{AT}}^{s}(\omega ){c}_{\mathit{eff}}^{s}(\omega ,{\tau}^{s}){\lambda}_{s}(\omega ){\omega}_{\mathit{tot}}^{s}(\omega )}{b{K}_{f}}-{y}^{s}(\omega ,{\tau}^{s})\frac{{\omega}^{s}}{b}\right)\) avec

\({c}_{\mathit{eff}}^{s}(\omega ,{\tau}^{s})=\left(1-\frac{\text{<}{\tau}_{\mathit{eff}}^{s}(\omega ,{\tau}^{s})\text{>}}{{\tau}_{0}}\right)\) \(\frac{1}{{y}^{s}}(\omega ,{\tau}^{s})=\frac{1}{{y}_{\mathit{AT}}^{s}}+\frac{2\pi {\tau}_{\mathit{eff}}^{s}(\omega ,{\tau}^{s})}{\mu b}\)

Dans le cas DD_CC_IRRA, \(\dot{{\omega}_{\mathit{irr}}^{s}}=-\xi {\omega}_{\mathit{irr}}^{s}∣\dot{{\gamma}^{s}}(\omega )∣\) , soit : \({\omega}_{\mathit{irr}}^{s}(t+\Delta t)={\omega}_{\mathit{irr}}^{s}(t)\exp(-\xi ∣\Delta {\gamma}_{s}∣)\)

Il faut donc résoudre un système de 12 équations en \(\Delta {\omega}^{s}\)

Remarque: Si on veut résoudre complètement les équations ci-dessus sans l’approximation \(\Delta {G}_{\mathit{app}}\) , il faut alors ajouter les 12 inconnues \({\tau}_{\mathit{eff}}\) . En effet la dépendance à \({\tau}_{\mathit{eff}}\) est implicite dans l’expression de \({R}_{s}\) si on remplace \(\Delta {G}_{\mathit{app}}\) par \(\Delta G\) .

La matrice jacobienne du système peut être alors calculée pour l’intégration par la méthode de Newton .

Matrice jacobienne

Avec l’approximation choisie, le système d’équations à résoudre est:

\({R}_{i}^{(1)}(\sigma ,\omega )={({\Lambda}^{-1}\sigma )}_{i}-{({({\Lambda}^{-})}^{\text{-1}}{\sigma}^{-})}_{i}-\Delta {\varepsilon}_{i}+\Delta {\varepsilon}_{i}^{\text{th}}+\sum_{s}({({\mu}_{s})}_{i}\Delta {p}_{s}(\sigma ,\omega ).\frac{{\tau}_{s}}{\mid {\tau}_{s}\mid })=0\) ou encore

\({R}_{i}^{(1)}\left(\Delta {\epsilon}^{\mathrm{el}},\omega \right)=\Delta {\epsilon}_{i}^{\mathit{el}}-\Delta {\epsilon}_{i}+\Delta {\epsilon}_{i}^{\text{th}}+\sum_{s}\left({\mu}_{i}^{s}\Delta {\gamma}^{s}\left({\tau}^{s},\omega \right)\right)=0\) avec \({\tau}^{s}=\sigma :{\mu}^{s}=\Lambda {\epsilon}^{\mathit{el}}:{\mu}^{s}\)

\({R}_{s}^{(2)}\left({\tau}^{s},\omega \right)=\Delta {\omega}_{s}-\Delta {p}_{s}\left(\omega ,{\tau}^{s}\right){H}_{s}\left(\omega ,{\tau}^{s}\right)=0\) avec \(\Delta {p}_{s}(\omega ,{\tau}^{s})=∣(\Delta {\gamma}_{s}(\omega ,{\tau}^{s}))∣\) et

\({H}^{s}(\omega ,{\tau}^{s})=\left(\frac{b}{{d}_{\mathit{lath}}}+\frac{{c}_{\mathit{eff}}^{s}(\omega ,{\tau}^{s})\sqrt{{h}^{\mathit{ss}}{\omega}^{s}}}{{K}_{\mathit{self}}}+\frac{{\alpha}_{\mathit{AT}}^{s}(\omega ){c}_{\mathit{eff}}^{s}(\omega ,{\tau}^{s}){\lambda}_{s}(\omega ){\omega}_{\mathit{tot}}^{s}(\omega )}{b{K}_{f}}-{y}^{s}(\omega ,{\tau}^{s})\frac{{\omega}^{s}}{b}\right)\)

La matrice jacobienne du système s’écrit sous la forme:

\(J=\left[\begin{array}{cc}{\left(\frac{\partial {R}_{i}^{(1)}}{\partial {\epsilon}_{j}^{\mathit{el}}}\right)}_{i=1,6;j=1,6}& {\left(\frac{\partial {R}_{i}^{(1)}}{\partial \Delta {\omega}_{s}}\right)}_{i=1,6;s=1,12}\\ {\left(\frac{\partial {R}_{s}^{(2)}}{\partial {\epsilon}_{j}^{\mathit{el}}}\right)}_{s=1,12;j=1,6}& {\left(\frac{\partial {R}_{r}^{(2)}}{\partial \Delta {\omega}_{s}}\right)}_{r=1,12;s=1,12}\end{array}\right]\) de dimension \(18\times 18\)

\(\frac{\partial {R}_{i}^{(1)}}{\partial {\epsilon}_{j}^{\mathit{el}}}={\delta}_{ij}+\sum_{s}{\mu}_{i}^{s}\frac{\partial {\tau}^{s}}{\partial {\sigma}_{j}}={\delta}_{ij}+\sum_{s}\frac{\partial \Delta {\gamma}^{s}}{\partial {\tau}^{s}}{\mu}_{i}^{s}{\mu}_{k}^{s}{\Lambda}_{\mathit{kj}}\)

avec \(\frac{\partial \Delta {\gamma}^{s}}{\partial {\tau}^{s}}={(\Delta {\gamma}^{s})}^{2}\left[\frac{1}{{(\Delta {\gamma}_{\mathit{nuc}}^{s})}^{2}}\frac{\partial \Delta {\gamma}_{\mathit{nuc}}^{s}}{\partial {\tau}_{s}}+\frac{1}{{(\Delta {\gamma}_{\mathit{prob}}^{s})}^{2}}\frac{\partial \Delta {\gamma}_{\mathit{prob}}^{s}}{\partial {\tau}_{s}}\right]\)

\(\frac{\partial \Delta {\gamma}_{\mathit{nuc}}^{s}}{\partial {\tau}_{s}}=\frac{\Delta {\gamma}_{\mathit{nuc}}^{s}\Delta {G}_{0}}{2{k}_{b}T\sqrt{{\tau}_{0}{\tau}_{\mathit{eff}}}}\mathit{sgn}({\tau}^{s})\) \(\mathit{si}{\tau}_{\mathit{eff}}>0\) , 0 sinon

\(\frac{\partial \Delta {\gamma}_{\mathit{prob}}^{s}}{\partial {\tau}_{s}}=\frac{n\Delta {\gamma}_{\mathit{prob}}^{s}}{∣{\tau}_{s}∣}\)

\(\frac{\partial {{R}^{s}}^{(2)}}{\partial {\epsilon}_{i}^{\mathit{el}}}=-\left[\frac{\partial \Delta {p}^{s}}{\partial {\tau}^{s}}{H}^{s}+\Delta {p}^{s}\frac{\partial {H}^{s}}{\partial {\tau}^{s}}\right]{\mu}_{k}^{s}{\Lambda}_{\mathit{kj}}\)

avec \(\frac{\partial \Delta {p}_{s}}{\partial {\tau}_{s}}=\frac{\partial \Delta {p}^{s}}{\partial \Delta {\gamma}^{s}}\frac{\partial \Delta {\gamma}^{s}}{\partial {\tau}_{s}}=\frac{\Delta {\gamma}^{s}}{∣\Delta {\gamma}^{s}∣}\frac{\partial \Delta {\gamma}^{s}}{\partial {\tau}_{s}}\) \(\frac{\partial \Delta {\gamma}^{s}}{\partial {\tau}_{s}}\) déjà calculé précédemment.

et \(\frac{\partial {H}^{s}}{\partial {\tau}^{s}}=\left(\frac{\partial {c}_{\mathit{eff}}^{s}}{\partial {\tau}^{s}}\frac{\sqrt{{h}^{\mathit{ss}}{\omega}^{s}}}{{K}_{\mathit{self}}}+\frac{\partial {c}_{\mathit{eff}}^{s}}{\partial {\tau}^{s}}\frac{{\alpha}_{\mathit{AT}}^{s}(\omega ){\lambda}_{s}(\omega ){\omega}_{\mathit{tot}}^{s}(\omega )}{b{K}_{f}}-\frac{\partial {y}^{s}}{\partial {\tau}^{s}}\frac{{\omega}^{s}}{b}\right)\)

\(\frac{\partial {c}_{\mathit{eff}}^{s}}{\partial {\tau}^{s}}=\frac{-\partial {\tau}_{\mathit{eff}}^{s}}{\partial {\tau}^{s}}\frac{1}{{\tau}_{0}}=\frac{-{\tau}^{s}}{∣{\tau}^{s}∣}\frac{1}{{\tau}_{0}}=\frac{-\mathit{sgn}({\tau}^{s})}{{\tau}_{0}}\) si \({\tau}_{\mathit{eff}}^{s}>0\) , 0 sinon

\(\frac{\partial {y}^{s}}{\partial {\tau}^{s}}=-{({y}^{s})}^{2}\frac{2\pi }{\mu b}\frac{\partial {\tau}_{\mathit{eff}}^{s}}{\partial {\tau}^{s}}=-{({y}^{s})}^{2}\frac{2\pi }{\mu b}\mathit{sgn}({\tau}^{s})\)

\(\frac{\partial {R}_{i}^{(1)}}{\partial \Delta {\omega}^{s}}=\sum_{r=1,12}{\mu}_{i}^{r}\frac{\partial \Delta {\gamma}^{r}}{\partial \Delta {\omega}^{s}}\)

avec \(\frac{\partial \Delta {\gamma}^{r}}{\partial \Delta {\omega}^{s}}={(\Delta {\gamma}^{r})}^{2}\left[\frac{1}{{(\Delta {\gamma}_{\mathit{nuc}}^{r})}^{2}}\frac{\partial \Delta {\gamma}_{\mathit{nuc}}^{r}}{\partial \Delta {\omega}^{s}}+\frac{1}{{(\Delta {\gamma}_{\mathit{prob}}^{r})}^{2}}\frac{\partial \Delta {\gamma}_{\mathit{prob}}^{r}}{\partial \Delta {\omega}^{s}}\right]\)

Les différents termes intervenant ci-dessus sont:

\(\frac{\partial \Delta {\gamma}_{\mathit{nuc}}^{r}}{\partial \Delta {\omega}^{s}}=\Delta {\gamma}_{\mathit{nuc}}^{r}\left[\frac{1}{{l}^{r}}\frac{\partial {l}^{r}}{\partial \Delta {\omega}^{s}}-\frac{1}{{k}_{b}T}\frac{\partial \Delta G}{\partial \Delta {\omega}^{s}}\right]\)

avec

\(\frac{\partial {l}^{r}}{\partial \Delta {\omega}^{s}}=\frac{\partial {\lambda}^{r}}{\partial \Delta {\omega}^{s}}-2\frac{\partial {\alpha}_{\mathit{AT}}^{r}}{\partial \Delta {\omega}^{s}}{R}^{r}-2{\alpha}_{\mathit{AT}}^{r}\frac{\partial {R}^{r}}{\partial \Delta {\omega}^{s}}\mathit{si}{\lambda}^{s}(\omega )-2{\alpha}_{\mathit{AT}}^{s}(\omega ){R}^{s}(\omega )>{l}_{c};0\mathit{sinon}\)

Le calcul de \(\frac{\partial {\lambda}^{r}}{\partial \Delta {\omega}^{s}}\) se décline en deux possibilités :

si \(min\left(\frac{\sqrt{{\omega}_{\mathit{tot}}^{r}(\omega )}}{b};(D+2{R}^{r}(\omega ))\frac{{\omega}_{\mathit{tot}}^{r}(\omega )}{{b}^{2}}\right)=\frac{\sqrt{{\omega}_{\mathit{tot}}^{r}(\omega )}}{b}\) alors \(\frac{\partial {\lambda}^{r}}{\partial \Delta {\omega}^{s}}=\frac{-b}{2{({\omega}^{r})}^{3/2}}(1-{\delta}_{\mathit{rs}})\)

sinon \(\frac{\partial {\lambda}^{r}}{\partial \Delta {\omega}^{s}}=\frac{-{b}^{2}}{{(D+2{R}^{r})}^{2}{({\omega}_{\mathit{tot}}^{r})}^{2}}\left[2{\omega}_{\mathit{tot}}^{r}\frac{\partial {R}^{r}}{\partial \Delta {\omega}^{s}}+(D+2{R}^{r})(1-{\delta}_{\mathit{rs}})\right]\)

\(\frac{\partial {\alpha}_{\mathit{AT}}^{r}}{\partial \Delta {\omega}_{s}}=\frac{1}{2{\alpha}_{\mathit{AT}}^{r}}\frac{1-{\delta}_{\mathit{rs}}}{{\omega}_{\mathit{tot}}^{r}}\left[{h}^{\mathit{rs}}-{({\alpha}_{\mathit{AT}}^{r})}^{2}\right]\)

\(\frac{\partial {R}^{r}}{\partial \Delta {\omega}_{s}}=\frac{2{R}^{r}}{\Delta {G}_{0}-\Delta {G}_{\mathit{app}}}\frac{\partial \Delta {G}_{\mathit{app}}}{\partial \Delta {\omega}_{s}}\)

\(\frac{\partial \Delta {G}_{\mathit{app}}}{\partial \Delta {\omega}_{s}}=\frac{-kT(1-{\delta}_{\mathit{rs}})}{2{\omega}_{\mathit{tot}}^{r}}\)

\(\frac{\partial \Delta G}{\partial \Delta {\omega}^{s}}=\frac{\Delta {G}_{0}}{2\sqrt{{\tau}_{0}{\tau}_{\mathit{eff}}}}\frac{\partial {\tau}_{c}^{r}}{\partial \Delta {\omega}^{s}}\)

\(\frac{\partial \Delta {\gamma}_{\mathit{prob}}^{r}}{\partial {\tau}^{s}}=-\Delta {\gamma}_{\mathit{prob}}^{r}\frac{n}{∣{\tau}_{c}^{r}∣}\frac{\partial {\tau}_{c}^{r}}{\partial \Delta {\omega}^{s}}\)

avec \(\frac{\partial {\tau}_{c}^{r}}{\partial \Delta {\omega}^{s}}=\frac{{\tau}_{\mathit{LR}}^{r}\frac{\partial {\tau}_{\mathit{LR}}^{r}}{\partial \Delta {\omega}^{s}}+{\tau}_{\text{LT}}^{r}\frac{\partial {\tau}_{\text{LT}}^{r}}{\partial \Delta {\omega}^{s}}}{{\tau}_{c}^{r}-{\tau}_{F}}\)

\(\frac{\partial {\tau}_{\mathit{LR}}^{r}}{\partial \Delta {\omega}^{s}}={\delta}_{\mathit{rs}}\frac{\mu}{2}\sqrt{\frac{{h}^{\mathit{rr}}}{{\omega}^{r}}}=\frac{{\delta}_{\mathit{rs}}}{2}\frac{{\tau}_{\mathit{LR}}^{r}}{{\omega}^{r}}\)

\(\frac{\partial {\tau}_{\text{LT}}^{r}}{\partial \Delta {\omega}^{s}}=\frac{\partial {\alpha}_{\mathit{AT}}^{r}}{\partial \Delta {\omega}_{s}}\frac{{\tau}_{\text{LT}}^{r}}{{\alpha}_{\mathit{AT}}^{r}}-\mu b\frac{{\alpha}_{\mathit{AT}}^{r}}{{({\lambda}^{r})}^{2}}\frac{\partial {\lambda}^{r}}{\partial \Delta {\omega}^{s}}+2\mu b\frac{{\alpha}_{\mathit{AT}}^{r}}{{(2{\alpha}_{\mathit{AT}}^{r}{R}^{r}+{L}_{c})}^{2}}\left(\frac{\partial {\alpha}_{\mathit{AT}}^{r}}{\partial \Delta {\omega}_{s}}{R}^{r}+{\alpha}_{\mathit{AT}}^{r}\frac{\partial {R}^{r}}{\partial \Delta {\omega}_{s}}\right)\) si \({\tau}_{\text{LT}}^{s}>0\) , 0

sinon

\(\frac{\partial {R}_{r}^{(2)}}{\partial \Delta {\omega}^{s}}={\delta}_{\mathit{rs}}-\Delta {p}_{r}\frac{\partial {H}_{r}}{\partial \Delta {\omega}_{s}}-{H}_{r}\frac{\partial \Delta {p}_{r}}{\partial \Delta {\omega}_{s}}\)

\(\frac{\partial \Delta {p}_{r}}{\partial \Delta {\omega}_{s}}=\frac{{\gamma}^{r}}{∣{\gamma}^{r}∣}\frac{\partial \Delta {\gamma}^{r}}{\partial \Delta {\omega}_{s}}\) \(\frac{\partial \Delta {\gamma}^{s}}{\partial \Delta {\omega}_{s}}\) a déjà été calculé précédemment.

Il reste à calculer:

\(\begin{array}{c}\frac{\partial {H}_{r}}{\partial \Delta {\omega}_{s}}=\frac{\partial {c}_{\mathit{eff}}^{r}}{\partial \Delta {\omega}_{s}}\frac{\sqrt{{h}^{\mathit{rr}}{\omega}^{r}}}{{K}_{\mathit{self}}}+\frac{{c}_{\mathit{eff}}^{r}\sqrt{{h}^{\mathit{rr}}}{\delta}_{\mathit{rs}}}{2\sqrt{{\omega}^{r}}{K}_{\mathit{self}}}-\frac{\partial {y}^{r}}{\partial \Delta {\omega}_{s}}\frac{{\omega}^{s}}{b}-\frac{{y}^{r}}{b}{\delta}_{\mathit{rs}}\\ \frac{+1}{b{K}_{f}}\left(\frac{\partial {c}_{\mathit{eff}}^{r}}{\partial \Delta {\omega}_{s}}{\alpha}_{\mathit{AT}}^{r}{\lambda}^{r}{\omega}_{\mathit{tot}}^{r}+{c}_{\mathit{eff}}^{r}\frac{\partial {\alpha}_{\mathit{AT}}^{r}}{\partial \Delta {\omega}_{s}}{\lambda}^{r}{\omega}_{\mathit{tot}}^{r}+{c}_{\mathit{eff}}^{r}{\alpha}_{\mathit{AT}}^{r}\frac{\partial {\lambda}^{r}}{\partial \Delta {\omega}_{s}}{\omega}_{\mathit{tot}}^{r}+{c}_{\mathit{eff}}^{r}{\alpha}_{\mathit{AT}}^{r}{\lambda}^{r}(1-{\delta}_{\mathit{rs}})\right)\end{array}\)

avec:

\(\frac{\partial {c}_{\mathit{eff}}^{s}}{\partial \Delta {\omega}^{s}}=\frac{1}{{\tau}_{0}}\frac{\partial {\tau}_{c}^{r}}{\partial \Delta {\omega}^{s}}\)

\(\frac{\partial {y}^{r}}{\partial \Delta {\omega}^{s}}={({y}^{r})}^{2}\frac{2\pi }{\mu b}\frac{\partial {\tau}_{c}^{r}}{\partial \Delta {\omega}^{s}}\)

Algorithme d’intégration implicite en grandes déformations

A partir de \({U}_{n}\) et \(\Delta U\) , on peut calculer \({F}_{n}\) , \(\Delta F\) et \({F}_{n+1}=\Delta F{F}_{n}\) . de plus, on connait les résultats de l’incrément de temps \(n\) : contraintes de Cauchy \({\sigma}_{n}\) variables internes.

Il s’agit de déterminer \({\sigma}_{n+1}\) et les variables internes \({({\gamma}_{s})}_{n+1}\) (ou \({({\alpha}_{s})}_{n+1}\) ) solutions du système de \(6+{n}_{s}\) équations :

\({S}_{n+1}=\Lambda .\frac{1}{2}({{F}_{n+1}^{e}}^{T}{F}_{n+1}^{e}–{I}_{d})\)

et pour chaque système de glissement :

\({(\Delta {\gamma}_{s})}_{n+1}=\Delta tg({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\) ou \(\Delta {\alpha}_{s}=\Delta th({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\)

avec \({\tau}_{s}(S)={M}_{n+1}:{m}_{s}\otimes {n}_{s}=({{F}^{e}}_{n+1}^{T}{{F}^{e}}_{n+1}{(S)}_{n+1}):{m}_{s}\otimes {n}_{s}=\left[(2{\Lambda}^{\text{-1}}S+{I}_{d})S\right]:{m}_{s}\otimes {n}_{s}\)

en effet suivant 7 et 13 , on considère ici des contraintes de Piola-Kirchhoff mesurées à partir de la configuration intermédiaire, donc seule la partie élastique du gradient de transformation intervient.

\({F}_{n+1}^{\text{e}}={F}_{n+1}{({F}_{n+1}^{p})}^{\text{-1}}=\Delta F{F}_{n}^{\text{e}}{(\Delta {F}^{p})}^{\text{-1}}\)

Plusieurs choix sont possibles pour le calcul de \({(\Delta {F}^{p})}^{\text{-1}}\) :

• \({(\Delta {F}^{p})}^{\text{-1}}=\exp(-\sum\Delta {\gamma}_{s}(S,{\omega}_{s}){m}_{s}\otimes {n}_{s})\) 12 qui a l’avantage de conserver la trace nulle des déformations plastiques, mais qui peut présenter un surcoût ;

• \({(\Delta {F}^{p})}^{\text{-1}}={({I}_{d}+\sum\Delta {\gamma}_{s}(S,{\omega}_{s}){m}_{s}\otimes {n}_{s})}^{\text{-1}}{\left(det({I}_{d}+\sum\Delta {\gamma}_{s}(S,{\omega}_{s}){m}_{s}\otimes {n}_{s})\right)}^{-\frac{1}{3}}\) 13

• \({(\Delta {F}^{p})}^{\text{-1}}=({I}_{d}-\sum\Delta {\gamma}_{s}(S,{\omega}_{s}){m}_{s}\otimes {n}_{s}){\left(det({I}_{d}-\sum\Delta {\gamma}_{s}(S,{\omega}_{s}){m}_{s}\otimes {n}_{s})\right)}^{-\frac{1}{3}}\)

Les gradients \({F}^{\text{e}}\) sont stockés comme variables internes (en ôtant le tenseur identité)

En post-traitement on calcule les contraintes de Cauchy : \({\sigma}_{n+1}=\frac{1}{\det{F}_{n+1}^{e}}{F}_{n+1}^{e}{(S)}_{n+1}{{F}_{n+1}^{e}}^{\text{T}}\)

ou les contraintes de Kirchhoff : \({\tau}_{n+1}=\det{F}_{n+1}{\sigma}_{n+1}={F}_{n+1}^{e}{(S)}_{n+1}{{F}_{n+1}^{e}}^{\text{T}}\)

car \(\det{F}_{n+1}=\det{F}_{n+1}^{e}\det{F}_{n+1}^{p}=\det{F}_{n+1}^{e}\)

Le problème à résoudre en chaque point d’intégration a donc pour but, connaissant les variables internes. \({({\gamma}_{s})}_{n}\) \({({\alpha}_{s})}_{n}\) \({({p}_{s})}_{n}\) \({F}_{n}^{e}\) de déterminer \({\sigma}_{n+1}\) et les variables internes \({({\gamma}_{s})}_{n+1}\) ( \({({\alpha}_{s})}_{n+1}\) , \({F}_{n+1}^{e}\) ) solutions du système à \(6+{n}_{s}\) équations :

\({R}_{1}={(\Lambda )}^{\text{-1}}{S}_{n+1}-\frac{1}{2}({{F}_{n+1}^{e}}^{T}{F}_{n+1}^{e}–{I}_{d})=0\)

et pour chaque système de glissement :

\({R}_{2}={(\Delta {\gamma}_{s})}_{n+1}-\Delta tg({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})=0\) ou \({R}_{2}=\Delta {\alpha}_{s}-\Delta t\cdot h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})=0\)

avec \({\tau}_{s}(S)=\left[(2{\Lambda}^{\text{-1}}S+{I}_{d})S\right]:{m}_{s}\otimes {n}_{s}\)

\({F}_{n+1}^{e}=\Delta F{F}_{n}^{e}{(\Delta {F}^{p})}^{\text{-1}}\)

Le résolution de ce système d’équation non linéaires peut se faire classiquement par l’algorithme de Newton [R5.03.14].

Le vecteur des inconnues est composé de : \({S}_{n+1}\) (6 composantes) et \(\Delta {\alpha}_{s}\) (\({n}_{s}\) composantes).

Le calcul de la matrice jacobienne est décrit en annexe.

Initialisation ,

L’algorithme est initialisé de la façon suivante :

On pose \({F}_{n+1}^{\text{e-trial}}={F}_{n+1}{({F}_{n}^{p})}^{\text{-1}}\) ce qui revient à supposer que \({F}_{n+1}^{p}=\Delta {F}^{p}{F}_{n}^{p}={I}_{d}{F}_{n}^{p}\)

On calcule \({R}_{1}\) avec \({F}_{n+1}^{\text{e}}={F}_{n+1}^{\text{e-trial}}\) puis toutes les équations \({R}_{2}\) Si aucun seuil n’est atteint, alors \({F}_{n+1}^{\text{e}}={F}_{n+1}^{\text{e-trial}}\) . Sinon il faut résoudre le système \({R}_{1}\) \({R}_{2}\) .

On peut choisir de garder comme variable interne le tenseur \({F}_{n+1}^{e}\) . pour calculer plus aisément \({F}_{n+1}^{e}=\Delta F{F}_{n}^{e}{(\Delta {F}^{p})}^{\text{-1}}\)

En pratique on stocke \({F}^{e}-{I}_{d}\) dans les variables internes \(V(6+3\ast {n}_{s}+1)\) à \(V(6+3\ast {n}_{s}+9)\)

On conserve aussi \({F}^{p}-{I}_{d}\) dans les variables internes \(V(6+3\ast {n}_{s}+10)\) à \(V(6+3\ast {n}_{s}+19)\)

En post-traitement, les tenseurs \({m}_{s}\) et \({n}_{s}\) peuvent être mis à jour en utilisant la rotation issue de la décomposition polaire \({F}^{e}={R}^{e}{U}^{e}\) : \({m}_{s}={R}^{e}{{m}_{s}}_{0}\) \({n}_{s}={R}^{e}{{n}_{s}}_{0}\)

Pour le calcul de la décomposition polaire, un algorithme exact est proposé dans 12 .

Suivant 13 , il peuvent aussi être transportés dans la configuration actuelle par :

\({m}_{s}={F}^{e}{{m}_{s}}_{0}\) et \({n}_{s}={F}^{e-T}{{n}_{s}}_{0}\)

Critères de convergence utilisés pour les résolutions implicites

Le critère d’arrêt des itérations porte sur la nullité relative du résidu: \(\frac{{∣∣R(Y)∣∣}_{h}}{{∣∣{R}_{\mathrm{ref}}∣∣}_{h}}<{\varepsilon}_{c}\) . où \({\varepsilon}_{c}\) représente la tolérance donnée par RESI_INTE_RELA. Le problème consiste à bien choisir \({R}_{\mathrm{ref}}\) ainsi que la norme \({∣∣.∣∣}_{h}\)

Par homogénéité avec l’algorithme global de Newton utilisé dans STAT_NON_LINE[R5.03.02], et pour éviter des itérations inutiles quand le résidu initial est très petit, on choisit :

\(\frac{\underset{i=1,6}{\max}∣{R}_{i}(Y)∣}{\underset{i=1,6}{\max}∣{E}_{i}^{\mathrm{tr}}∣}<{\varepsilon}_{c}\) avec \({E}^{\mathrm{tr}}={\Lambda}^{\text{-1}}{\Sigma}^{-}+\Delta E-\Delta {E}^{\mathrm{th}}\)

\(\frac{\underset{i=7,n}{\max}∣{R}_{i}(Y)∣}{\underset{i=7,n}{\max}∣{Y}_{i}^{-}+\Delta {Y}_{i}∣}<{\varepsilon}_{c}\)

Les deux critères précédents doivent être vérifiés pour obtenir la convergence.

Si la convergence n’est pas atteinte après le nombre maximum d’itérations, on teste également la stationnarité de la solution:

\(\mid {Y}_{k+1}-{Y}_{k}\mid <\varepsilon\)

La méthode utilisée permet un re-découpage local du pas de temps, soit systématique, soit en cas de non convergence (mot-clé ITER_INTE_PAS) : ceci permet une intégration plus facile, mais en perdant alors le bénéfice d’un opérateur tangent cohérent à l’issue de l’intégration du comportement.

Il est souvent préférable, en cas de non convergence locale au bout de ITER_INTE_MAXI itérations, de procéder à un redécoupage global du pas de temps (en utilisant DEFI_LIST_INST), afin de conserver une matrice tangente cohérente.

Résolution explicite

Une autre méthode de résolution, très simple à mettre en œuvre pour résoudre les équations différentielles du comportement monocristallin est la résolution explicite. Pour qu’elle soit efficace numériquement, il est indispensable de lui associer un contrôle de pas automatique. Comme dans [R5.03.14], on utilise la méthode de Runge et Kutta. Le calcul des variables internes à l’instant

n’est fonction que des valeurs de leurs dérivés \(\frac{\text{dY}}{\text{dt}}=F(Y,t)\) :

\(\Delta Y=\lbrace \begin{array}{}\Delta {\alpha}_{s}\\ \Delta {\gamma}_{s}\\ \Delta {p}_{s}\\ \Delta {E}^{\mathrm{vp}}\end{array}=\lbrace \begin{array}{}h({\tau}_{s}^{-},{\alpha}_{s}^{-},{\gamma}_{s}^{-},{p}_{s}^{-})\\ g({\tau}_{s}^{-},{\alpha}_{s}^{-},{\gamma}_{s}^{-},{p}_{s}^{-})\\ ∣\Delta {\gamma}_{s}∣=f({\tau}_{s}^{-},{\alpha}_{s}^{-},{\gamma}_{s}^{-},{p}_{s}^{-})\\ \sum_{s}{m}_{s}\Delta {\gamma}_{s}\end{array}\)

avec \({\tau}_{s}=\Sigma :{\mu}_{s}=({\Sigma}^{-}+\Lambda (\Delta E-\Delta {E}^{\mathrm{th}}-\Delta {E}^{\mathrm{vp}})):{\mu}_{s}\)

On intègre selon le schéma suivant :

\({Y}_{t+h}={Y}^{\text{(2)}}\) si le critère de précision est satisfait

\({Y}^{(2)}=Y+\frac{h}{2}\left[F(Y,t)+F({Y}^{(1)},t+h)\right]\) avec \({Y}^{(1)}=Y+hF(Y,t)\)

La différence entre \({Y}^{(2)}\) (schéma d’ordre 2) et \({Y}^{(1)}\) (schéma d’ordre 1, Euler) fournit une estimation de l’erreur d’intégration et permet de contrôler la taille du pas de temps

qui est initialisé à \(\Delta {t}_{i}\) pour la première tentative. La stratégie du contrôle du pas est définie à partir d’une norme de l’écart entre les deux méthodes d’intégration : \(∥{Y}^{(2)}-{Y}^{(1)}∥\) et de la précision requise par l’utilisateur \(\eta\) (mot-clé: RESI_INTE_RELA). Le critère retenu est le suivant, où l’on note

:

\(\delta Y(t)=\underset{j=1,N}{\sup}\left\lbrace \frac{∣{y}_{i}^{(2)}-{y}_{i}^{(1)}∣}{\max\left[\varepsilon ,∣{y}_{j}^{(t)}∣\right]}\right\rbrace <\eta\)

Le paramètre \(\varepsilon\) est fixé à 0,001. La précision d’intégration souhaitée

doit être cohérente avec le niveau de précision requis pour l’étape globale.

Si le critère n’est pas vérifié, le pas de temps est re-découpé selon une méthode heuristique (nombre de sous-pas défini par l’utilisateur via le mot-clé ITER_INTER_PAS). Lorsque le pas de temps devient trop faible (h < 1.10–20), le calcul est arrêté avec un message d’erreur.

Algorithme d’intégration explicite en grandes déformations

A chaque instant et en chaque point d’intégration, on peut calculer \({F}_{n}\) , \(\Delta F\) et \({F}_{n+1}=\Delta F{F}_{n}\) . On connait les contraintes de Cauchy \({\sigma}_{n}\) et les variables internes au temps \(n\)

Il s’agit de déterminer \(S\) et les variables internes \(({\gamma}_{s})\) , \(({\alpha}_{s})\) \({F}^{p}\) au temps \(n+1\) . On peut utiliser pour ce faire une méthode de Runge-Kutta [cf. R5.03.14] sur le système d’équations différentielles :

\({\dot{F}}^{p}=(\sum_{s}\dot{{\gamma}_{s}}{m}_{s}\otimes {n}_{s}){F}^{p}\)

\(\dot{{\gamma}_{s}}=g({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\) pour chaque système de glissement

\(\dot{{\alpha}_{s}}=h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\) pour chaque système de glissement

avec : \({\tau}_{s}=M:{m}_{s}\otimes {n}_{s}=({{F}^{e}}^{T}{F}^{e}S):{m}_{s}\otimes {n}_{s}=\left[(2{\Lambda}^{\text{-1}}S+{I}_{d})S\right]:{m}_{s}\otimes {n}_{s}\)

\((S)=\Lambda .\frac{1}{2}({{F}^{e}}^{T}{F}^{e}–{I}_{d})\) et \({F}^{e}=F{({F}^{p})}^{\text{-1}}\)

Variables internes

Cas du monocristal

Les variables internes dans Code_Aster sont dénommées \({V}_{1}\) , \({V}_{2}\) , … \({V}_{p}\) .

Les six premières sont les 6 composantes de la déformation visco-plastique : \({E}_{ij}^{\mathrm{vp}}\) :

\({E}^{\mathrm{vp}}=\sum_{t}(\Delta {E}^{\mathrm{vp}})\) avec \(\Delta {E}^{\mathrm{vp}}=\sum_{s}{\mu}_{s}{\mathrm{\Delta \gamma }}_{s}\)

\({V}_{1}={E}_{xx}^{\mathrm{vp}}\) , \({V}_{2}={E}_{yy}^{\mathrm{vp}}\) , \({V}_{3}={E}_{zz}^{\mathrm{vp}}\) , \({V}_{4}=\sqrt{(2)}{E}_{xy}^{\mathrm{vp}}\) , \({V}_{5}=\sqrt{(2)}{E}_{xz}^{\mathrm{vp}}\) , \({V}_{6}=\sqrt{(2)}{E}_{yz}^{\mathrm{vp}}\)

\({V}_{7}\) , \({V}_{8}\) , \({V}_{9}\) sont les valeurs de \({\alpha}_{1}\) \({\gamma}_{1}\) \({p}_{1}\) pour le système de glissement \(s=1\)

\({V}_{10}\) , \({V}_{11}\) , \({V}_{12}\) correspondent au système \(s=2\) , et ainsi de suite,où:

• \({\alpha}_{s}\) représente la variable cinématique du système \(s\) dans le cas des modèles phénoménologiques, et la densité de dislocations dans un modèle issu de la DD;

• \({\gamma}_{s}\) représente le glissement plastique du système \(s\)

• \({p}_{1}\) représente le glissement plastique cumulé du système \(s\)

Prise en compte de l’irradiation :

• dans le cas DD_CC_IRRA, il faut ajouter \({n}_{\mathrm{irra}}=12\) variables internes : \({V}_{6+{\mathrm{3n}}_{s}+1}\) à \({V}_{6+{\mathrm{3n}}_{s}+12}\) contiennent pour chaque système de glissement la densité de dislocations liée à l’irradiation \({\rho}_{s}^{\mathrm{irr}}\)

• dans le cas DD_CFC_IRRA, il faut ajouter \({n}_{\mathrm{irra}}=24\) variables internes : \({V}_{6+{\mathrm{3n}}_{s}+1}\) à \({V}_{6+{\mathrm{3n}}_{s}+12}\) contiennent pour chaque système de glissement \({\rho}_{s}^{\mathrm{loops}}\ast {b}^{2}\) \({V}_{6+{\mathrm{3n}}_{s}+13}\) à \({V}_{6+{\mathrm{3n}}_{s}+24}\) contiennent pour chaque système de glissement \({\phi }_{s}^{\mathrm{voids}}\)

On stocke ensuite les cissions pour chaque système de glissement: \({\tau}_{1}\) , … \({\tau}_{{n}_{s}}\)

Dans le cas où on prend en compte la rotation du réseau cristallin, il faut ajouter \({n}_{\mathrm{rota}}=16\) variables internes :

• \({V}_{6+{\mathrm{3n}}_{s}+1}\) à \({V}_{6+{\mathrm{3n}}_{s}+9}\) sont les 9 composantes de la matrice de rotation \(Q\) ,

• \({V}_{6+{\mathrm{3n}}_{s}+10}\) à \({V}_{6+{\mathrm{3n}}_{s}+12}\) sont les 3 composantes de \(\Delta {\omega}^{p}\) ,

• \({V}_{6+{\mathrm{3n}}_{s}+13}\) à \({V}_{6+{\mathrm{3n}}_{s}+15}\) sont les 3 composantes de \(\Delta {\omega}^{e}\) ,

• \({V}_{6+{\mathrm{3n}}_{s}+16}\) représente \(\Theta\)

L’antépénultième variable interne est la contrainte de clivage: \(\underset{s}{\max}(\Sigma .n):n\)

L’avant dernière variable interne contient la déformation plastique cumulée globale, définie par :

\({V}_{p-1}=\sum\Delta {E}_{\mathrm{eq}}^{\mathrm{vp}}\) avec \(\Delta {E}_{\mathrm{eq}}^{\mathrm{vp}}=\sqrt{\frac{2}{3}(\Delta {E}^{\mathrm{vp}}:\Delta {E}^{\mathrm{vp}})}\)

La dernière variable interne, \(\mathit{Vp}\) , (\(p=6+{\mathrm{3n}}_{s}+{n}_{\mathrm{rota}}+3\) , \({n}_{s}\) étant le nombre total de systèmes de glissement) est un indicateur de plasticité (seuil dépassé en au moins un système de glissement au pas de temps courant). S’il est nul, il n’y a pas eu d’accroissement de variables internes à l’instant courant. Sinon, il contient le nombre d’itérations de Newton local (pour une résolution implicite) qui ont été nécessaires pour obtenir la convergence.

Cas du polycristal

Les variables internes dans Code_Aster sont dénommées \({V}_{1}\) , \({V}_{2}\) , … \({V}_{p}\) .

Le nombre de variables internes est \(p=7+\mathrm{6m}+\sum_{g=1,m}({\mathrm{3n}}_{s}(g))+\mathrm{6m}+1\) , \(m\) étant le nombre de phases et \({n}_{s}(g)\) étant le nombre de systèmes de glissement de la phase \(g\) ) .

Dans le cas où on prend en compte l’irradiation, le nombre total de variables internes est:

\(p=7+\mathrm{6m}+\sum_{g=1,m}({\mathrm{3n}}_{s}(g))+\mathrm{12m}+\mathrm{6m}+1\)

• Les six premières variables internes sont les composantes de la déformation viscoplastique macroscopique \({E}^{\mathrm{vp}}\) :

\({V}_{1}={E}_{xx}^{\mathrm{vp}}\) , \({V}_{2}={E}_{yy}^{\mathrm{vp}}\) , \({V}_{3}={E}_{zz}^{\mathrm{vp}}\) , \({V}_{4}=\sqrt{(2)}{E}_{xy}^{\mathrm{vp}}\) , \({V}_{5}=\sqrt{(2)}{E}_{xz}^{\mathrm{vp}}\) , \({V}_{6}=\sqrt{(2)}{E}_{yz}^{\mathrm{vp}}\) ;

• la septième est la déformation viscoplastique équivalente cumulée macroscopique \(P\) :

\({V}_{7}=\sum\Delta {E}_{\mathrm{eq}}^{\mathrm{vp}}\) avec \(\Delta {E}_{\mathrm{eq}}^{\mathrm{vp}}=\sqrt{\frac{2}{3}(\Delta {E}^{\mathrm{vp}}:\Delta {E}^{\mathrm{vp}})}\) ;

• puis, pour chaque phase, on trouve les 6 composantes des déformations viscoplastiques ou du tenseur \(\beta\) de la phase : \({\left\lbrace {\varepsilon}_{xx}^{\mathrm{vp}}(g),{\varepsilon}_{yy}^{\mathrm{vp}}(g),{\varepsilon}_{zz}^{\mathrm{vp}}(g),\sqrt{(2)}{\varepsilon}_{xy}^{\mathrm{vp}}(g),\sqrt{(2)}{\varepsilon}_{xz}^{\mathrm{vp}}(g),\sqrt{(2)}{\varepsilon}_{yz}^{\mathrm{vp}}(g)\right\rbrace }_{g=1,m}\) ;

• ensuite, pour chaque phase, et pour chaque système de glissement de la phase, on trouve les valeurs de \({\alpha}_{s}\) \({\gamma}_{s}\) \({p}_{s}\) ;

• si l’irradiation est prise en compte, pour chaque phase, on trouve les 12 densités de dislocation liées à l’irradiation \({\rho}_{\mathit{irr}}^{s}\) ;

• puis, pour chaque phase, on trouve les 6 composantes des contraintes de la phase : \({\left\lbrace {\sigma}_{xx}(g),{\sigma}_{yy}(g),{\sigma}_{zz}(g),\sqrt{(2)}{\sigma}_{xy}(g),\sqrt{(2)}{\sigma}_{xz}(g),\sqrt{(2)}{\sigma}_{yz}(g)\right\rbrace }_{g=1,m}\) ;

• la dernière variable interne est un indicateur de plasticité (seuil dépassé en au moins un système de glissement au pas de temps courant).

Implantation numérique dans Code_Aster

D’une façon générale, les comportements monocristallins sont intégrés aux méthodes de Runge-Kutta pour l’intégration explicite, et à l’environnement «plasti» pour l’intégration implicite [R5.03.14]. Les tenseurs d’orientation des systèmes de glissement sont quant à eux tous définis dans la routine LCMMSG fournissant le tenseur en repère global pour le nième système de la famille fournie dont le nom est fourni par la routine appelante.

Pour ajouter un nouveau comportement de monocristal, ou simplement une nouvelle loi d’écoulement ou d’écrouissage, il convient de définir ses paramètres dans DEFI_MATERIAU. Suivant le cas (écoulement, écrouissage isotrope ou cinématique), il faut ajouter la lecture de ces paramètres dans les routines LCMAFL, LCMAEI, LCMAEC. Pour l’intégration, il suffit d’écrire la définition des accroissements de variables internes dans les routines LCMMFE (écoulement), LCMMFC(écrouissage cinématique) et LCMMFI(écrouissage isotrope), pour que l’intégration explicite fonctionne.

L’intégration implicite utilise également les routines LCMMFE, LCMMFC et LCMMFI. Elle demande en plus de définir les dérivées des équations par rapport aux différentes variables. Les dérivées sont à écrire dans les routines LCMMJF (dérivées de l ‘équation d’écoulement), LCMMJI (dérivées de la relation d’écrouissage isotrope) et LCMMJC (dérivées de la relation d’écrouissage cinématique).

Pour plus de détails sur l’architecture de la résolution des comportements cristallins, voir [D5.04.02].

Utilisation

Ces modèles sont accessibles dans Code_Aster en 3D , déformations planes ( D_PLAN ), contraintes planes ( C_PLAN ,par la méthode de De Borst [R5.03.03]) et axisymétrie ( AXIS ).

Cas du monocristal

Dans le cas de microstructures maillées, les différents grains d’un monocristal étant représentés par des groupes de mailles, il faut affecter les paramètres des matériaux et les comportements des monocristaux ainsi que leurs orientations aux différents grains.

Les valeurs des paramètres des relations de comportement sont fournies à l’aide de la commande DEFI_MATERIAU. Ceci se définit à partir des mots clés MONO_VISC1, MONO_VISC2, MONO_DD_KR, MONO_DD_CFC pour l’écoulement, MONO_ISOT1, MONO_ISOT2, MONO_DD_CFC pour l’écrouissage isotrope et MONO_CINE1, MONO_CINE2 pour l’écrouissage cinématique [U4.43.01]. Par exemple [V6.04.172]:

MATER1=DEFI_MATERIAU(

/ ELAS =_F( E=..., NU=...,

/ ELAS_ORTH=_F( E_L=192500,

E_T=128900,

NU_LT=0.23,

G_LT=74520,),


# RELATIONS D’ECOULEMENT

/ MONO_VISC1=_F(N=10, K=40,C=6333),

/ MONO_VISC2=_F(N=10, K=40,C=6333,D=37,A=121),

/ MONO_DD_KR=_F(...),


# ECROUISSAGE ISOTROPE

/ MONO_ISOT1=_F(R_0=75.5,Q=9.77,B=19.34,H=2.54),

/ MONO_ISOT2=_F(R_0=75.5,Q1=9.77,B1=19.34,H=2.54,Q2=-33.27, B2=5.345,),


# ECROUISSAGE CINEMATIQUE


/ MONO_CINE1=_F(D=36.68),

/ MONO_CINE2=_F(D=36.68, GM=, PM=,),

);

On peut ainsi dissocier, au niveau des données, l’écoulement de l’écrouissage isotrope et de l’écrouissage cinématique.

Il faut maintenant définir le (ou les) type de monocristal étudié. Pour cela, on définit le comportement de façon externe à STAT_NON_LINE, par l’intermédiaire de l’opérateur DEFI_COMPOR, par exemple:

MONO1=DEFI_COMPOR(MONOCRISTAL =(_F( MATER=MATER1,

ECOULEMENT=MONO_VISC1,

MONO_ISOT=MONO_ISOT1,

MONO_CINE=MONO_CINE1,

FAMI_SYST_GLIS=(“CUBIQUE1”,),

_F( MATER=MATER1,

ECOULEMENT=MONO_VISC1,

MONO_ISOT=MONO_ISOT2,

MONO_CINE=MONO_CINE2,

FAMI_SYST_GLIS=”CUBIQUE2”,),

),

_F( MATER=MATER2,

ECOULEMENT=MONO_VISC1,

MONO_ISOT=MONO_ISOT2,

MONO_CINE=MONO_CINE2,

FAMI_SYST_GLIS=”PRISMATIQUE”,),

_F( MATER=MATER2,

ECOULEMENT=MONO_VISC1,

MONO_ISOT=MONO_ISOT2,

MONO_CINE=MONO_CINE2,

FAMI_SYST_GLIS=”BCC24”,), ),

)

L’opérateur DEFI_COMPOR calcule le nombre total de variables internes associé au monocristal.

Enfin, pour réaliser un calcul de micro-structure, il faut donner, grain par grain, ou par groupe de mailles (représentant des ensembles de grains) une orientation, à l’aide du mot-clé MASSIF de AFFE_CARA_ELEM. Par exemple:

ORIELEM = AFFE_CARA_ELEM ( MODELE = MO_MECA,MASSIF =(

_F ( GROUP_MA='GRAIN1', ANGL_REP=(348.0,24.0,172.0),),

_F ( GROUP_MA='GRAIN2', ANGL_REP=( 327.0, 126.0, 335.0),),

_F ( GROUP_MA='GRAIN3', ANGL_REP=( 235.0, 7.0, 184.0),),

_F ( GROUP_MA='GRAIN4', ANGL_REP=( 72.0, 338.0, 73.0), ,

…)

Remarque 1 :

Les orientations des systèmes de glissement peuvent être renseignées en angles nautiques sous le mot clé de ANGL_REPou en angles d’Euler sous le mot clé ANGL_EULER.

Remarque 2:

Pour un même monocristal, les valeurs des paramètres peuvent être différentes d’une famille de systèmes de glissement à l’autre. C’est pourquoi on peut définir un matériau différent par occurrence du mot-clé facteur MONOCRISTAL . Mais dans ce cas, comment fournir transmettre à STAT_NON_LINE l’information stipulant qu’en un point de gauss (tous ceux du groupe de mailles concerné), on a plusieurs matériaux présents? Ceci est possible grâce à une évolution de AFFE_MATERIAU [U4.43.03] et de la structure de données matériau [D4.06.18]): MAT=AFFE_MATERIAU(MAILLAGE=MAIL, AFFE =_F(GROUP_MA=’GRAIN1’, MATER=(MATER1,MATER2),), );

Les autres données du calcul sont identiques à un calcul de structure habituel.

Enfin, dans STAT_NON_LINE, le comportement issu de DEFI_COMPOR est fourni, sous le mot clé COMPORTEMENT via le mot clé COMPOR, obligatoire avec le mot-clé RELATION=’MONOCRISTAL’.

COMPORTEMENT = _F(     RELATION =’MONOCRISTAL’, COMPOR = COMP1

Précisions que pour l’intégration explicite, (ALGO_INTE=’RUNGE_KUTTA’), il est inutile de demander la réactualisation de la matrice tangente puisque celle-ci n’est pas calculée. Pour débuter des itérations de Newton de l’algorithme global, il peut être utile de spécifier PREDICTION=’EXTRAPOLE’ [U4.51.03].

On pourra trouver des exemples d’utilisation dans les tests: SSNV171, SSNV172, SSNV194.

Cas du polycristal

Dans le cas de polycristaux multiphasés, chaque phase correspond à un monocristal. On utilisera donc les paramètres définis précédemment dans DEFI_MATERIAU pour le monocristal. Ici, il s’agit de définir, pour chaque phase, l’orientation, la fraction volumique, le monocristal utilisé, et le type de loi de localisation. Ceci est effectué sous le mot-clé facteur POLYCRISTAL de DEFI_COMPOR.

MONO1=DEFI_COMPOR(MONOCRISTAL=_F(MATER=MATPOLY,

ECOULEMENT=”MONO_VISC2”,

MONO_ISOT=”MONO_ISOT2”,

MONO_CINE=”MONO_CINE1”,

ELAS=”ELAS”,

FAMI_SYST_GLIS=”OCTAEDRIQUE”,),);

POLY1=DEFI_COMPOR(POLYCRISTAL=(_F(MONOCRISTAL=MONO1,

FRAC_VOL=0.025,

ANGL_REP=(-149.676,15.61819,154.676,),),

_F(MONOCRISTAL=MONO1,

FRAC_VOL=0.025,

ANGL_REP=(-150.646,33.864,55.646,),),

_F(MONOCRISTAL=MONO1,

FRAC_VOL=0.025,

ANGL_REP=(-137.138,41.5917,142.138,),),

……

_F(MONOCRISTAL=MONO1,

FRAC_VOL=0.025,

ANGL_REP=(-481.729,35.46958,188.729,),),),

LOCALISATION=”BETA”,

MU_LOCA=80000,,

DL=321.5,

DA=0.216,);

Le mot-clé POLYCRISTAL permet de définir chaque phase par la donnée d’une orientation, (fournie par ANGL_REP ou ANGL_EULER), d’une fraction volumique, d’un mono-cristal (c’est-à-dire un modèle de comportement et des systèmes de glissement).

Le mot-clé LOCALISATION permet de choisir la méthode de localisation pour l’ensemble des phases du polycristal.

Le mot clé MASSIF permet de choisir une orientation globale du polycristal sous le mot clé ANGL_REP si c’est en angles nautiques ou bien sous le mot clé ANGL_EULER si c’est en angles d’Euler.

Enfin, dans STAT_NON_LINE, le comportement issu de DEFI_COMPOR est fourni, sous le mot clé COMPORTEMENT via le mot clé COMPOR, obligatoire avec le mot-clé RELATION=’POLYCRISTAL’.

COMPORTEMENT = _F(     RELATION =’POLYCRISTAL’,

COMPOR = COMP1 )

Exemple

A titre d’exemple de mis en œuvre, on présente ici succinctement un calcul d’agrégat, de forme cubique (volume élémentaire) comprenant 100 grains monocristallins, définis chacun par un groupe de mailles. Le nombre total d’éléments est 86751. Avec des mailles d’ordre 1 (TETRA4) il comporte 15940 nœuds. Avec des mailles d’ordre 2 (TETRA10), il en comporte 121534.

Le chargement consiste en une déformation homogène, appliquée par l’intermédiaire d’un déplacement normal imposé sur une face du cube (direction \(z\) ). On atteint une déformation de 4% en 1s et 50 incréments.

ACIER=DEFI_MATERIAU(ELAS=_F(E =145200.0,NU=0.3,),

MONO_VISC1=_F( N=10., K=40.,C=10.,),

MONO_ISOT2=_F(R_0=75.5,

B1 =19.34,

B2 =5.345,

Q1 =9.77,

Q2 =33.27,

H=0.5),

MONO_CINE1=_F(D=36.68,),);

COEF=DEFI_FONCTION(NOM_PARA ='INST', VALE =(0.0,0.0,1.0,1.0,),);

MAT=AFFE_MATERIAU(MAILLAGE=MAIL, AFFE=_F(TOUT ="OUI", MATER=(ACIER),),);


COMPORT=DEFI_COMPOR(MONOCRISTAL=( _F(MATER =ACIER,

ECOULEMENT="MONO_VISC1",

MONO_ISOT ="MONO_ISOT2",

MONO_CINE ="MONO_CINE1",

ELAS="ELAS",

FAMI_SYST_GLIS='OCTAEDRIQUE',),

),);


ORIELEM = AFFE_CARA_ELEM ( MODELE = MO_MECA, MASSIF =(

_F ( GROUP_MA='GRAIN1', ANGL_REP=(348.0,24.0,172.0), ),

_F ( GROUP_MA='GRAIN2', ANGL_REP=( 327.0, 126.0, 335.0), ),

_F ( GROUP_MA='GRAIN3', ANGL_REP=( 235.0, 7.0, 184.0), ),

……………..

_F ( GROUP_MA='GRAIN99', ANGL_REP=( 201.0, 198.0, 247.0), ),

_F ( GROUP_MA='GRAIN100', ANGL_REP=( 84.0, 349.0, 233.0), ),

))


FO_UZ = DEFI_FONCTION ( NOM_PARA = 'INST',VALE = ( 0.0, 0.0, 1.0, 0.04, ),)


CHME4=AFFE_CHAR_MECA_F(MODELE=MO_MECA, DDL_IMPO=_F(GROUP_NO='HAUT', DZ=FO_UZ,),)


LINST = DEFI_LIST_REEL (DEBUT= 0.,INTERVALLE =(_F ( JUSQU_A = 1.,NOMBRE= 50 ),) )


SIG=STAT_NON_LINE(MODELE =MO_MECA,CARA_ELEM=ORIELEM,CHAM_MATER =MAT,

EXCIT=(_F(CHARGE=CHME1),

_F(CHARGE=CHME2),

_F(CHARGE=CHME3),

_F(CHARGE=CHME4),),

COMPORTEMENT=_F(    RELATION ='MONOCRISTAL',COMPOR =COMPORT),

INCREMENT=_F(LIST_INST=LINST),))

Les figures suivantes représentent des isovaleurs des déformations les contraintes suivant \(z\) . On note la non homogénéité des valeurs, et on peut même discerner le contour des grains.

Pour pouvoir exploiter ce type de résultats, on peut par exemple calculer des champs moyens par grains. Sur la figure suivante, on a représenté les contraintes équivalentes en fonction des déformations plastiques équivalentes pour l’ensemble des grains.

Bibliographie

1. MERIC L., CAILLETAUD G. : «single crystal modeling fort structural calculations» in Journal of Engineering Material and Technology, janvier 1991, vol 113, pp171-182.

2. LECLERCQ S., DIARD, O., PROIX J.M. : “Etude d’impact de l’implantation d’une bibliothèque de lois de comportement et de règles de transition d’échelles» Note EDF R&D HT-26/03/053/A.

3. CAILLETAUD G.: “A micromechanical approach to inelastic behaviour of metals”, Int. J. of Plasticity, 8, pp. 55-73, 1992.

4. PILVIN P.: “The contribution of micromechanical approaches to the modelling of inelastic behaviour of polycrystals”, Int. Conf. on Biaxial / Multiaxial fatigue, France, ESIS/SF2M, pp.31‑46, 1994.

5. BERVEILLER M., ZAOUI A.: “An extension of the self-consistant scheme to plasticity flowing polycrystal” J. Mech. Phys. Solids, 6, pp. 325-344, 1979.

6. G.MONNET “A crystalline plasticity law for austenitic stainless steels”, Note EDF R&D H-B60-2008-04690-EN

7. N.Rupin «implementation of a new constitutive law based on dislocation dynamics for FCC materials» Note EDF-R&D : Note EDF R&D HT24-2010-01128-EN.

8. C.Petry “Procédure de génération de maillages d’agrégats polycristallins” Note EDF R&D H-T24-2008-03481

9. A.ZEGHADI: «Effet de la morphologie 3D et de la taille de grain sur le comportement mécanique d’agrégats polycristallins». Thèse de l’Ecole des Mines de Paris, 2005.

10. N.Rupin:«Déformation à chaud de métaux biphasés. Modélisations théoriques et confrontations expérimentales». Thèse de l’Ecole Polytechnique, 2007.

11. M.Fivel, S.Forest “Plasticité cristalline et transition d’échelle : cas du monocristal”. Techniques de l’ingénieur M4 016

12. «Computational methods for plasticity» Wiley 2008, EA de Souza Neto, D. Peric, DRJ. Owen

13. «A computational procedure for rate-independent crystal plasticity» L.ANAND & M.KOTHARI, J.Mech.Phys.Solids, vol 44, n°4, pp.525-558, 1996

14. 3. Stolz. Milieux continus en transformations finies : hyperélasticité, rupture, élasto-plasticité. Les éditions de l’Ecole Polytechnique, 2002.

15. Homogénéisation en mécanique des matériaux”, tome 2. M.Bornert, P.Gilormini, T.Bertheau

16. 10. SCHWARTZ : “Approche non locale en plasticité cristalline : application à l’étude du comportement mécanique de l’acier AISI 316LN en fatigue oligocyclique”. Thèse de l’Ecole Centrale de Paris, Juin 2011.

17. G.Monnet : “DELIVERABLE D1-2.9. Crystal plasticity constitutive law for irradiated RPV steel” Note EDF R&D H-T27-2011-02738-EN, Décembre 2011.

Historique des versions du document

Version Aster	Auteur(s) ou contributeur(s), organisme	Description des modifications
\(7.4\)	J.M.PROIX, O.DIARD, T.KANIT EDF/R&D	Version initiale
\(8.4\)	J.EL-GHARIB,J.M.PROIX, EDF/R&D	Compléments
\(9.2\)	J.EL-GHARIB,J.M.PROIX, EDF/R&D	Ajout du comportement MONO_DD_KR.
\(9.4\)	J.EL-GHARIB,J.M.PROIX, EDF/R&D	Optimisation du comportement MONO_DD_KR.
\(10.2\)	J.M.PROIX	remarque sur la matrice d’interaction
\(10.3\)	N.RUPIN,J.M.PROIX, EDF/R&D	Ajout du comportement DD_CFC
\(10.5\)	N.RUPIN,J.M.PROIX, EDF/R&D	Ajout de la rotation du réseau cristallin
\(11.1\)	N.RUPIN, F.LATOURTE, J.M.PROIX, EDF/R&D	Ajout des grandes déformations
\(11.1\)	N.RUPIN,J.M.PROIX, EDF/R&D	Fiches 16373 (critère de convergence) et 17422 : correction variables internes du POLYCRISTAL.
\(11.1\)	DE BONNIERES, J.M. PROIX	Fiche 14586 : ajout du comportement DD_FAT
\(11.2\)	J.M. PROIX	Fiche 18398: résorption de ALGO_C_PLAN
\(11.2\)	J.M. PROIX	Fiche 18692 ajout de DD_CC
\(11.3\)	J.M. PROIX	Fiche 19021 ajout des variables internes pour DD_CC_IRRA
\(12.2\)	J.M. PROIX	Fiches 22563: simplification de la loi DD_CC et 22491 coefficient \({\mu}^{\mathit{loca}}\) pour l’homogénéisation,

Expression du Jacobien des équations élasto-visco-plastiques intégrées

Le système à résoudre est de la forme:

\(R(Y)=R(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})=\lbrace \begin{array}{}s(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ e(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ {n}_{s}\left\lbrace \begin{array}{}a(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ g(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ p(\Sigma ,\Delta {E}^{\text{vp}},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\end{array}\right\rbrace \end{array}\)

\(=\lbrace \begin{array}{}{\Lambda}^{\text{-1}}\Sigma -({\Lambda}_{-}^{\text{- 1}}){\Sigma}^{-}-(\Delta E-\Delta {E}^{\text{th}}-\Delta {E}^{\text{vp}})\\ \Delta {E}^{\text{vp}}-{\sum}_{s}{\mu}_{s}\Delta {\gamma}_{s}\\ {n}_{s}\left\lbrace \begin{array}{}\Delta {\alpha}_{s}-h({\tau}_{s}^{+},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ \Delta {\gamma}_{s}-g({\tau}_{s}^{+},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ {\mathrm{\Delta p}}_{s}-f({\tau}_{s}^{+},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\end{array}\right\rbrace \end{array}\)

\(=0\text{avec}{\tau}_{s}^{+}={\Sigma}^{+}:{\mu}_{s}\)

Soit donc à évaluer les termes de l’hypermatrice jacobienne \(J\) à l’instant \(t+\Delta t\)

\(J=\left[\begin{array}{ccccc}\frac{\partial s}{\partial \text{ΔΣ}}& \frac{\partial s}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial s}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial s}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial s}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial e}{\partial \text{ΔΣ}}& \frac{\partial e}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial e}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial e}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial e}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial a}{\partial \text{ΔΣ}}& \frac{\partial a}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial a}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial a}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial a}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial g}{\partial \text{ΔΣ}}& \frac{\partial g}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial g}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial g}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial g}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial p}{\partial \text{ΔΣ}}& \frac{\partial p}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial p}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial p}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial p}{\partial {\mathrm{\Delta p}}_{s}}\end{array}\right]\)

En ce qui concerne la première ligne de la matrice, indépendamment des équations d’écrouissage et d’écoulement, on a:

\(\begin{array}{cccc}\frac{\partial s}{\partial \text{ΔΣ}}={D}^{-1}& \frac{\partial s}{\partial {\mathrm{\Delta E}}^{\text{vp}}}=\text{Id}& \frac{\partial s}{\partial {\mathrm{\Delta \alpha }}_{s}}=& \frac{\partial s}{\partial {\mathrm{\Delta \gamma }}_{s}}=\frac{\partial s}{\partial {\mathrm{\Delta p}}_{s}}\end{array}=0\)

La deuxième ligne peut s’écrire également indépendamment de l’écoulement et des écrouissages:

\(\begin{array}{ccc}\begin{array}{ccc}\frac{\partial e}{\partial \text{ΔΣ}}=0& \frac{\partial e}{\partial {\mathrm{\Delta E}}^{\text{vp}}}=\text{Id}& \frac{\partial e}{\partial {\mathrm{\Delta \alpha }}_{s}}=0\end{array}& \frac{\partial e}{\partial {\mathrm{\Delta \gamma }}_{s}}=-{\mu}_{s}& \frac{\partial e}{\partial {\mathrm{\Delta p}}_{s}}=0\end{array}\)

La première colonne des lignes correspondant aux équations (a), (g) et (p) s’écrit:

\(\begin{array}{c}\frac{\partial a}{\partial \text{ΔΣ}}=\frac{\partial a}{\partial {\mathrm{\Delta \tau }}_{s}}\frac{{\mathrm{\Delta \tau }}_{s}}{\text{ΔΣ}}\\ \frac{\partial g}{\partial \text{ΔΣ}}=\frac{\partial g}{\partial {\mathrm{\Delta \tau }}_{s}}\frac{{\mathrm{\Delta \tau }}_{s}}{\text{ΔΣ}}\\ \frac{\partial p}{\partial \text{ΔΣ}}=\frac{\partial p}{\partial {\mathrm{\Delta \tau }}_{s}}\frac{{\mathrm{\Delta \tau }}_{s}}{\text{ΔΣ}}\end{array}\)

avec

\(\frac{{\mathrm{\Delta \tau }}_{s}}{\text{ΔΣ}}={({\mu}_{s})}^{T}\)

La deuxième colonne est identiquement nulle (du fait de l’équation (e): les relations d’écoulement et d’écrouissage ne peuvent s’exprimer qu’en fonction de \(\Delta {\gamma}_{s}\) et non pas de \({\mathrm{\Delta E}}^{\text{vp}}\) .

Le dernier bloc d’équations, dépend quant à lui des comportements choisis:

\(\begin{array}{ccc}\frac{\partial a}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial a}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial a}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial g}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial g}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial g}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial p}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial p}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial p}{\partial {\mathrm{\Delta p}}_{s}}\end{array}\)

Exemple

Choisissons la relation d’écoulement viscoplastique MONO_VISC1

\(\begin{array}{}(g){\mathrm{\Delta \gamma }}_{s}-{\mathrm{\Delta p}}_{s}\frac{{\tau}_{s}-{\mathrm{c\alpha }}_{s}}{\mid {\tau}_{s}-{\mathrm{c\alpha }}_{s}\mid }=0\\ (p){\mathrm{\Delta p}}_{s}-\mathrm{\Delta t}.{\langle \frac{\mid {\tau}_{s}-{\mathrm{c\alpha }}_{s}\mid -{R}_{s}({p}_{s})}{k}\rangle }^{n}=0\end{array}\)

avec l’écrouissage isotrope MONO_ISOT1: \({R}_{s}({p}_{s})={R}_{0}+Q(\sum_{r=1}^{N}{h}_{\text{sr}}(1-{e}^{-{\text{bp}}_{r}}))\) ,

et un écrouissage cinématique défini par MONO_CINE1

\((a){\mathrm{\Delta \alpha }}_{s}-{\mathrm{\Delta \gamma }}_{s}+{\mathrm{d\alpha }}_{s}{\mathrm{\Delta p}}_{s}=0\)

alors:

\(\begin{array}{c}\frac{\partial a}{\partial {\mathrm{\Delta \tau }}_{s}}=0\\ \frac{\partial g}{\partial {\mathrm{\Delta \tau }}_{s}}=0\\ \frac{\partial p}{\partial {\mathrm{\Delta \tau }}_{s}}=\frac{-\mathrm{n\Delta t}}{{K}^{n}}{\langle \mid {\tau}_{s}-{\mathrm{c\alpha }}_{s}\mid -{R}_{s}({p}_{s})\rangle }^{n-1}\frac{{\tau}_{s}-\mathrm{c\alpha }}{\mid {\tau}_{s}-{\mathrm{c\alpha }}_{s}\mid }\end{array}\)

\(\begin{array}{ccc}\frac{\partial a}{\partial {\mathrm{\Delta \alpha }}_{s}}& =& 1+{\mathrm{d\Delta p}}_{s}\\ \frac{\partial g}{\partial {\mathrm{\Delta \alpha }}_{s}}& =& 0\\ \frac{\partial p}{\partial {\mathrm{\Delta \alpha }}_{s}}& =& \frac{\text{nc}\mathrm{\Delta t}}{{K}^{n}}{\langle \mid {\tau}_{s}-{\mathrm{c\alpha }}_{s}\mid -{R}_{s}({p}_{s})\rangle }^{n-1}\frac{{\tau}_{s}-{\mathrm{c\alpha }}_{s}}{\mid {\tau}_{s}-{\mathrm{c\alpha }}_{s}\mid }\end{array}\)

\(\begin{array}{ccc}\frac{\partial a}{\partial {\mathrm{\Delta \gamma }}_{s}}& =& -1\\ \frac{\partial g}{\partial {\mathrm{\Delta \gamma }}_{s}}& =& 1\\ \frac{\partial p}{\partial {\mathrm{\Delta \gamma }}_{s}}& =& 0\end{array}\)

\(\begin{array}{}\begin{array}{ccc}\frac{\partial a}{\partial {\mathrm{\Delta p}}_{s}}& =& {\mathrm{d\alpha }}_{s}\\ \frac{\partial g}{\partial {\mathrm{\Delta p}}_{s}}& =& \frac{{\tau}_{s}-{\mathrm{c\alpha }}_{s}}{\mid {\tau}_{s}-{\mathrm{c\alpha }}_{s}\mid }\\ \frac{\partial p}{\partial {\mathrm{\Delta p}}_{s}}& =& 1+\frac{\mathrm{n\Delta t}}{{K}^{n}}{\langle \mid {\tau}_{s}-{\mathrm{c\alpha }}_{s}\mid -{R}_{s}({p}_{s})\rangle }^{n-1}\frac{{\text{dR}}_{s}({p}_{s})}{{\mathrm{d\Delta p}}_{s}}\end{array}\\ \frac{{\text{dR}}_{s}({p}_{s})}{{\mathrm{d\Delta p}}_{s}}={\text{Qbh}}_{\text{ss}}{e}^{-{\text{bp}}_{s}}\end{array}\)

et, concernant l’interaction entre systèmes de glissement, il y un seul terme non nul:

\(\begin{array}{}\frac{\partial p}{\partial {\mathrm{\Delta p}}_{r}}=1+\frac{\mathrm{n\Delta t}}{{K}^{n}}{\langle \mid {\tau}_{s}-{\mathrm{c\alpha }}_{s}\mid -{R}_{s}({p}_{s})\rangle }^{n-1}\frac{{\text{dR}}_{s}({p}_{s})}{{\mathrm{d\Delta p}}_{r}}\\ \frac{{\text{dR}}_{s}({p}_{s})}{{\mathrm{d\Delta p}}_{r}}={\text{Qbh}}_{\text{sr}}{e}^{-{\text{bp}}_{r}}\end{array}\)

Évaluation de l’opérateur tangent cohérent

Il s’agit de trouver l’opérateur tangent cohérent, c’est à dire calculé à partir de la solution de \((R(Y)=0)\) en fin d’incrément. Pour une petite variation de \(R\) , en considérant cette fois \(\Delta E\) comme variable et non comme paramètre, on obtient :

\(\frac{\partial R}{\partial \Sigma }\delta \Sigma +\frac{\partial R}{\partial \Delta E}\delta \Delta E+\frac{\partial R}{\partial \Delta {\gamma}_{s}}\delta \Delta {\gamma}_{s}=0\)

Ce système peut s’écrire:

\(\frac{\partial R}{\partial Y}\delta (Y)=X\) avec \(Y=\left[\begin{array}{c}\Sigma \\ \Delta {\gamma}_{s}\end{array}\right]\) et \(X=\left[\begin{array}{c}\delta \Delta E\\ 0\end{array}\right]\)

En écrivant la matrice jacobienne sous la forme:

\(J.\delta Y=\left[\begin{array}{cc}{Y}_{0}& {Y}_{1}\\ {Y}_{2}& {Y}_{3}\end{array}\right]\left[\begin{array}{c}\Sigma \\ \Delta Z\end{array}\right]\)

Avec:

\(\begin{array}{}{Y}_{0}={\Lambda}^{-1}\\ \Delta Z=\left\lbrace \Delta {\gamma}_{s}\right\rbrace \times {n}_{s}\end{array}\)

Les sous-matrices ont pour dimensions:

\(\begin{array}{}\text{dim}({Y}_{0}={\Lambda}^{-1})=\left[6,6\right]\\ \text{dim}{Y}_{1}=\left[6,{n}_{s}\right]\\ \text{dim}{Y}_{2}=\left[{n}_{s},6\right]\\ \text{dim}{Y}_{3}=\left[{n}_{s},{n}_{s}\right]\end{array}\)

En opérant par éliminations et substitutions successives, le troisième bloc du système d’équation donne:

\(\begin{array}{c}\Delta Z=-{({Y}_{3})}^{-1}{Y}_{2}\Sigma \\ ({Y}_{0}-{Y}_{1}{({Y}_{3})}^{-1}{Y}_{2})\Sigma =\Delta E\end{array}\)

l’opérateur tangent recherchépeut donc s’écriredirectement:

\({(\frac{\partial \Sigma }{\partial E})}_{t+\Delta t}={(\frac{\partial \Sigma }{\partial \Delta E})}_{t+\Delta t}={({Y}_{0}-{Y}_{1}{Y}_{3}^{-1}{Y}_{2})}^{-1}\)