étude, le matériau ne voit pas suffisamment de cycles de charge pour pouvoir faire cette hy-
pothèse, et le(s) premier(s) cycle(s) est(sont) même d'importance majeure. Nous choisisons
donc de représenter l'effet Mullins en utilisant le modèle d'OGDEN & ROXBURGH 1999. Cela
permettra de modéliser l'effet pour des cycles à amplitude de déformation croissante :

$$
\sigma = \eta \bar{\sigma}_{eq} + \bar{\sigma}_{neq} + \hat{\sigma} \qquad (III.57)
$$

avec η la variable pilotant l'adoucissement de la contrainte :

$$
\eta = 1 - \frac{1}{r} \operatorname{erf} \left( \frac{\rho_0 \bar{\psi}^{\text{max}} - \rho_0 \bar{\psi}}{m + \beta \rho_0 \bar{\psi}^{\text{max}}} \right) \qquad (\text{III.58})
$$

où r, β et m sont des paramètres à identifier. $\bar{\psi}^{max}$ est l'énergie libre correspondant à la déformation maximale atteinte sur l'histoire des déformations. La variable d'endommagement s'applique ici sur la partie hyperélastique de la contrainte déviatorique.

Par la suite, ce cadre pourra être étendu à la partie sphérique $\hat{\sigma}$ de la contrainte. En effet, le matériau présente ce phénomène d'adoucissement sous chargement sphérique ou quasiment (cf. Figure II.18b).

III.C.3 Choix des énergies libres

III.C.3.a Partie isochore

L'énergie libre choisie (après tests d'autres formes) pour la partie isochore à l'équilibre $\bar{\psi}_{eq}$ est de type polynomiale avec $C_{10}$, $C_{01}$ et $C_{30}$ :

$$
\rho_0 \bar{\psi}_{eq} = C_{10} (I_{\bar{B}} - 3) + C_{01} (I_{\bar{B}} - 3) + C_{30} (I_{\bar{B}} - 3)^3 \quad (III.59)
$$

Pour la partie hors équilibre, par soucis de simplicité et pour ne pas introduire de paramètres superflus, une énergie libre à un paramètre, noté $C_p$, de type Neo-Hooke est utilisée :

$$
\rho_0 \bar{\psi}_{neq} = C_p (I_{\bar{B}_e} - 3) \qquad (III.60)
$$

III.C.3.b Partie volumique

La première énergie libre considérée pour la partie volumique est celle d'EIPPER 1998 :

$$
\rho_0 \hat{\psi}_{eq}(J) = \frac{\Lambda \left[ J^{\gamma} - 1 - \gamma \ln \frac{J - n_0^S}{1 - n_0^S} + \gamma n_0^S \frac{J-1}{1-n_0^S} \right]}{\gamma \left( \gamma - 1 + \frac{1}{(1-n_0^S)^2} \right)} \quad (III.61)
$$

avec le paramètre déterministe $n_0^S$, ne demandant a priori pas d'identification, car il représente
la part volumique initiale de RTV. $\gamma$ et $\Lambda$ des paramètres à identifier.
La contrainte hydrostatique s'écrit alors :

$$ \hat{\sigma}_{eq} = \frac{\Lambda \left[ J^{\gamma-1} - \frac{1}{J-n_0^S} + \frac{n_0^S}{1-n_0^S} \right]}{\left( \gamma - 1 + \frac{1}{(1-n_0^S)^2} \right)} \mathbf{1} \qquad (III.62) $$

L'énergie libre de MOERMAN et al. 2020 semble également très adaptée pour représenter la forte compressibilité des RTV/Expancel, mais aussi la non-linéarité du comportement sphérique. Elle s'écrit :

$$ \rho_0 \hat{\psi}_{eq}(J) = \kappa \left[ -(1-q)a^2 \ln \left( \cos \left( \frac{J-1}{a} \right) \right) + q b^2 \ln \left( \cosh \left( \frac{J-1}{b} \right) \right) \right] \qquad (III.63) $$

$$ \text{avec } a = \frac{2}{\pi} \begin{cases} (J_1 - 1) & \text{si } J \ge 1 \\ (J_2 - 1) & \text{si } J < 1 \end{cases}, \quad b = \frac{1}{\kappa} \begin{cases} s_1 & \text{si } J \ge 1 \\ s_2 & \text{si } J < 1 \end{cases} \text{ et } q = \begin{cases} q_1 & \text{si } J \ge 1 \\ q_2 & \text{si } J < 1 \end{cases} \qquad (III.64) $$

La contrainte hydrostatique issue de cette énergie libre s'écrit :

$$ \hat{\sigma}_{eq}(J) = \kappa \left[ (1-q)a \tan \left( \frac{J-1}{a} \right) + q b \tanh \left( \frac{J-1}{b} \right) \right] \mathbf{1} \qquad (III.65) $$

L'inconvénient est le grand nombre de paramètres à identifier (7) en raison de leur dépendance à la sollicitation. Une étude de l'influence de ces paramètres pourrait être intéressante pour voir si certains peuvent convenir à la fois en traction et en compression.

Quant à l'énergie libre pour la partie volumique hors équilibre, elle peut prendre la forme du deuxième modèle proposé par MOERMAN et al. 2020. Il s'agit d'une forme simplifiée de celle utilisée pour la partie à l'équilibre.

$$ \rho_0 \hat{\psi}_{neq}(J_e) = -\kappa_e \left[ a^2 \ln \left( \cos \left( \frac{J_e - 1}{a} \right) \right) \right] \qquad (III.66) $$

$$ \text{avec } a_e = \frac{2}{\pi} \begin{cases} (J_{1e} - 1) & \text{si } J_{1e} \ge 1 \\ (J_{2e} - 1) & \text{si } J_{2e} < 1 \end{cases} \qquad (III.67) $$

Cette énergie présente 3 paramètres mais on pourrait faire l'hypothèse que $J_{1e} = J_1$ et $J_{2e} = J_2$. Il ne resterait alors plus que le paramètre $\kappa_e$ à identifier. La contrainte en découlant s'écrit :

$$ \hat{\sigma}_{neq}(J_e) = \kappa_e \left[ a_e \tan \left( \frac{J_e - 1}{a_e} \right) \right] \mathbf{1} \qquad (III.68) $$

### III.C.4 Choix des potentiels de dissipation et lois d'évolution

La loi créée a pour but de modéliser le comportement dissipatif des matériaux R23 et R33. Ces derniers présentent également une déformation résiduelle lorsqu'ils sont sollicités en traction, compression confinée ou cisaillement. Afin de modéliser ces phénomènes, il a été choisi de construire un modèle élasto-visco-plastique à partir des équations (III.51) à (III.54). Cela est possible par le choix des potentiels de dissipation qui entrent en jeu dans les lois d'évolution
des variables $\bar{\mathbf{B}}_e$ et $J_e$.

### III.C.4.a Partie isochore

La loi d'évolution (III.53) peut être réécrite :

$$ \bar{\mathbf{D}}_a^o = \frac{\partial \bar{\varphi}_a^*}{\partial \bar{\boldsymbol{\sigma}}_{neq}^D} \qquad (III.69) $$

où $\bar{\varphi}_a^*$ est le pseudo-potentiel de dissipation dual de $\bar{\varphi}_a$ obtenu par transformation de Legendre-Fenchel (LEMAÎTRE et al. 2009) exprimé par :

$$ \bar{\varphi}_a^* = \frac{1}{2} \frac{\langle \|\bar{\boldsymbol{\sigma}}_{neq}^D\| - \bar{\chi} \rangle^2}{\bar{\nu}} \qquad (III.70) $$

avec $\bar{\nu}$ un paramètre de viscosité, $\bar{\chi}$ le seuil de plasticité, $\bar{\boldsymbol{\sigma}}_{neq}^D$ la partie déviatorique de $\bar{\boldsymbol{\sigma}}_{neq}$, $\| \cdot \|$ la norme tensorielle et $\langle \cdot \rangle$ les crochets de Mac-Cauley permettant de définir la limite d'écoulement $\langle \|\bar{\boldsymbol{\sigma}}_{neq}\| - \bar{\chi} \rangle$. Le taux de déformation isochore anélastique objectif s'écrit donc :

$$ \bar{\mathbf{D}}_a^o = \frac{\langle \|\bar{\boldsymbol{\sigma}}_{neq}^D\| - \bar{\chi} \rangle \bar{\boldsymbol{\sigma}}_{neq}^D}{\bar{\nu} \|\bar{\boldsymbol{\sigma}}_{neq}^D\|} \qquad (III.71) $$

Pour la première loi d'écoulement, on part de l'expression (III.40) de $\dot{\bar{\mathbf{B}}}_e$ et on remplace $\bar{\mathbf{D}}_a^o$ par son expression de l'équation (III.71) :

$$ \begin{align} \dot{\bar{\mathbf{B}}}_e &= -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{B}}_e + \frac{2}{3} [\mathbf{1} : \mathbf{D}_a^o] \bar{\mathbf{B}}_e + \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2\bar{\mathbf{V}}_e \cdot \mathbf{D}_a^o \cdot \bar{\mathbf{V}}_e \nonumber \\ \dot{\bar{\mathbf{B}}}_e &= -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{B}}_e + \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2\bar{\mathbf{V}}_e \cdot \bar{\mathbf{D}}_a^o \cdot \bar{\mathbf{V}}_e \tag{III.72} \\ \dot{\bar{\mathbf{B}}}_e &= -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{B}}_e + \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2\bar{\mathbf{V}}_e \cdot \left[ \frac{\langle \|\bar{\boldsymbol{\sigma}}_{neq}^D\| - \bar{\chi} \rangle}{\bar{\nu}} \frac{\bar{\boldsymbol{\sigma}}_{neq}^D}{\|\bar{\boldsymbol{\sigma}}_{neq}^D\|} \right] \nonumber \end{align} $$

La résolution de cette équation différentielle permet d'obtenir la valeur de $\bar{\mathbf{B}}_e$ à injecter dans l'expression de la contrainte déviatorique $\bar{\boldsymbol{\sigma}}$ (III.51).

Le potentiel de dissipation $\bar{\varphi}_0$ est considéré nul.

### III.C.4.b Partie volumique

On ajoute également de la visco-plasticité à la partie sphérique étant donné la déformation résiduelle présente en compression confinée, qui est mixte, mais majoritairement sphérique. De la même manière que dans la partie précédente, on transforme la loi d'évolution (III.54) par transformée de Legendre-Fenchel, et on obtient :

$$ \hat{\mathbf{D}}_a^o = \frac{\partial \hat{\varphi}_a^*}{\partial \hat{\boldsymbol{\sigma}}_{neq}^D} \qquad (III.73) $$
où $\hat{\varphi}_a^*$ est le pseudo-potentiel de dissipation dual de $\hat{\varphi}_a$ obtenu par transformation de Legendre-Fenchel exprimé par :

$$ \hat{\varphi}_a^* = \frac{1}{2} \frac{\langle \|\hat{\boldsymbol{\sigma}}_{neq} \| - \hat{\chi} \rangle^2}{\hat{\nu}} \qquad (III.74) $$

avec $\hat{\nu}$ un paramètre de viscosité, $\hat{\chi}$ le seuil de plasticité, $\| \cdot \|$ la norme tensorielle et $\langle \cdot \rangle$ les crochets de Mac-Cauley permettant de définir la limite d'écoulement $\langle \|\hat{\boldsymbol{\sigma}}_{neq} \| - \hat{\chi} \rangle$. Le taux de déformation volumique anélastique objectif s'écrit donc :

$$ \hat{\mathbf{D}}_a^o = \frac{\langle \|\hat{\boldsymbol{\sigma}}_{neq} \| - \hat{\chi} \rangle}{\hat{\nu}} \frac{\hat{\boldsymbol{\sigma}}_{neq}}{\|\hat{\boldsymbol{\sigma}}_{neq}\|} \qquad (III.75) $$

Or ce tenseur a également pour expression :

$$ \hat{\mathbf{D}}_a^o = \frac{1}{3} J_a^{-1} \dot{\mathbf{J}}_a \mathbf{1} \qquad (III.76) $$

On peut alors écrire :

$$ J_a^{-1} \dot{\mathbf{J}}_a = \frac{\langle \|\hat{\boldsymbol{\sigma}}_{neq} \| - \hat{\chi} \rangle (\hat{\boldsymbol{\sigma}}_{neq} : \mathbf{1})}{\hat{\nu} \|\hat{\boldsymbol{\sigma}}_{neq} \|} \qquad (III.77) $$

Cette équation est résolue assez rapidement grâce à un schéma d'Euler explicite par exemple. On obtient $J_a$ qu'on peut injecter dans l'expression de $\hat{\boldsymbol{\sigma}}$ (III.52), avec $J_e = J J_a^{-1}$.

Le potentiel de dissipation $\hat{\varphi}_0$ est considéré nul.

Les équations finales de ce modèle élasto-visco-plastique sont alors :

$$ \left\{
\begin{aligned}
\bar{\boldsymbol{\sigma}} &= 2J^{-1}\bar{\eta} \left[ \left(C_{10} + 3C_{30}(I_{\bar{\mathbf{B}}} - 3)^2\right) \bar{\mathbf{B}} - C_{01}\bar{\mathbf{B}}^{-1} \right] : \mathbb{P}_1 + 2J^{-1} \left[ C_p (I_{\bar{\mathbf{B}}_e} - 3) \right] : \mathbb{P}_1 \\
\hat{\boldsymbol{\sigma}} &= \hat{\eta}\kappa \left[ (1-q)a \tan\left(\frac{J-1}{a}\right) + qb \tanh\left(\frac{J-1}{b}\right) \right] \mathbf{1} + J_a^{-1}\kappa_e \left[ a \tan\left(\frac{JJ_a^{-1}-1}{a}\right) \right] \mathbf{1} \\
\dot{\bar{\mathbf{B}}}_e &= -\frac{2}{3} [\mathbf{1}:\mathbf{D}] \bar{\mathbf{B}}_e + \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2\bar{\mathbf{B}}_e \cdot \left[ \frac{\langle \|\bar{\boldsymbol{\sigma}}_{neq}^D \| - \bar{\chi} \rangle}{\bar{\nu}} \frac{\bar{\boldsymbol{\sigma}}_{neq}^D}{\|\bar{\boldsymbol{\sigma}}_{neq}^D\|} \right] \\
J_a^{-1} \dot{\mathbf{J}}_a &= \frac{\langle \|\hat{\boldsymbol{\sigma}}_{neq} \| - \hat{\chi} \rangle (\hat{\boldsymbol{\sigma}}_{neq} : \mathbf{1})}{\hat{\nu} \|\hat{\boldsymbol{\sigma}}_{neq} \|}
\end{aligned}
\right. \qquad (III.78) $$
### III.D Implémentation dans un code de calculs aux éléments finis

#### III.D.1 Description, généralités et spécificités liées à Abaqus

Les calculs aux éléments finis nécessitent la connaissance du comportement du matériau. On compte de nombreux modèles dans les codes comme l'élasticité linéaire, quelques modèles hyperélastiques ou encore l'élasto-plasticité pour les métaux par exemple. Cependant, il est souvent nécessaire de recourir à des lois spécifiques et, dans ce cas, la plupart des codes de calculs offrent la possibilité à l'utilisateur d'intégrer sa propre loi de comportement.

Abaqus est l'un des codes commerciaux les plus connus et dont l'implémentation de lois de comportement est bien documentée. Néanmoins, il existe d'autres solutions comme MFront permettant de décrire des comportements de matériaux pour plusieurs codes de calculs comme Cast3M, Abaqus, Code\_Aster et d'autres encore (HELFER et al. 2015).

L'objectif de cette thèse est de simuler le comportement des joints RTV/Expancel en utilisant Abaqus. Pour plus de simplicité, nous avons développé une subroutine pour Abaqus en langage Fortran. Pour les comportements purement hyperélastiques, il est possible de programmer une routine UHYPER. Cette dernière demande à l'utilisateur de renseigner les différentes dérivées de l'énergie libre en fonction des trois invariants de **B**. L'avantage est de pouvoir utiliser des énergies libres non natives d'Abaqus, notamment de nouvelles expressions pour la partie volumique. Abaqus ne propose que la plus basique pour les matériaux quasi-incompressibles.

**Algorithmme 1** Principaux calculs effectués dans une subroutine UMAT décrivant un comportement hyperélastique dissipatif

À l'incrément de temps $t_{n+1}$, Abaqus retourne la valeur du tenseur gradient de déformation **F**$_{n+1}$ (DFGRD1)

1. Calcul du gradient de déformation isochore

$$ \bar{\mathbf{F}}_{n+1} = J_{n+1}^{-1/3} \mathbf{F}_{n+1} $$

2. Calcul du tenseur des dilatations de Cauchy-Green gauche isochore et de sa partie déviatorique

$$ \bar{\mathbf{B}}_{n+1} = \bar{\mathbf{F}}_{n+1} \cdot \bar{\mathbf{F}}_{n+1}^T $$

$$ \bar{\mathbf{B}}_{n+1}^D = \bar{\mathbf{B}}_{n+1} - \frac{1}{3} \text{tr}(\bar{\mathbf{B}}_{n+1}) \mathbf{1} $$

3. Calcul des invariants de $\bar{\mathbf{B}}_{n+1}$ et de l'énergie de déformation déviatorique $\bar{\psi}_{n+1}^{dev}$

$$ I_{\bar{\mathbf{B}}} = \text{tr}(\bar{\mathbf{B}}_{n+1}) $$

$$ \mathbb{I}_{\bar{\mathbf{B}}} = \frac{1}{2} [I_{\bar{\mathbf{B}}}^2 - \text{tr}(\bar{\mathbf{B}}_{n+1}^2)] $$

4. Calcul de la variable d'endommagement $\eta$ de l'effet Mullins et de sa dérivée

$$ \eta_{n+1} = f(\bar{\psi}_{n+1}^{dev}, \bar{\psi}_{dev}^m) \text{ et } \left( \frac{\partial \eta}{\partial \bar{\psi}_{dev}} \right)_{n+1} = f(\bar{\psi}_{n+1}^{dev}, \bar{\psi}_{dev}^m) $$

5. Résolution de l'équation différentielle de la loi d'évolution :

Résolution de $\dot{\bar{\mathbf{B}}}_e = -\frac{2}{3}\text{tr}(\mathbf{D})\bar{\mathbf{B}}_e + \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2\bar{\mathbf{B}}_e \cdot \bar{\boldsymbol{\sigma}}_{neq} \Rightarrow$ Obtention de $(\bar{\mathbf{B}}_e)_{n+1}$

6. Calcul du tenseur des contraintes de Cauchy STRESS

$$ \boldsymbol{\sigma}_{n+1} = f(\eta_{n+1}, \bar{\mathbf{B}}_{n+1}^D, (\bar{\mathbf{B}}_e)_{n+1}, J_{n+1}) $$

7. Calcul de l'opérateur tangent eulérien DDSDDE
L'intégration de phénomènes dissipatifs ou de l'effet Mullins au comportement peut se faire grâce à la subroutine UMAT. L'utilisateur fournit l'expression de la contrainte et des variables internes nécessaires à son calcul, ainsi que la matrice tangente. Toutes ces grandeurs sont calculées à partir d'une estimation de la déformation à un incrément de temps $t_{n+1}$ fournie par Abaqus au début de l'incrément. Les étapes à suivre pour la programmation d'un modèle hyperélastique avec effet Mullins sont décrites schématiquement dans l'Algorithme 1. La partie la plus délicate à mettre en œuvre pour l'utilisateur est le calcul de l'opérateur tangent. Ce tenseur d'ordre 4 est indispensable à la convergence du calcul. Abaqus nécessite sa définition dans la configuration eulérienne ce qui rend le calcul difficile, car l'opérateur tangent s'obtient par les variations de la contrainte $\sigma$ en fonction des variations de $\mathbf{D}$.

La suite du mémoire présente la méthode de calcul analytique de l'opérateur tangent, puis une méthode numérique.

### III.D.2 Calcul exact en hyperélasticité compressible

Abaqus définit l'opérateur tangent consistant à partir de la variation objective de Jaumann de la contrainte de Kirchhoff $\tau = J\sigma$ tel que :

$$ \delta(J\sigma) = J(\mathbb{C} : \delta\mathbf{D} + \delta\mathbf{W} \cdot \boldsymbol{\sigma} - \boldsymbol{\sigma} \cdot \delta\mathbf{W}) \quad (III.79) $$

Calculer $\mathbb{C}$ revient à exprimer les variations de $J\sigma$ grâce à la différentiation :

$$ \delta(J\sigma) = \frac{\partial J\sigma}{\partial \mathbf{F}} : d\mathbf{F} \quad (III.80) $$

Le développement de cette expression peut être remis en forme tel que le second terme de (III.79) afin d'identifier l'expression de l'opérateur tangent $\mathbb{C}$. Cette étape peut s'avérer fastidieuse, car elle nécessite beaucoup de développements.

Afin de simplifier la démarche, la documentation d'Abaqus (ABAQUS 2020b) propose les expressions des variations indispensables au calcul de $\mathbb{C} : \delta\bar{\mathbf{B}}$, $\delta\bar{\mathbf{B}}^2$, $\delta\mathbf{I}_{\bar{\mathbf{B}}}$, $\delta\mathbf{I}_{\bar{\mathbf{B}}}^2$ et $\delta J$.
On définit les variations des dérivées d'énergies de déformation par :

$$ \delta\left(\frac{\partial\bar{\psi}}{\partial\mathbf{I}_{\bar{\mathbf{B}}}}\right) = \frac{\partial^2\bar{\psi}}{\partial\mathbf{I}_{\bar{\mathbf{B}}}^2}\delta\mathbf{I}_{\bar{\mathbf{B}}} + \frac{\partial^2\bar{\psi}}{\partial\mathbf{I}_{\bar{\mathbf{B}}}\partial\mathbf{I}_{\bar{\mathbf{B}}}}\delta\mathbf{I}_{\bar{\mathbf{B}}} \quad (III.81) $$

$$ \delta\left(\frac{\partial\bar{\psi}}{\partial\mathbf{I}_{\bar{\mathbf{B}}}}\right) = \frac{\partial^2\bar{\psi}}{\partial\mathbf{I}_{\bar{\mathbf{B}}}\partial\mathbf{I}_{\bar{\mathbf{B}}}}\delta\mathbf{I}_{\bar{\mathbf{B}}} + \frac{\partial^2\bar{\psi}}{\partial\mathbf{I}_{\bar{\mathbf{B}}}^2}\delta\mathbf{I}_{\bar{\mathbf{B}}} \quad (III.82) $$

En se servant de ces expressions, il ne reste plus qu'à calculer celle de $J\sigma$ :

$$ \delta(J\sigma) = \delta\left(2\rho_0 \left[\frac{\partial\bar{\psi}}{\partial\mathbf{I}_{\bar{\mathbf{B}}}}\bar{\mathbf{B}} + \frac{\partial\bar{\psi}}{\partial\mathbf{I}_{\bar{\mathbf{B}}}}(\mathbf{I}_{\bar{\mathbf{B}}}\bar{\mathbf{B}} - \bar{\mathbf{B}}^2)\right] : \mathbb{P}_1 + \rho_0 J \frac{\partial\bar{\psi}}{\partial J}\mathbf{1}\right) \quad (III.83) $$
Le développement de cette expression permet d'écrire le produit doublement contracté d'un tenseur d'ordre 4 et des variations de **D**. Cela permet donc de définir :

$$
\begin{align}
\mathcal{C} = \bar{\mathcal{C}} + \hat{\mathcal{C}} &= \frac{2}{J} \left[ \frac{\partial \bar{\psi}}{\partial \mathbf{I}_{\bar{\mathbf{B}}}} + \mathbf{I}_{\bar{\mathbf{B}}} \frac{\partial \bar{\psi}}{\partial \mathbf{\Pi}_{\bar{\mathbf{B}}}} \right] \left[ \mathbb{H}_1 - \frac{1}{3} (\mathbb{H}_1 : \mathbf{1}) \otimes \mathbf{1} \right] - \frac{2}{J} \frac{\partial \bar{\psi}}{\partial \mathbf{\Pi}_{\bar{\mathbf{B}}}} \left[ \mathbb{H}_2 - \frac{1}{3} (\mathbb{H}_2 : \mathbf{1}) \otimes \mathbf{1} \right] \nonumber \\
&\quad + \frac{4}{J} \left[ \frac{\partial^2 \bar{\psi}}{\partial \mathbf{I}_{\bar{\mathbf{B}}}^2} + \frac{\partial \bar{\psi}}{\partial \mathbf{\Pi}_{\bar{\mathbf{B}}}} + 2 \mathbf{I}_{\bar{\mathbf{B}}} \frac{\partial^2 \bar{\psi}}{\partial \mathbf{I}_{\bar{\mathbf{B}}} \partial \mathbf{\Pi}_{\bar{\mathbf{B}}}} + \mathbf{I}_{\bar{\mathbf{B}}}^2 \frac{\partial^2 \bar{\psi}}{\partial \mathbf{\Pi}_{\bar{\mathbf{B}}}^2} \right] \left[ \bar{\mathbf{B}}^D \otimes \bar{\mathbf{B}}^D \right] + \frac{4}{J} \frac{\partial^2 \bar{\psi}}{\partial \mathbf{\Pi}_{\bar{\mathbf{B}}}^2} \left[ \left(\bar{\mathbf{B}}^2\right)^D \otimes \left(\bar{\mathbf{B}}^2\right)^D \right] \nonumber \\
&\quad - \frac{4}{J} \left[ \frac{\partial^2 \bar{\psi}}{\partial \mathbf{I}_{\bar{\mathbf{B}}} \partial \mathbf{\Pi}_{\bar{\mathbf{B}}}} + \mathbf{I}_{\bar{\mathbf{B}}} \frac{\partial^2 \bar{\psi}}{\partial \mathbf{\Pi}_{\bar{\mathbf{B}}}^2} \right] \left[ \bar{\mathbf{B}}^D \otimes \left(\bar{\mathbf{B}}^2\right)^D + \left(\bar{\mathbf{B}}^2\right)^D \otimes \bar{\mathbf{B}}^D \right] \nonumber \\
&\quad - \frac{4}{3J} \frac{\partial \bar{\psi}}{\partial \mathbf{I}_{\bar{\mathbf{B}}}} \left[ \mathbf{1} \otimes \bar{\mathbf{B}}^D \right] + \frac{4}{3J} \frac{\partial \bar{\psi}}{\partial \mathbf{\Pi}_{\bar{\mathbf{B}}}} \left[ \mathbf{1} \otimes \left(\bar{\mathbf{B}}^2\right)^D \right] + \left[ J \frac{\partial^2 \hat{\psi}}{\partial J^2} + \frac{\partial \hat{\psi}}{\partial J} \right] \left[ \mathbf{1} \otimes \mathbf{1} \right] \tag{III.84}
\end{align}
$$

avec

$$ (\mathbb{H}_1)_{ijkl} = \frac{1}{2} (\delta_{ik}\bar{B}_{jl} + \bar{B}_{ik}\delta_{jl} + \delta_{il}\bar{B}_{jk} + \bar{B}_{il}\delta_{jk}) \quad (\text{III.85}) $$

$$ (\mathbb{H}_2)_{ijkl} = \frac{1}{2} (\delta_{ik}\bar{B}_{jp}\bar{B}_{pl} + \bar{B}_{ip}\bar{B}_{pk}\delta_{jl} + \delta_{il}\bar{B}_{jp}\bar{B}_{pk} + \bar{B}_{ip}\bar{B}_{pl}\delta_{jk}) + \bar{B}_{ik}\bar{B}_{jl} + \bar{B}_{il}\bar{B}_{jk} \quad (\text{III.86}) $$

**Partie issue de l'effet Mullins** Dans le cas du modèle d'Ogden et Roxburgh, Abaqus définit l'opérateur tangent comme étant :

$$ \mathcal{C} = \eta \bar{\mathcal{C}} + J \frac{\partial \eta}{\partial \bar{\psi}} \tilde{\bar{\sigma}} \otimes \tilde{\bar{\sigma}} + \hat{\mathcal{C}} \quad (\text{III.87}) $$

avec $\tilde{\bar{\sigma}}$ la contrainte déviatorique vierge, i.e. le premier cycle de chargement, donc sans en-dommagement. $\bar{\mathcal{C}}$ et $\hat{\mathcal{C}}$ sont respectivement les parties déviatorique et sphérique de l'opérateur tangent. La dérivée de $\eta$ par rapport à l'énergie de déformation s'écrit :

$$ \frac{\partial \eta}{\partial \bar{\psi}} = - \frac{2}{r\sqrt{\pi}(m + \beta \rho_0 \bar{\psi}^{\text{max}})} \exp \left( \frac{\rho_0 \bar{\psi}^{\text{max}} - \rho_0 \bar{\psi}}{m + \beta \rho_0 \bar{\psi}^{\text{max}}} \right)^2 \quad (\text{III.88}) $$

### III.D.3 Calcul numérique de l'opérateur tangent

La méthode de calcul numérique de l'opérateur tangent est basée sur des perturbations de **F**. Elle a notamment été introduite par MIEHE 1995 et on la retrouve sous plusieurs formes dans les travaux de SUN et al. 2008 et ZHU et al. 2022.

On considère une forme linéarisée des variations de Jaumann de la contrainte de Kirchhoff $\tau$ :

$$ \Delta\tau - \Delta W\tau - \tau\Delta W^T = \mathcal{C}^\tau : \Delta D \quad (\text{III.89}) $$

La différence $\Delta D$ s'exprime en fonction de $\Delta F$ et c'est ce dernier qui va subir l'opération de perturbation afin d'obtenir :

$$ \Delta \mathbf{D}^{(ij)} = \frac{1}{2} \left[ \Delta \mathbf{F}^{(ij)} \mathbf{F}^{-1} + \left( \Delta \mathbf{F}^{(ij)} \mathbf{F}^{-1} \right)^T \right] \quad (\text{III.90}) $$
avec $\Delta \mathbf{F}^{(ij)}$ les perturbations de $\mathbf{F}$ :

$$ \Delta \mathbf{F}^{(ij)} = \frac{\varepsilon}{2} [\mathbf{e}_i \otimes \mathbf{e}_j \cdot \mathbf{F} + \mathbf{e}_j \otimes \mathbf{e}_i \cdot \mathbf{F}] \quad \text{avec } i, j = 1, 2, 3 \qquad (\text{III.91}) $$

$\mathbf{e}_i \otimes \mathbf{e}_j$ sont les tenseurs de la base canonique des tenseurs d'ordre 2 et $\varepsilon$ est un scalaire définissant la perturbation.

Cela donne en tout six perturbations du gradient de déformation qui, lorsqu'elles sont insérées dans l'équation (III.89), donnent l'expression :

$$ \Delta \boldsymbol{\tau} = \mathcal{C}^{\tau(ij)} : \frac{\varepsilon}{2} (\mathbf{e}_i \otimes \mathbf{e}_j + \mathbf{e}_j \otimes \mathbf{e}_i) \qquad (\text{III.92}) $$

et le gradient perturbé s'écrit :

$$ \mathbf{F}^{(ij)} = \mathbf{F} + \Delta \mathbf{F}^{(ij)} \qquad (\text{III.93}) $$

$\Delta \boldsymbol{\tau}$ peut être estimée par les différences finies comme $\Delta \boldsymbol{\tau} \approx \boldsymbol{\tau}(\mathbf{F}^{(ij)}) - \boldsymbol{\tau}(\mathbf{F})$.

Enfin, après des considérations de symétrie, l'opérateur tangent peut s'écrire de manière indicielle :

$$ \mathcal{C}_{ijkl}^{\tau} \approx \frac{1}{\varepsilon} [J^{(kl)} \sigma_{ij}(\mathbf{F}^{(kl)}) - J \sigma_{ij}(\mathbf{F})] \qquad (\text{III.94}) $$

L'expression respectant le formalisme d'Abaqus est :

$$ \mathcal{C}^{abq} = \frac{1}{J} \mathcal{C}^{\tau} \qquad (\text{III.95}) $$

Plus pratiquement, l'implémentation dans la UMAT se fait comme dans l'Algorithme 2. Le paramètre $\varepsilon$ varie en général entre $10^{-8}$ et $10^{-4}$.

Dans le cas d'un modèle dissipatif, effectuer les variations de $\mathbf{F}$ revient à résoudre les équations différentielles des lois d'évolution en considérant chacune des six perturbations.
**Algorithmme 2** Implémentation d'un opérateur tangent numérique dans une subroutine UMAT d'un comportement hyperélastique dissipatif

À l'incrément de temps $t_{n+1}$, Abaqus retourne la valeur du tenseur gradient de déformation $\mathbf{F}_{n+1}$. Pour simplifier la notation, on le notera simplement $\mathbf{F}$ et les grandeurs associées $\bar{\mathbf{B}}$, $\bar{\mathbf{B}}^D$, $\bar{\mathbf{B}}_e$, etc.

1. Calcul du gradient de déformation isochore

$$ \bar{\mathbf{F}} = J^{-1/3} \mathbf{F} $$

2. Calcul du tenseur des dilatations de Cauchy-Green gauche isochore et de sa partie déviationique

$$ \begin{aligned} \bar{\mathbf{B}} &= \bar{\mathbf{F}} \cdot \bar{\mathbf{F}}^T \\ \bar{\mathbf{B}}^D &= \bar{\mathbf{B}} - \frac{1}{3} \mathrm{tr}(\bar{\mathbf{B}}) \mathbf{1} \end{aligned} $$

3. Calcul des invariants de $\bar{\mathbf{B}}$ et de l'énergie de déformation déviatorique $\bar{\psi}^D$

$$ \begin{aligned} I_{\bar{\mathbf{B}}} &= \mathrm{tr}(\bar{\mathbf{B}}) \\ \mathbf{I}_{\bar{\mathbf{B}}} &= \frac{1}{2} [I_{\bar{\mathbf{B}}}^2 - \mathrm{tr}(\bar{\mathbf{B}}^2)] \end{aligned} $$

4. Calcul de la variable d'endommagement $\eta$ de l'effet Mullins

$$ \eta = f(\bar{\psi}^D, \bar{\psi}_D^m) $$

5. Résolution de l'équation différentielle de la loi d'évolution

$$ \text{Résolution de } \dot{\bar{\mathbf{B}}}_e = f(\bar{\mathbf{B}}_e) \Rightarrow \text{Obtention de } (\bar{\mathbf{B}}_e) $$

6. Calcul du tenseur des contraintes de Cauchy STRESS

$$ \boldsymbol{\sigma} = f(\eta, \bar{\mathbf{B}}^D, \bar{\mathbf{B}}_e, J) $$

7. Calcul des 6 gradients de déformation perturbés $\mathbf{F}^{(kl)}$ et leur déterminant $J^{(kl)}$

$$ \mathbf{F}^{(kl)} = \mathbf{F} + \frac{\varepsilon}{2} [\mathbf{e}_k \otimes \mathbf{e}_l \mathbf{F} + \mathbf{e}_l \otimes \mathbf{e}_k \mathbf{F}] $$

8. Calcul des 6 tenseurs des dilatations de Cauchy-Green gauche isochore

$$ \bar{\mathbf{B}}^{(kl)} = \bar{\mathbf{F}}^{(kl)} \cdot \bar{\mathbf{F}}^{(kl)T} $$

9. Calcul des invariants de $\bar{\mathbf{B}}^{(kl)}$ et de l'énergie libre $\bar{\psi}^{(kl)}$

10. Calcul de la variable d'endommagement $\eta^{(kl)}$ de l'effet Mullins

11. Résolution de l'équation différentielle de la loi d'évolution pour les 6 perturbations

$$ \text{Résolution de } \dot{\bar{\mathbf{B}}}_e^{(kl)} = f(\bar{\mathbf{B}}_e^{(kl)}) \Rightarrow \text{Obtention de } \bar{\mathbf{B}}_e^{(kl)} $$

12. Calcul de la contrainte $\boldsymbol{\sigma}^{(kl)}$

$$ \boldsymbol{\sigma}^{(kl)} = f(\bar{\mathbf{B}}^{(kl)}, \det(\mathbf{F}^{(kl)}), \bar{\mathbf{B}}_e^{(kl)}, \eta^{(kl)}) $$

13. Calcul de l'opérateur tangent $\mathcal{C}_{ijkl}$

$$ \mathcal{C}_{ijkl} = \frac{1}{J\varepsilon} [J^{(kl)} \sigma_{ij}^{(kl)}(\mathbf{F}^{(kl)}) - J \sigma_{ij}(\mathbf{F})] $$
# Conclusions du chapitre

Ce chapitre a permis de détailler le formalisme adopté pour la modélisation des matériaux RTV/microsphères qui sont fortement compressibles. Les différentes étapes à retenir sont :

* La compressibilité est intégrée, par l'hypothèse d'une décomposition isochore/volumique du gradient de déformation et par l'utilisation d'énergies libres volumiques adaptées pour la modélisation du comportement non-linéaire de mousses d'élastomères. Il s'agit d'une hypothèse forte de ce travail.

* l'introduction de configurations intermédiaires pour représenter les variables internes permet de modéliser plusieurs phénomènes dissipatifs comme la visco-élasticité, la plasticité ou une combinaison de ces phénomènes,

* les configurations intermédiaires sont introduites sur la partie isochore, mais aussi sur la partie volumique, permettant la modélisation de phénomènes dissipatifs qui seraient liés aux microsphères,

* ce formalisme induit un découplage de la contrainte de Cauchy en contributions déviato-rique et sphérique couplée à une décomposition équilibre (hyperélastique) / hors équilibre (hyperélasto-dissipative),

* l'effet Mullins est modélisé à la fois sur la partie isochore et volumique de la transformation, car l'endommagement est présent sur des essais de traction et cisaillement simple, mais aussi sur des essais de compression confinée qui a une part majoritaire sphérique,

* l'implémentation de ce modèle se fait par le biais d'une subroutine UMAT dans le code de calculs Abaqus,

* l'opérateur tangent peut être calculé de manière analytique pour la partie hyperélastique du modèle. La partie dissipative est plus complexe à calculer. C'est pourquoi un opérateur tangent numérique est utilisé, permettant une plus grande souplesse de développement (au détriment d'un temps de calcul un peu plus important).

Globalement les avantages de ce modèle sont de pouvoir intégrer différents types de phénomènes dissipatifs dans un formalisme étendu aux grandes compressibilités, et ce pour tous types de sollicitations.
# Chapitre IV

## Stratégie d'identification de modèles hyperélastiques compressibles

Ce chapitre introduit la méthode d'identification utilisée dans la thèse. Cette dernière assez innovante est basée sur la décomposition orthogonale aux valeurs propres (POD) couplée à une interpolation (dans notre cas krigeage) qui permet une réduction du modèle. Les données sont alors représentées selon leur influence pour la sollicitation mécanique considérée. Donc en plus d'identifier la valeur des paramètres, il est possible d'estimer leur identifiabilité. En effet, un paramètre n'influençant pas le comportement a très peu de chance d'être identifié. C'est l'avantage de cette méthode par rapport aux méthodes inverses, qui permettraient d'identifier tous les paramètres, mais sans garantie sur l'unicité de la solution.

Cette méthode a d'abord été utilisée pour identifier les paramètres du modèle hyperélastique compressible. Cela a pour avantage de mettre en place la chaîne d'identification pour des équations plus simples, ne nécessitant pas l'intégration d'une ou plusieurs lois d'évolution.

Grâce au découplage de la contrainte, il est possible d'identifier les paramètres en plusieurs temps sur des sollicitations différentes comme la compression hydrostatique, la traction uniaxiale et la compression confinée. Les paramètres permettant de représenter l'effet Mullins seront également identifiés.

La stabilité du modèle avec ces paramètres (hors ceux décrivant l'effet Mullins) est ensuite vérifiée pour s'assurer que les simulations éléments finis ne rencontrent pas de problèmes.
## PLAN DU CHAPITRE IV

<table><tr><td>IV.A Méthode d'identification par réduction de modèle</td><td>95</td></tr><tr><td>IV.A.1 Préparation des données d'entrée</td><td>95</td></tr><tr><td>IV.A.2 Réduction des données d'entrée</td><td>96</td></tr><tr><td>IV.A.3 Interpolation des snapshots dans l'espace α et identification des paramètres</td><td>97</td></tr><tr><td>IV.B Application de la méthode</td><td>98</td></tr><tr><td>IV.B.1 Stratégie adoptée</td><td>98</td></tr><tr><td>IV.B.2 Identification sur la compression hydrostatique</td><td>99</td></tr><tr><td>IV.B.3 Identification sur des chargements mixtes déviatorique/sphérique</td><td>107</td></tr><tr><td>IV.B.4 Identification de la partie liée à l'effet Mullins</td><td>111</td></tr><tr><td>IV.B.4.a Partie déviatorique</td><td>112</td></tr><tr><td>IV.B.4.b Partie sphérique</td><td>114</td></tr><tr><td>IV.C Stabilité du modèle</td><td>116</td></tr><tr><td>Conclusions du chapitre</td><td>121</td></tr></table>
# IV.A Méthode d'identification par réduction de modèle

La méthode d'identification employée est basée sur une technique de décomposition aux valeurs propres, plus communément appelée POD (Proper Orthogonal Decomposition) (EZZAAMARI et al. 2021; LE QUILLIEC et al. 2015). Elle consiste à construire un sous-espace de dimension réduite où les données de l'espace initial sont projetées. La direction des vecteurs de cette base orthogonale est déterminée par la variance des données d'entrée. Cela permet donc de garder le maximum d'informations par ordre d'importance de la variance sur les données tout en ayant une dimensionnalité réduite. Nous appliquons cette méthode à l'identification des paramètres. Les données d'entrée sont, dans notre cas, des courbes contrainte-déformation générées numériquement à partir du modèle créé dans le chapitre précédent.

## IV.A.1 Préparation des données d'entrée

Les courbes contrainte-déformation sont obtenues à partir d'un plan d'expériences réalisé sur les paramètres recherchés du modèle. Ce dernier fait varier les valeurs des paramètres sur une plage de valeur estimée. Deux outils d'échantillonnage ont été utilisés à cet effet : le LHS (Latin Hypercube Sampling) et le plan factoriel complet (PFC). Le LHS est une méthode d'échantillonnage quasi-aléatoire. Chaque plage de valeurs est divisée en intervalles de probabilité égale où un point y est placé aléatoirement. Le PFC permet de créer un plan d'expériences combinant toutes les combinaisons possibles de paramètres. Ainsi, en divisant chaque plage de valeurs en n intervalles, on obtient un plan d'expérience de $n^m$ points, où m est le nombre de paramètres.

<figure><img src="image_45.png" /><figcaption>FIGURE IV.1 – Représentation des couples de paramètres $C_{10}$ et $C_{01}$ obtenus par LHS (a) et plan factoriel complet (b).</figcaption></figure>

Prenons l'exemple d'un modèle de Mooney-Rivlin incompressible : les plans d'expériences de 64 points des paramètres $C_{10}$ et $C_{01}$ sont présentés sur la Figure IV.1 où l'échelle de couleurs correspond aux variations de $C_{10}$.

En traction, la réponse monodimensionnelle (i.e. la contrainte de Cauchy) s'écrit :

$$ \sigma = 2(C_{10} + \frac{C_{01}}{\lambda})(\lambda^2 - \frac{1}{\lambda}) \qquad (\text{IV.1}) $$

La contrainte peut ensuite être calculée pour chaque point du plan d'expériences. L'évolution
de chacune des contraintes est sauvegardée à intervalles réguliers d'élongation dans un vecteur. On obtient alors $n$ vecteurs contrainte correspondant aux données d'entrée. Ces derniers sont appelés plus couramment des *snapshots* et sont échantillonnés en un nombre de points qu'on notera $m$. Ils sont représentés sur la Figure IV.2.

<figure><img src="image_46.png" /><figcaption>FIGURE IV.2 – Courbes contrainte-élongation calculées avec un modèle de Mooney-Rivlin dont les paramètres correspondent au LHS présenté sur la Figure IV.1a.</figcaption></figure>

## IV.A.2 Réduction des données d'entrée

Á partir de cette étape, on utilise la POD, qui va permettre de représenter les snapshots, notés $\varphi_i$, dans une base de dimension réduite, notée $\Psi$.

On commence par créer $\mathbf{S}$ ($m \times n$) la matrice des snapshots centrés autour du snapshot moyen $\varphi_0$.

Á partir de $\mathbf{S}$, on calcule $\mathbf{C}$ ($n \times n$), la matrice de covariance, qui a pour expression :

$$ \mathbf{C} = \mathbf{S}^T \mathbf{S} \qquad (\text{IV.2}) $$

La base est définie par l'expression :

$$ \mathbf{\Psi} = \mathbf{S} \mathbf{\phi} \mathbf{A}^{-1} \qquad (\text{IV.3}) $$

avec $\mathbf{\phi}$ les vecteurs propres de $\mathbf{C}$ et $\mathbf{A}$ la matrice diagonale des valeurs propres de $\mathbf{C}$. La base est orthogonale et ses directions sont classées par ordre d'importance.

Dans cette base, les snapshots s'écrivent :

$$ \varphi_i = \varphi_0 + \mathbf{\Psi} \alpha \qquad (\text{IV.4}) $$

avec $\alpha = \mathbf{A} \mathbf{\phi}^T (n \times n)$ leurs coordonnées classées de manière décroissante.

En utilisant les trois premières dimensions de cette base, on obtient un nuage de $n$ points dans le repère $(0; \alpha_1, \alpha_2, \alpha_3)$. Il est représenté sur la Figure IV.3. On remarque que la variance est portée seulement sur les deux premières directions de la base réduite. En effet, l'échelle de valeurs de $\alpha_3$ est très proche de zéro.
<figure><img src="image_47.png" /><figcaption>FIGURE IV.3 – Données représentées dans l'espace $\alpha$ de dimension réduite créé par POD.</figcaption></figure>

### IV.A.3 Interpolation des snapshots dans l'espace $\alpha$ et identification des paramètres

L'identification des paramètres se fait grâce à une méthode d'interpolation. Ici on choisit le krigeage qui consiste à créer une surface de réponse (continue) à partir du nuage de points (discret) résultant du plan d'expériences de départ.

Le krigeage donne accès à une surface de réponse dont la dimension donne le nombre de paramètres identifiables. Cette dernière s'appelle une variété de formes (ou *shape manifold* en anglais). Sur la base des 3 premiers $\alpha$, la dimension de cette variété est visuellement plus claire à déterminer. Elle peut être construite à partir d'une grille régulière très fine des plages de valeurs du plan d'expériences (si elle est 0D, 1D ou 2D). Les coordonnées $\alpha_i$ des points du maillage sont prédites à partir du méta-modèle de krigeage et permettent la visualisation présentée sur la Figure IV.4a. Le krigeage a été réalisé avec le package Python « pyKriging » (PAULSON & RAGKOUSIS 2015).

Si la dimension de la variété est égale à 3, on peut construire un volume à partir d'une succession d'hypersurfaces représentées à intervalles réguliers. Pour des dimensions supérieures à trois, la méthode graphique devient inutilisable. Il existe cependant des algorithmes capables de déterminer la dimensionnalité locale, celui de FUKUNAGA & OLSEN 1971 par exemple.

Une courbe de contrainte-elongation expérimentale peut être projetée dans la base réduite. On obtient alors un nouveau point dans le nuage. Afin d'identifier les valeurs de paramètres correspondant aux coordonnées $\alpha_i$ de ce point, un algorithme d'optimisation minimise l'erreur entre les coordonnées $\alpha_i$ du point expérimental et celles prédites par le modèle de krigeage pour des valeurs de $C_{10}$ et $C_{01}$. Si les paramètres sont identifiables, on doit pouvoir trouver les valeurs donnant une courbe représentant le mieux le jeu de données expérimental. Dans cet exemple, la manifold a deux dimensions. Les deux paramètres seraient donc identifiables.

Pour l'exemple du modèle de Mooney-Rivlin, on construit un snapshot expérimental pour des valeurs $C_{10} = C_{01} = 0,5$ MPa. La projection de ce dernier dans la base des $\alpha$ est présentée sur
<figure><img src="image_48.png" /><figcaption>FIGURE IV.4 – Données interpolées par krigeage représentées dans l'espace α de dimension réduite créé par POD. (b) Le snapshot expérimental est ajouté en rouge dans l'espace réduit.</figcaption></figure>

la Figure IV.4b par le point de couleur rouge. L'application de la méthode permet d'identifier les valeurs $C_{10} = 0,49997843$ MPa et $C_{01} = 0,50003554$ MPa.

Les principales étapes de la méthode d'identification sont résumées sur le schéma de la Figure IV.5.

# IV.B Application de la méthode

## IV.B.1 Stratégie adoptée

La méthode d'identification a été développée et appliquée dans un premier temps sur une partie des paramètres de la loi de comportement hyperélastique compressible. La contrainte dans ce cas simplifié s'écrit comme suit :

$$ \sigma = \bar{\sigma} + \hat{\sigma} = 2\rho_0 J^{-1} \bar{\eta} \left[ \frac{\partial \bar{\psi}}{\partial I_{\bar{B}}} \bar{B} - \frac{\partial \bar{\psi}}{\partial \Pi_{\bar{B}}} \bar{B}^{-1} \right] : \mathbb{P}_1 + \rho_0 \hat{\eta} \frac{\partial \hat{\psi}}{\partial J} \mathbf{1} \quad (IV.5) $$

Comme évoqué dans le chapitre précédent, la contrainte de Cauchy est découplée en une partie déviatorique et une partie sphérique.

Partant de cette observation, nous avons fait l'hypothèse que les paramètres de chaque contribution pouvaient être identifiés sur des essais correspondants.

Dans l'ordre, nous allons donc identifier :

- les paramètres de $\hat{\sigma}$ avec la compression hydrostatique pour les matériaux R23 et R33,

- les paramètres de $\bar{\sigma}$ avec la traction uniaxiale et la compression confinée pour les matériaux R23 et R33 en utilisant les paramètres identifiés en compression hydrostatique pour la partie sphérique.
<figure><img src="image_49.png" /><figcaption>FIGURE IV.5 – Procédure d'identification des paramètres du modèle.</figcaption></figure>

Dans un premier temps, on néglige les variables de dommage de l'effet Mullins $\bar{\eta}$ et $\hat{\eta}$ en identifiant les paramètres sur des courbes monotones. Par la suite, nous reviendrons sur l'identification des paramètres décrivant leur évolution.

## IV.B.2 Identification sur la compression hydrostatique

Dans le cas de la partie sphérique, les données d'entrée de la POD ont été calculées analytiquement grâce à l'expression de la pression hydrostatique :

$$P = -\frac{1}{3}\text{tr}\left(\rho_0 \frac{\partial \hat{\psi}}{\partial J} \mathbf{1}\right) \qquad (\text{IV.6})$$

Les courbes expérimentales sur lesquelles la pression est identifiée disposent, de par la méthode de porosimétrie à mercure, seulement d'une trentaine de points de mesure. Afin d'améliorer l'identification, les courbes ont été interpolées pour permettre un ré-échantillonnage avec 100 points régulièrement espacés.
Modèle de Eipper (EIPPER 1998)

On considère, dans un premier temps, l'énergie libre d'Eipper :

$$
\rho_0 \hat{\psi}(J) = \frac{\Lambda \left[ J^{\gamma} - 1 - \gamma \ln \frac{J - n_0^S}{1 - n_0^S} + \gamma n_0^S \frac{J-1}{1-n_0^S} \right]}{\gamma \left( \gamma - 1 + \frac{1}{(1-n_0^S)^2} \right)} \quad (IV.7)
$$

La pression hydrostatique, dépendant de 2 paramètres Λ et γ, s'écrit :

$$
P = - \frac{\Lambda}{\left( \gamma - 1 + \frac{1}{(1 - n_0^S)^2} \right)} \left[ J^{(\gamma-1)} - \frac{1}{J - n_0^S} + \frac{n_0^S}{1 - n_0^S} \right] \quad (IV.8)
$$

On se place dans le cas du matériau R23 donc $n_0^S = 0,77$. Un plan d'expérience de 75 points est généré avec un LHS. La Figure IV.6a présente les jeux de paramètres pour $\Lambda$ et $\gamma$. Les snapshots obtenus par calcul avec l'équation (IV.8) sont représentés sur la Figure IV.6b (l'échelle de couleur représente l'évolution du paramètre $\gamma$). Visuellement, la manifold obtenue a deux dimensions. Cela signifie que les deux paramètres pourraient être identifiables. Cependant, la Figure IV.6c montre que la projection dans l'espace des $\alpha$ du vecteur pression expérimental (le point noir) est très éloignée de celle des snapshots du plan d'expérience. La conséquence est une incertitude lors de l'identification conduisant même à des valeurs erronées. La distance élevée entre le point noir et la manifold est en général attribuée au fait que la loi de comportement n'est pas apte à représenter les données expérimentales. Une autre raison pourrait être un domaine de variation des paramètres du plan d'expériences pas suffisamment étendu. Le graphique IV.6d montre les courbes de pression-variation de volume expérimentale et identifiée. On remarque que le modèle décrit très bien la position de l'asymptote verticale, mais ne semble pas apte à modéliser la non-linéarité du pied de courbe. Dans le cas de notre étude, la bonne représentation du comportement sur les premiers 20% de variation de volume est un objectif majeur (plus que la bonne représentation du comportement asymptotique). Afin d'atteindre cet objectif, il faut trouver un modèle de comportement avec plus de non-linéarités.

**Modèle de Moerman** Le modèle de MOERMAN *et al.* 2020 est spécialement développé pour décrire le comportement hydrostatique de matériaux poreux. Il décrit en effet les trois phénomènes typiques de la compression des mousses et autres matériaux cellulaires, à savoir : la déformation élastique des pores, leur flambement et enfin la densification du matériau. L'énergie libre de ce modèle s'écrit :

$$
\rho_0 \hat{\psi}(J) = \kappa \left[ -(1-q)a^2 \ln \left( \cos \left( \frac{J-1}{a} \right) \right) + q b^2 \ln \left( \cosh \left( \frac{J-1}{b} \right) \right) \right] \quad (IV.9)
$$

$$
\text{avec } a = \frac{2}{\pi} \begin{cases} (J_1 - 1) & \text{si } J \ge 1 \\ (J_2 - 1) & \text{si } J < 1 \end{cases}, \quad b = \frac{1}{\kappa} \begin{cases} s_1 & \text{si } J \ge 1 \\ s_2 & \text{si } J < 1 \end{cases} \text{ et } q = \begin{cases} q_1 & \text{si } J \ge 1 \\ q_2 & \text{si } J < 1 \end{cases} \quad (\text{IV.10})
$$

La pression hydrostatique a pour expression :
<figure><img src="image_50.png" /><figcaption>FIGURE IV.6 – Étapes d'identification des paramètres de la pression hydrostatique sur le matériau R23. Valeurs identifiées : $\Lambda=3,87$ MPa et $\gamma=18,49$.</figcaption></figure>

$$P = -\kappa \left[ (1-q)a \tan\left(\frac{J-1}{a}\right) + qb \tanh\left(\frac{J-1}{b}\right) \right] \quad (IV.11)$$

Le paramètre $J_2$ fixe la position de l'asymptote verticale. Il correspond à $n_0^S$ dans le modèle de Eipper. Sa valeur est donc connue et est fixée à 0,77 pour le R23 et 0,67 pour le R33. Il reste à identifier $\kappa$, $s_2$ et $q_2$. Le paramètre $\kappa$ correspond au module en pied de courbe (i.e. le module de compressibilité du matériau) tel que représenté sur la Figure IV.7. Le paramètre $s_2$ fixe la position du premier point d'inflexion et $q_2$ l'inclinaison du plateau.

### REMARQUE IV.1

Dans le cas de l'étude de la compression hydrostatique, il sera bien évidemment impossible d'identifier $J_1$, $s_1$ et $q_1$. Nous reviendrons sur ces paramètres ultérieurement.

**Identification des trois paramètres** - On tente d'identifier, dans un premier temps, ces paramètres pour chaque matériau dans la gamme complète de déformation, c'est-à-dire jusqu'à
<figure><img src="image_51.png" /><figcaption>FIGURE IV.7 – Influence des paramètres du modèle de Moerman MOERMAN et al. 2020</figcaption></figure>

23% pour le R23 et 33% pour le R33. Le plan d'expériences à trois dimensions est construit à partir d'un LHS. On fait varier $s_2$ et $q_2$ entre 0 et 1 et $\kappa$ entre 60 et 100 MPa.

Après avoir créé les snapshots et appliqué la POD, on peut les représenter dans l'espace des $\alpha$ et effectuer le krigeage. Les résultats pour le R23 sont représentés sur le bas de la Figure IV.8.

<figure><img src="image_52.png" /><figcaption>FIGURE IV.8 – Identification des paramètres pour le R23 : représentation des snapshots dans l'espace réduit. (a) L'échelle de couleurs porte sur les variations de $q_2$. (b) L'échelle de couleurs porte sur les variations de $\kappa$.</figcaption></figure>

Ces derniers sont tracés selon deux échelles de couleurs : le graphique IV.8a prend ses couleurs des variations de $q_2$ et le graphique IV.8b des variations de $\kappa$. Cette représentation est très
utile, car elle permet d'observer l'influence des paramètres. En effet, on voit que la manifold a bien trois dimensions avec cependant une différence de valeurs assez importante entre $α_1$, $α_2$ et $α_3$ ($α_2 < α_1/10$ et $α_3 < α_1/100$). La plus grande variance du modèle est donc principalement décrite selon $α_1$. Si on regarde bien l'échelle de couleurs (liée à $q_2$) sur la Figure IV.8a, on remarque qu'elle varie selon $α_1$. Cela signifie que le paramètre $q_2$ est celui qui influe le plus le comportement. La Figure IV.8b (couleurs liées à $κ$) indique que $κ$ varie suivant l'axe $α_3$ dont la variance est très faible comparée aux deux autres axes. Il sera plus difficile d'identifier $κ$ car cela signifie qu'il a une faible influence sur le comportement en compression hydrostatique, sur la plage considérée allant jusqu'au comportement asymptotique. La procédure appliquée au matériau R33 livre les mêmes conclusions, les graphiques ne seront donc pas représentés.

<table><caption>TABLEAU IV.1 – Paramètres de la pression hydrostatique identifiés pour le R23 et le R33.</caption><thead><tr><td>Matériau</td><td>s<sub>2</sub></td><td>q<sub>2</sub></td><td>κ [MPa]</td></tr></thead><tbody><tr><td>R23</td><td>0,593</td><td>0,918</td><td>48,599</td></tr><tr><td>R33</td><td>0,217</td><td>0,982</td><td>236,368</td></tr></tbody></table>

La conclusion de toutes ces observations est que $q_2$ et $s_2$ pourront être identifiés, mais la valeur obtenue pour $κ$ ne sera peut-être pas très représentative. Les valeurs de paramètres identifiées sont listées dans le Tableau IV.1 et les courbes pression-variation de volume correspondantes sont représentées sur la Figure IV.9.

<figure><img src="image_53.png" /><figcaption>FIGURE IV.9 – Résultats d'identification sur la compression hydrostatique pour le R23 (a) et le R33 (b).</figcaption></figure>

On remarque que les paramètres identifiés pour le R23 représentent globalement très bien le
comportement expérimental, même si le module de compressibilité $\kappa$ semble légèrement faible. Quant aux résultats sur le R33 sur la Figure IV.9b, on remarque que les valeurs identifiées sont moins aptes à modéliser le comportement. En effet, l'écart avec les données expérimentales est assez élevé au pied de courbe ainsi qu'à l'asymptote verticale, marquant le retour à un comportement incompressible. En pied de courbe, le premier point d'inflexion, qui prend sa position grâce au paramètre $s_2$, apparaît très tôt et en effet, la valeur de $s_2$ paraît assez faible. La valeur de $\kappa$ identifiée semble aberrante, car beaucoup trop grande et contre-intuitive puisqu'elle est plus grande que celle obtenue pour le R23. L'inclinaison du plateau (0,75<J<0,95) est plutôt bien reproduite. Cela signifie que l'identification du paramètre $q_2$ est pertinente. On remarque d'ailleurs que la valeur est du même ordre de grandeur pour les deux matériaux. Les hypothèses émises à l'observation de la Figure IV.8 semblent donc se confirmer après l'identification des paramètres, puisqu'il est difficile d'identifier $\kappa$. D'ailleurs en observant les valeurs et les agrandissements des pieds de courbes sur la Figure IV.9, on voit dans le cas du R23, par exemple, que le module $\kappa$ est sous-estimé et que l'inflexion arrive donc à une plus grande variation de volume. Au contraire pour le R33, la valeur de $\kappa$ anormalement élevée entraîne une grande diminution de $s_2$ et donc un point d'inflexion arrivant à de faibles valeurs de $J$.

Le fait que $\kappa$ soit difficile à identifier peut venir du fait que les courbes expérimentales présentent peu de points aux faibles pressions et donc un manque d'information dans cette zone.

**Identification de $q_2$ et $s_2$** - Afin d'améliorer la qualité de l'identification, l'idée a été ensuite de déterminer le module $\kappa$ en pied de courbe manuellement. Les valeurs sélectionnées sont la moyenne pour chaque essai de reproductibilité. On trouve ainsi des valeurs plus logiques, dans le sens où le module décrit en fonction du pourcentage de microsphères (72 MPa pour le R23, 54 MPa pour le R33). Ayant les valeurs du module pour nos deux matériaux, il est possible de mener l'identification de seulement 2 paramètres, $s_2$ et $q_2$. De plus, pour avoir un résultat optimal répondant aux applications visées, l'identification a été réalisée pour des déformations volumiques réduites à 15%.

La Figure IV.10a montre le plan d'expériences à deux paramètres généré pour le R23 et les snapshots correspondants à ces valeurs. L'échelle colorimétrique représente l'évolution de $q_2$, pour ces derniers. On remarque sur la Figure IV.10b que le plateau de la courbe est caractéristique des valeurs de $q_2$ élevées.

La Figure IV.10c est la représentation des snapshots dans l'espace des $\alpha$. La manifold obtenue après krigeage est bi-dimensionnelle et indique que les deux paramètres sont identifiables. De plus le snapshot expérimental (point noir) est dans la zone de la manifold. Les valeurs de $\alpha_1$ et $\alpha_2$ montrent comme dans le cas précédent que le paramètre $q_2$ a le plus d'influence sur le comportement. Sa valeur identifiée devrait donc être optimale.
<table><thead><tr><th>Matériau</th><th>s<sub>2</sub></th><th>q<sub>2</sub></th><th>κ [MPa]</th></tr></thead><tbody><tr><td>R23</td><td>0,472</td><td>0,933</td><td>72</td></tr><tr><td>R33</td><td>0,456</td><td>0,93</td><td>54</td></tr></tbody></table>

TABLEAU IV.2 – Paramètres s₂ et q₂ de la pression hydrostatique identifiés pour le R23 et le R33. La paramètre κ, module en pied de courbe, a été calculé pour chaque essai et la valeur moyenne a été sélectionnée.

Pour le cas du R33 (Figure IV.11), le même plan d'expériences a été généré. Le comportement de la manifold est donc similaire et les valeurs des coefficients α₁ et α₂ sont du même ordre de grandeur que pour le R23. Le snapshot expérimental du R33 est également très proche de la manifold. Les paramètres doivent être identifiables.

Les valeurs identifiées pour chaque matériau sont contenues dans le Tableau IV.2.

<figure><img src="image_54.png" /><figcaption>FIGURE IV.10 – Étapes d'identification des paramètres s₂ et q₂ de la pression hydrostatique sur le matériau R23.</figcaption></figure>
Les Figures IV.10d et IV.11d superposent les courbes expérimentales et identifiées respectivement pour le R23 et le R33. Les résultats sont représentés sur la pleine échelle de compressibilité des matériaux. La modélisation du comportement des matériaux dans la plage de déformation considérée est de très bonne qualité et pour des déformations supérieures, l'asymptote verticale est bien positionnée.

<figure><img src="image_55.png" /><figcaption>FIGURE IV.11 – Étapes d'identification des paramètres $s_2$ et $q_2$ de la pression hydrostatique sur le matériau R33.</figcaption></figure>
IV.B.3 Identification sur des chargements mixtes déviatorique/sphérique

REMARQUE IV.2

Dans la partie précédente, les snapshots étaient calculés analytiquement grâce à l'expression de la pression hydrostatique (IV.11) qui est simple à programmer et ne comporte aucune inconnue. Les parties suivantes impliquent des chargements plus complexes comme la traction uniaxiale dont la contrainte comporte, dans le cas de nos matériaux compressibles, une inconnue : la déformation transversale de l'éprouvette. Il devient donc moins aisé de calculer analytiquement la contrainte. C'est pourquoi à partir de maintenant, les snapshots seront calculés grâce à des simulations Abaqus effectuées sur un élément cubique C3D8 à partir du modèle implémenté dans une UMAT.

La théorie de MOERMAN et al. 2020 précise que les paramètres de la contrainte hydrostatique prennent des valeurs différentes selon le type de chargement induisant la variation de volume. Ainsi, les valeurs que nous venons d'identifier conviennent normalement dans le cas de la compression hydrostatique. Ces valeurs pourraient cependant ne pas être adaptées dans le cas de la traction uniaxiale par exemple, ou tout autre chargement où le volume augmente. On rappelle l'expression de la contrainte de Cauchy totale :

$$ \sigma = 2J^{-1} \left[ \begin{aligned} & [C_{10} + 3C_{30} (I_{\bar{B}} - 3)^2] \bar{B} - C_{01} \bar{B}^{-1} \\ & + \kappa \left[ (1-q)a \tan \left( \frac{J-1}{a} \right) + qb \tanh \left( \frac{J-1}{b} \right) \right] \end{aligned} \right] : \mathbb{P}_1 \qquad (IV.12) $$

$$ \text{avec } a = \frac{2}{\pi} \begin{cases} (J_1 - 1) & \text{si } J \ge 1 \\ (J_2 - 1) & \text{si } J < 1 \end{cases}, \quad b = \frac{1}{\kappa} \begin{cases} s_1 & \text{si } J \ge 1 \\ s_2 & \text{si } J < 1 \end{cases} \quad \text{et} \quad q = \begin{cases} q_1 & \text{si } J \ge 1 \\ q_2 & \text{si } J < 1 \end{cases} \qquad (IV.13) $$

Le paramètre $\kappa$, s'apparentant au module de compressibilité, doit logiquement rester identique. $s_1$ et $q_1$ pourraient avoir une valeur différente de $s_2$ et $q_2$. Cependant, le rôle du paramètre $J_1$ n'est pas clair. Il doit marquer la position d'une autre asymptote verticale de la contrainte hydrostatique dans le cas d'une augmentation de volume, ce qui semble difficile à imaginer sous sollicitation sphérique. Afin de vérifier l'influence de ces paramètres, des simulations de traction uniaxiale (pilotage en déplacement) ont été réalisées selon deux cas :

* $J_1$ varie entre 1 et 10, les autres paramètres prennent une valeur fixe,

* $s_1$ et $q_1$ varient respectivement dans les intervalles (0;2) et (0;1), les autres paramètres prennent une valeur fixe.

La Figure IV.12a montre les deux parties de la contrainte de Cauchy dans le premier cas. On remarque que $J_1$ n'a aucune influence sur le comportement de la partie sphérique de $\sigma$. La valeur de $J_1$ sera donc fixée arbitrairement à 4 afin de s'assurer de la stabilité du modèle (cet aspect sera détaillé dans la suite du chapitre). La Figure IV.12b montre une légère influence de la partie sphérique en traction au-delà de 80% de déformation lorsque les paramètres $s_1$ et $q_1$ varient. Dans les gammes de déformations considérées, le comportement reste quasiment inchangé. Les
<figure><img src="image_56.png" /><figcaption>FIGURE IV.12 – Simulations de traction uniaxiale pour des valeurs fixes des $C_{ij}$ ($C_{10}=-0,1$ MPa, $C_{01}=0,6$ MPa et $C_{30}=0,01$ MPa). (a) Les paramètres de la partie volumique de l'énergie libre sont fixes, seul $J_1$ varie de 1 à 10. (b) Les paramètres de la partie volumique de l'énergie libre varient respectivement dans les intervalles (0;2) et (0;1) pour $s_1$ et $q_1$ et entre 40 et 80 MPa pour $\kappa$, $J_1=4$.</figcaption></figure>

paramètres n'ont donc pas d'influence marquée en traction. Nous allons donc faire l'hypothèse que les paramètres $s_1$ et $q_1$ sont respectivement égaux à $s_2$ et $q_2$ précédemment identifiés en compression hydrostatique. Ainsi seuls les paramètres $C_{10}$, $C_{01}$ et $C_{30}$ seront nécessaires à la détermination de la réponse en traction.

Dans le cas de la compression confinée, la forme typique non-linéaire de la courbe contraint variation de volume est clairement pilotée par la valeur de la contrainte sphérique (voir la courbe bleue sur la Figure IV.13a). Ceci ne sous-entend pas pour autant que la réponse est purement sphérique, puisque les courbes rouges sur cette même figure montrent que la partie déviatorique de la réponse n'est pas nulle et qu'elle évolue selon les valeurs prises par $C_{10}$, $C_{01}$ et $C_{30}$. Sur la courbe de droite IV.13b, la contrainte de Piola-Kirchhoff I totale est représentée. On voit que le pied de courbe reste le même. Il est donc entièrement déterminé par les paramètres de la contrainte hydrostatique. Ensuite, le comportement est influencé par les valeurs prises par $C_{10}$, $C_{01}$ et $C_{30}$.

Ces observations laissent penser qu'identifier seulement les valeurs des paramètres déviato-riques tout en conservant ceux de la partie sphérique est possible. Nous tenterons donc d'iden-tifier les paramètres $C_{10}$, $C_{01}$ et $C_{30}$ simultanément sur des essais de traction simple et de compression confinée.

**Identification simultanée traction/compression confinée** Afin de réaliser l'identification de $C_{10}$, $C_{01}$ et $C_{30}$, on réalise des simulations de compression confinée et de traction uniaxiale selon le plan d'expériences de la Figure IV.14. Il est obtenu par la méthode de LHS dont on a conservé que les points dont la somme de $C_{10}$ et $C_{01}$ est positive pour assurer la stabilité du modèle. En effet, des tests préliminaires ont montré que pour reproduire le comportement de nos matériaux $C_{10}$ était négatif, avec $C_{10} + C_{01} > 0$.
<figure><img src="image_57.png" /></figure>

FIGURE IV.13 – Simulations de compression confinée sur le R23. (a) .(b) . Les paramètres de la partie isochore de l'énergie varient : -0,5 < C₁₀ < 1 MPa, 0 < C₀₁ < 1 MPa et 0 < C₃₀ < 0,1 MPa. Les paramètres de la partie volumique de l'énergie sont fixes s₁=0,472, q₁=0,933 et κ=72 MPa.

<figure><img src="image_58.png" /></figure>

FIGURE IV.14 – Plan d'expériences tri-dimensionnel sur C₁₀, C₀₁ et C₃₀.

On obtient les snapshots de l'évolution de la contrainte de Piola-Kirchhoff I représentés sur la Figure IV.15 (en fonction de la variation de volume et de l'élongation longitudinale respectivement pour la compression confinée et la traction uniaxiale). Les courbes en croix noires correspondent aux snapshots expérimentaux. Les snapshots en bleu (sur la Figure IV.15b) sont négatifs en raison de faibles valeurs des paramètres Cᵢⱼ mais ne gêneront pas leur identification. L'étape suivante consiste à concatener les snapshots de compression et de traction afin de réaliser la POD sur l'ensemble. Le snapshot expérimental est concatené de manière analogue.

La Figure IV.16 permet comme précédemment d'estimer l'influence des paramètres sur le modèle pour le matériau R23. La manifold obtenue après l'application de la POD est représentée avec une échelle de couleur dépendant de C₁₀, C₀₁ ou C₃₀. Le paramètre C₁₀ varie principalement
<figure><img src="image_59.png" /><figcaption>FIGURE IV.15 – Snapshots de la contrainte de Piola-Kirchhoff I pour la compression confinée (a) et la traction uniaxiale (b). Les courbes en croix noires correspondent aux snapshots expérimentaux.</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>FIGURE IV.16 – Représentation graphique du krigeage des snapshots dans la base de dimension réduite. L’échelle de couleurs porte sur les variations de : (a) $C_{10}$, (b) $C_{01}$ et (c) $C_{30}$.</figcaption></figure>

ment selon l’axe $\alpha_1$, de plus forte variance, il a donc plus d’importance/d’influence que les deux autres paramètres sur le comportement en traction et compression confinée. Les Figures IV.16b et IV.16c montrent que l’échelle de couleur varie dans le plan $(\alpha_2, \alpha_3)$. $C_{01}$ et $C_{30}$ doivent avoir une influence similaire.

La manifold est tri-dimensionnelle et le snapshot expérimental (étoile noire) se trouve à l’intérieur. Les trois paramètres doivent pouvoir être identifiés.

On retrouve la comparaison entre l’expérience et le modèle sur la Figure IV.17 pour le R23. La sollicitation en traction est plutôt bien reproduite par le modèle. En compression confinée, la contrainte en pied de courbe est surestimée. Néanmoins les résultats expérimentaux sont à prendre avec précaution, car, comme évoqué dans la partie II.B.4.a, il peut y avoir un certain temps de mise en place de l’éprouvette lorsque l’essai commence. Les premiers points de mesure peuvent être de la compression uniaxiale, le temps que l’éprouvette touche parfaitement les
<table><thead><tr><th>Matériau</th><th>C<sub>10</sub> [MPa]</th><th>C<sub>01</sub> [MPa]</th><th>C<sub>30</sub> [MPa]</th></tr></thead><tbody><tr><td>R23</td><td>-0,296</td><td>0,779</td><td>0,030</td></tr><tr><td>R33</td><td>-0,627</td><td>1,188</td><td>0,049</td></tr></tbody></table>

TABLEAU IV.3 – Paramètres C₁₀, C₀₁ et C₃₀ de la partie déviatorique de la contrainte de Cauchy identifiés en traction et compression confinée pour le R23 et le R33.

bords de la cavité cylindrique et que le confinement augmente.

La démarche est similaire dans le cas du R33. On retrouve les résultats dans le Tableau IV.3 et sur la Figure IV.18. La courbe identifiée en traction pour le R33 est moins satisfaisante que celle obtenue pour le R23. Cela peut être dû à l'incapacité du modèle à retranscrire la forte rigidité en pied de courbe suivie de la zone en plateau où la contrainte augmente lentement.

<figure><img src="image_61.png" /><figcaption>FIGURE IV.17 – Résultats d'identification du R23 en compression confinée (a) et en traction uniaxiale (b).</figcaption></figure>

IV.B.4 Identification de la partie liée à l'effet Mullins

Dans cette partie, nous allons identifier les paramètres $\bar{r}$, $\bar{m}$, $\bar{\beta}$, $\hat{r}$, $\hat{m}$ et $\hat{\beta}$ servant respectivement au calcul de la variable $\bar{\eta}$ et $\hat{\eta}$ de l'effet Mullins. Ce dernier est appliqué à la fois sur la partie déviatorique et sphérique de la contrainte de Cauchy. Cela fait donc 6 paramètres à identifier. On rappelle l'expression de $\bar{\eta}$ et $\hat{\eta}$ :

$$ \bar{\eta} = 1 - \frac{1}{\hat{r}}\operatorname{erf}\left(\frac{\rho_0 \bar{\psi}^{\text{max}} - \rho_0 \hat{\psi}}{\bar{m} + \bar{\beta} \bar{\rho}_0 \bar{\psi}^{\text{max}}}\right) \qquad \hat{\eta} = 1 - \frac{1}{\hat{r}}\operatorname{erf}\left(\frac{\rho_0 \hat{\psi}^{\text{max}} - \rho_0 \hat{\psi}}{\hat{m} + \hat{\beta} \hat{\rho}_0 \hat{\psi}^{\text{max}}}\right) \quad (\text{IV.14}) $$

L'idéal aurait été d'utiliser une sollicitation cyclée purement déviatorique ou purement sphérique pour identifier les paramètres séparément. Seulement, les essais de cisaillement dont nous disposons ne reflètent peut-être pas une sollicitation purement déviatorique à cause du montage
<figure><img src="image_62.png" /><figcaption>FIGURE IV.18 – Résultats d'identification du R33 en compression confinée (a) et en traction uniaxiale (b).</figcaption></figure>

(II.B.4.a). Quant aux essais de compression sphérique au mercure, ils ne permettent pas de gérer la vitesse de sollicitation, ni de réaliser un cyclage. De plus, il n'est pas possible d'avoir un temps associé aux mesures de variations de volume. Ils ne sont donc pas adaptés à l'identification de phénomènes dissipatifs. Les seuls essais à disposition sont alors des essais cyclés à amplitude croissante de déformation, en traction uniaxiale et compression confinée.

## IV.B.4.a Partie déviatorique

La démarche adoptée dans la partie IV.B.3 a montré que la contrainte sphérique n'avait que peu d'influence sur le comportement global en traction uniaxiale. C'est pourquoi, dans ce cas, les paramètres *r̂*, *m̂* et *β̂* sont respectivement choisis égaux à *r̄*, *m̄* et *β̄*. Le plan d'expérience créé pour cette partie est représenté sur la Figure IV.19a. Les bornes de celui-ci ont été choisies en fonction de valeurs trouvées dans la littérature (GRACIA et al. 2009; TOUPS et al. 2020). Ils doivent tous trois être positifs, *m* et *β* peuvent être nuls mais non simultanément. En général *m* est assez faible, il variera donc entre 0 et 0,1 MPa. Le cyclage utilisé pour les simulations est basé sur l'amplitude de déformation observée expérimentalement. Il y a 5 cycles de 10 à 50% de déformation. Les snapshots de contrainte de Piola-Kirchhoff I correspondants pour le R23 sont tracés sur la Figure IV.20a en fonction du temps d'essai. On remarque que certains jeux de paramètres conduisent à des valeurs de contrainte négatives. Il s'agit vraisemblablement d'une combinaison de valeurs trop faibles de *r* et *β*. Bien que non physiques, ces résultats ne devraient pas gêner l'identification des trois paramètres. La manifold obtenue après la POD est représentée sur la Figure IV.19b. Elle n'a plus une forme cubique comme le plan d'expériences mais plutôt conique. Les points situés en pointe de ce cône correspondent à des valeurs entre 1,9 et 2 (vert) pour les paramètres *r̄* et *β̄*. Cela veut dire que, pour ces valeurs, le comportement varie peu. En observant bien la couleur des snapshots sur la Figure IV.20a, on remarque que les courbes vertes correspondent aux contraintes les plus élevées et que leur forme
<figure><img src="image_63.png" /><figcaption>FIGURE IV.19 – (a) Plan d’expériences tri-dimensionnel sur $\bar{r}$, $\bar{m}$ et $\bar{\beta}$. (b) Snapshots de compression confinée sur le R23 représentés dans l’espace des $\alpha$ après la POD.</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>FIGURE IV.20 – Identification de $\bar{r}$, $\bar{m}$ et $\bar{\beta}$ pour le calcul de $\bar{\eta}$ en traction uniaxiale pour le R23. (a) Snapshots obtenus avec le plan d’expériences de la Figure. (b) Courbe identifiée pour le R23.</figcaption></figure>

varie peu, contrairement par exemple, aux courbes bleues/violettes qui ont un comportement différent. Malgré sa forme, la manifold reste tri-dimensionnelle et le snapshots expérimental du R23 n’en est pas très éloigné. Nous pouvons donc poursuivre avec le krigeage et l’identification qui conduisent à de bons résultats. La Figure IV.20b superpose les courbes expérimentale et identifiée avec une très bonne adéquation. Les paramètres correspondants peuvent être trouvés dans le Tableau IV.4.
<table><thead><tr><th>Matériau</th><th>r&#773; [-]</th><th>m&#773; [MPa]</th><th>&beta;&#773; [-]</th></tr></thead><tbody><tr><td>R23</td><td>1,580</td><td>0,089</td><td>0,467</td></tr><tr><td>R33</td><td>1,574</td><td>0,064</td><td>0,320</td></tr></tbody></table>

TABLEAU IV.4 – Paramètres $\bar{r}$, $\bar{m}$ et $\bar{\beta}$ de la partie déviatorique de la contrainte de Cauchy identifiés en traction pour le R23 et le R33.

Pour le R33, le plan d'expériences a été amélioré en tenant compte des observations sur la manifold du R23. Le paramètre $\bar{r}$ varie maintenant entre 1 et 2, $\bar{\beta}$ entre 0,4 et 1 (voir Figure IV.21a). Les snapshots contrainte-temps obtenus sont tracés sur la Figure IV.22a. Les formes instables présentes pour le R23 sont cette fois évitées. Quant à la manifold, représentée sur la Figure IV.21b, elle est bien tri-dimensionnelle et sa forme est améliorée grâce au nouveau plan d'expériences. L'échelle des $\alpha$ est plus grande, les paramètres sont donc à l'origine d'une plus grande variance que pour le R23. Le point expérimental est très proche des autres points, donc les paramètres correspondants devraient pouvoir être identifiés assez précisément. La courbe obtenue est en effet proche des données expérimentales, comme le montre la Figure IV.22b et les paramètres identifiés se trouvent dans le Tableau IV.4.

<figure><img src="image_65.png" /></figure>

FIGURE IV.21 – (a) Plan d'expériences tri-dimensionnel sur $r$, $m$ et $\beta$. (b) Snapshots de traction uniaxiale sur le R33 représentés dans l'espace des $\alpha$ après la POD.

### IV.B.4.b Partie sphérique

Une fois que l'identification des paramètres permettant le calcul de $\bar{\eta}$ est réalisée, ceux de la partie sphérique sont identifiés sur les essais de compression confinée. Les paramètres $\bar{r}$, $\bar{m}$ et $\bar{\beta}$ prennent les valeurs du Tableau IV.4 et $\hat{r}$, $\hat{m}$ et $\hat{\beta}$ varient selon le plan d'expériences de la Figure IV.19a. Dans cette partie, nous disposons d'essais de compression confinée à amplitude de déformation croissante jusqu'à 15%. Il y a trois amplitudes pour lesquelles on compte dix
<figure><img src="image_66.png" /><figcaption>FIGURE IV.22 – Identification de r, m et β pour le calcul de $\bar{\eta}$ en traction uniaxiale pour le R33. (a) Snapshots obtenus avec le plan d'expériences de la Figure. (b) Courbe identifiée pour le R33.</figcaption></figure>

cycles. La contrainte résultante est représentée en fonction du temps pour le plan d'expériences et la courbe expérimentale sur les Figures IV.24a et IV.25a pour respectivement le R23 et le R33.

<figure><img src="image_67.png" /><figcaption>FIGURE IV.23 – (a) Plan d'expériences tri-dimensionnel sur $\hat{r}$, $\hat{m}$ et $\hat{\beta}$. (b) Snapshots de compression confinée du R23 représentés dans l'espace des $\alpha$ après la POD.</figcaption></figure>
<table><thead><tr><th>Matériau</th><th>r̂ [-]</th><th>m̂ [MPa]</th><th>β̂ [-]</th></tr></thead><tbody><tr><td>R23</td><td>1,64</td><td>0,00001</td><td>1,05</td></tr><tr><td>R33</td><td>1,86</td><td>0,0015</td><td>0,643</td></tr></tbody></table>

TABLEAU IV.5 – Paramètres r̂, m̂ et β̂ de la partie sphérique de la contrainte de Cauchy identifiés en compression confinée pour le R23 et le R33.

<figure><img src="image_68.png" /></figure>

FIGURE IV.24 – Identification de r̂, m̂ et β̂ pour le calcul de η̂ en compression confinée pour le R23. (a) Snapshots. (b) Courbe identifiée pour le R23.

Les mêmes remarques que précédemment peuvent être faites quant à la forme de la manifold obtenue dans l'espace des α (voir Figure IV.23b). Elle sera représentée uniquement pour le R23 car le R33 est très similaire. On retient sa tri-dimensionnalité qui indique l'indentifiabilité des trois paramètres.

Après krigeage et optimisation, ces derniers ont pu être identifiés et sont explicités dans le Tableau IV.5. Les courbes comparatives sont tracées pour le R23 et le R33 sur les Figures IV.24b et IV.25b, et présentent une bonne adéquation.

## IV.C Stabilité du modèle

Le critère de Drucker, dans le cas d'un matériau hyperélastique compressible, assure la stabilité du modèle si la variation de la contrainte de Kirchhoff suivant une variation infinitésimale de déformation logarithmique satisfait l'inégalité (ABAQUS 2020a) :

$$d\tau : d\varepsilon > 0 \qquad (IV.15)$$

Or l'incrément de contrainte peut également s'écrire :

$$d\tau = D : d\varepsilon \qquad (IV.16)$$
<figure><img src="image_69.png" /><figcaption>FIGURE IV.25 – Identification de $\hat{r}$, $\hat{m}$ et $\hat{\beta}$ pour le calcul de $\hat{\eta}$ en compression confinée pour le R33. (a) Snapshots. (b) Courbe identifiée pour le R33.</figcaption></figure>

ou sous forme matricielle et en se plaçant dans l'espace des directions principales :

$$ \left\{ \begin{matrix} d\tau_1 \\ d\tau_2 \\ d\tau_3 \end{matrix} \right\} = \left[ \begin{matrix} D_{11} & D_{12} & D_{13} \\ D_{21} & D_{22} & D_{23} \\ D_{31} & D_{32} & D_{33} \end{matrix} \right] \left\{ \begin{matrix} d\epsilon_1 \\ d\epsilon_2 \\ d\epsilon_3 \end{matrix} \right\} \qquad (IV.17) $$

avec $D_{ij}(\lambda_1, \lambda_2, \lambda_3) = \frac{\partial \tau_i}{\partial \varepsilon_j}$, $\varepsilon$ représentant la déformation logarithmique telle que $\varepsilon_i = \ln(\lambda_i)$.

Les conditions nécessaires et suffisantes à la stabilité du modèle sont :

$$ \left\{ \begin{matrix} D_{11} + D_{22} + D_{33} > 0 \\ D_{11}D_{22} + D_{22}D_{33} + D_{33}D_{11} - D_{23}D_{32} - D_{13}D_{31} - D_{12}D_{21} > 0 \\ \det(\mathbf{D}) > 0 \end{matrix} \right. \qquad (IV.18) $$

L'obtention de la matrice **D** se fait grâce au calcul des variations des contraintes principales $\tau_i$ qui s'expriment :

$$ \tau_i = 2 \left[ J^{-\frac{2}{3}} \frac{\partial \bar{\psi}}{\partial I_{\bar{\mathbf{B}}}} \left( \lambda_i - \frac{1}{3} \mathrm{tr} \bar{\mathbf{B}} \right) - J^{\frac{2}{3}} \frac{\partial \bar{\psi}}{\partial \Pi_{\bar{\mathbf{B}}}} \left( \lambda_i^{-1} - \frac{1}{3} \mathrm{tr} \bar{\mathbf{B}}^{-1} \right) \right] + J \frac{\partial \bar{\psi}}{\partial J} \qquad (IV.19) $$

Leurs variations $d\tau_i$ sont calculées avec la relation (IV.20) qui fait apparaître les composantes de **D** :

$$ d\tau_i = \frac{\partial \tau_i}{\partial \lambda_j} d\lambda_j = \underbrace{\frac{\partial \tau_i}{\partial \lambda_j} \lambda_j}_{D_{ij}} \underbrace{\frac{d\lambda_j}{d\varepsilon_j}}_{d\varepsilon_j} \qquad (IV.20) $$

Une fois toutes les dérivées calculées, les conditions (IV.18), assurant le stabilité du modèle, peuvent être vérifiées.
Le développement des équations est cependant assez fastidieux et ne pourra pas être réalisé
dans le cadre de ce travail faute de temps.

Néanmoins, la formulation adoptée permet d'exprimer le critère en fonction des parties déviatoriques et sphériques :

$$
d\tau : d\epsilon > 0 \Leftarrow (d\bar{\tau} + d\hat{\tau}) : (d\bar{\epsilon} + d\hat{\epsilon}) > 0 \Leftarrow \begin{cases} d\bar{\tau} : d\bar{\epsilon} > 0 \\ d\hat{\tau} : d\hat{\epsilon} > 0 \end{cases} \qquad (IV.21)
$$

La stabilité du modèle peut donc être vérifiée avec la positivité des deux conditions de l'équation
(IV.21) et c'est une condition suffisante. La démarche adoptée sera donc en deux étapes :

■ La stabilité de la partie isochore du modèle est assurée par deux conditions :

— l'énergie libre $\bar{\psi}$ est nulle à l'origine, positive et convexe,

— le critère de Drucker est vérifié dans divers cas particuliers de sollicitations grâce au module d'évaluation d'Abaqus.

■ La stabilité de la partie volumique est assurée en vérifiant que :

— l'énergie libre $\hat{\psi}$ est nulle à l'origine, positive et convexe,

— le critère de Drucker est vérifié.

REMARQUE IV.3

La démarche permet de vérifier la stabilité du modèle pour des chargements mono-
tones. En effet, cela correspond à des valeurs de $\bar{\eta}$ et $\hat{\eta}$ égales à 1. Dès lors qu'un
chargement cyclé est considéré, $\bar{\eta}$ et $\hat{\eta}$ vont varier en fonction des valeurs prises par les
paramètres $\bar{r}$, $\bar{m}$, $\bar{\beta}$, $\hat{r}$, $\hat{m}$ et $\hat{\beta}$. La vérification de la stabilité dans ce cas serait possible
en s'assurant de la positivité de $\mathbf{D}$ comme explicité en (IV.20).

Les résultats fournis par Abaqus permettent de montrer que la partie isochore du modèle pour le R23 est stable pour toutes les sollicitations testées (listées dans le Tableau IV.6). Quant au R33, il devient instable en compression uniaxiale pour des déformations nominales inférieures à -28% et en traction biaxiale après 18% de déformation nominale. Il est peu probable, dans l'application du matériau, que celui-ci soit soumis à de la traction biaxiale, donc ce résultat ne devrait pas être un problème majeur. De même, au vu des faibles déformations appliquées au matériau, la compression uniaxiale devrait rester stable et supérieure à -28%.

Nous avons également vérifié, pour la partie isochore, que l'énergie libre $\bar{\psi}$ utilisée était positive et convexe. Sa valeur a été calculée pour toutes les valeurs de $\lambda_1$ et $\lambda_2$ comprises entre 0,6 et 2, $\lambda_3$ est choisi comme étant égal à $\lambda_1^{-1}\lambda_2^{-1}$ afin de satisfaire l'incompressibilité. Les surfaces obtenues pour le R23 et le R33, représentées sur la Figure IV.26, sont bien convexes et positives. Une autre condition nécessaire est vérifiée : l'hypothèse émise lors de l'identification des paramètres $C_{ij}$ (que $C_{10}$ pouvait être négatif, tant que $C_{10} + C_{01} > 0$ et $C_{30} > 0$). Cela permet de satisfaire la positivité de $\bar{\psi}$.
<table><caption>TABLEAU IV.6 – Stabilité de la partie isochore du modèle pour le R23 et le R33. Calcul des critères de Drucker à l'aide du module d'évaluation d'Abaqus.</caption><thead><tr><th>Sollicitation</th><th>R23</th><th>R33</th></tr></thead><tbody><tr><td>Traction uniaxiale</td><td>stable</td><td>stable</td></tr><tr><td>Compression uniaxiale</td><td>stable</td><td>stable jusqu'à -28% de déformation</td></tr><tr><td>Traction biaxiale</td><td>stable</td><td>stable jusqu'à 18% de déformation</td></tr><tr><td>Compression biaxiale</td><td>stable</td><td>stable</td></tr><tr><td>Traction plane</td><td>stable</td><td>stable</td></tr><tr><td>Compression plane</td><td>stable</td><td>stable</td></tr></tbody></table>

<figure><img src="image_70.png" /><figcaption>FIGURE IV.26 – Vérification de la positivité et convexité de l'énergie libre isochore pour les paramètres identifiés pour le matériau : (a) R23, (b) R33.</figcaption></figure>

L'énergie libre volumique utilisée dépend de fonctions trigonométriques et de leur logarithme népérien. La valeur $\hat{\psi}$ est alors calculable tant que les tangentes sont positives. Cela est assuré par les valeurs des paramètres $J_1$ et $J_2$. $\hat{\psi}$ retournera donc une valeur tant que $J_2 < J < J_1$. Pour le vérifier, nous avons fait varier $J$ dans cet intervalle et tracé l'énergie libre correspondante pour chaque matériau. Les résultats sur la Figure IV.27 montrent que $\hat{\psi}$ est positive, nulle à l'origine et convexe pour le R23 et le R33.

Le critère de Drucker appliqué à la contrainte hydrostatique s'écrit :

$$d\hat{\tau} : d\hat{\varepsilon} > 0 \qquad (IV.22)$$

Les valeurs de contraintes étant identiques sur la diagonale, et comme $\varepsilon = \bar{\varepsilon} + \hat{\varepsilon} = \ln \bar{V} + \frac{1}{3}\ln J$ cela revient à écrire :

$$d\hat{\tau} = \left( J \frac{\partial \hat{\psi}}{\partial J} + J^2 \frac{\partial^2 \hat{\psi}}{\partial J^2} \right) d\hat{\varepsilon} \qquad (IV.23)$$
<figure><img src="image_71.png" /></figure>

<figure><img src="image_72.png" /><figcaption>FIGURE IV.27 – Vérification de la positivité de l'énergie libre volumique pour les paramètres identifiés pour chaque matériau.</figcaption></figure>

Il convient alors de vérifier la positivité de :

$$ \frac{\partial \hat{\psi}}{\partial J} + J \frac{\partial^2 \hat{\psi}}{\partial J^2} \qquad (IV.24) $$

La Figure IV.28 montre que la valeur du critère reste positive quel que soit le matériau. On observe également que la valeur pour $J = 1$ n'est pas nulle et correspond à la dérivée seconde de $\hat{\psi}$, qui vaut $\kappa$ à déformation nulle. La partie volumique de notre modèle est stable entre $J_2$ et $J_1$, $J_2$ étant directement lié au pourcentage de microsphères dans les matériaux, et $J_1$ fixé arbitrairement à 4 comme expliqué dans la partie IV.B.3.

<figure><img src="image_73.png" /></figure>

<figure><img src="image_74.png" /><figcaption>FIGURE IV.28 – Vérification du critère de Drucker pour la partie sphérique de la contrainte, dans le cas de chaque matériau.</figcaption></figure>
# Conclusions du chapitre

Ce chapitre a permis l'identification des paramètres d'un modèle hyperélastique compres-
sible. Un processus d'identification robuste a été proposé et utilisé. Il présente de nombreux
avantages :

* La méthode d'identification par POD couplée au krigeage permet de mieux comprendre l'influence des paramètres sur les sollicitations considérées et cela est intéressant, car :

— si un ou plusieurs paramètres n'ont pas d'influence, la dimensionnalité dans la base réduite après POD sera inférieure au nombre de paramètres à identifier. L'unicité de l'identification n'est alors pas assurée (mais nous en sommes néanmoins conscients),

— au contraire, si la dimensionnalité et le nombre de paramètres sont égaux, ils peuvent être identifiés et l'unicité de la solution est assurée.

* Les données réduites par POD sont constituées de snapshots, ou vecteurs, contenant les valeurs de contrainte. Ces derniers sont obtenus par simulation de toute sollicitation qui peut s'avérer utile pour identifier les paramètres. Les vecteurs peuvent être concaténés pour permettre l'identification simultanée sur plusieurs cas de chargement. C'est même un avantage qui permet d'enrichir le modèle et augmente les chances d'identifiabilité des paramètres recherchés.

* Une fois la base de snapshots constituée et le krigeage réalisé, il est possible d'identifier les paramètres grâce aux snapshots expérimentaux de notre choix. Si un autre matériau ayant le même comportement était considéré, il suffirait de faire les essais expérimentaux idoines et de les projeter dans la base de POD. L'entraînement du modèle est certes couteux, mais n'est réalisé qu'une seule fois et réutilisable.

<figure><img src="image_75.png" /><figcaption>FIGURE IV.29 – Courbes identifiées contrainte-elongation pour le R23. (a) Traction uniaxiale. (b) Compression confinée.</figcaption></figure>

Finalement, après les développements de ce chapitre, le résultat final est la représentation
de la contrainte identifiée en fonction de l'élongation. Pour le R23 par exemple, on obtient les
graphiques de la Figure IV.29. On voit que le niveau de contrainte de l'enveloppe de charge est
bien respecté grâce à l'identification des $C_{ij}$, s, q et $κ$. Ensuite le positionnement des cycles est plutôt bien reproduit dans les deux cas par les valeurs de $\bar{r}$, $\bar{m}$, $\bar{\beta}$, $\hat{r}$, $\hat{m}$ et $\hat{\beta}$ identifiées.
L'hypothèse d'identification découplée selon la contrainte de Cauchy paraît au vu de ces résultats, satisfaisante. Il y a en tout 12 paramètres dans ce modèle hyperélastique compressible qu'il aurait été plus difficile d'identifier simultanément. Notre méthode a permis leur identification 3 à 3 en gardant une bonne cohérence globale. De plus, les paramètres identifiés assurent la stabilité pour le R23 et un domaine acceptable du point de vue de l'utilisation pour le R33.

Comme le montrent les courbes présentées Figure IV.29, il reste cependant, afin de reproduire
au mieux le comportement de ces matériaux et la déformation résiduelle qu’ils présentent, à
identifier des paramètres de la partie (visco-)plastique du modèle. Ce sera l’objet du prochain
chapitre.

Point sur les phénomènes pris en compte après le Chapitre 4

■ Phénomènes de 1<sup>er</sup> ordre :

✔ grandes déformations et non-linéarités,

✔ grande compressibilité,

✔ effet Mullins,

X déformation résiduelle,

▷ influence de la fraction volumique de microsphères.

■ Phénomènes de 2<sup>nd</sup> ordre :

X effet de la température,

X effet visqueux (Fluage, relaxation).
# Chapitre V

## Implémentation d'un modèle élasto-(visco-)plastique. Capacités et perspectives.

Ce chapitre présente l'implémentation du modèle complet élasto-(visco-)plastique avec effet Mullins dans le code de calculs aux éléments finis Abaqus. La première version de ce modèle est élasto-plastique car on choisit une faible viscosité afin de faciliter l'implémentation du caractère élasto-plastique. Les capacités de ce modèle sont simulées pour les cas de chargements cyclés qui nous intéressent à savoir la traction uniaxiale, les compressions confinée et hydrostatique et le cisaillement simple. Ces exemples mettent en évidence l'importance de la plasticité et la prise en compte de l'effet Mullins dans la modélisation du comportement de nos matériaux. C'est pourquoi la suite de ce chapitre est axée sur l'identification des 16 paramètres du modèle. Comme dans le chapitre précédent, nous voulons procéder par étapes, et identifier dans un premier temps les paramètres de la partie isochore à l'aide de simulations d'essais de traction cyclée à amplitude de déformation croissante. Cela réduit le nombre de paramètres à 8, qui pourront être partiellement identifiés. Cela permet de mettre en évidence un manque d'information du modèle de krigeage. Ce dernier gagnerait à contenir plus de types de sollicitations notamment purement déviationiques et/ou sphériques. Finalement, l'importance d'ajouter une prise en compte de plusieurs temps caractéristiques de viscosités est discutée, notamment au vu de l'application du matériau et de son comportement sous des sollicitations de type relaxation ou fluage. Nous proposerons une extension au modèle en ce sens.
PLAN DU CHAPITRE V

<table><tr><td>V.A</td><td>Implémentation d'un modèle élasto-(visco-)plastique</td><td>125</td></tr><tr><td>V.A.1</td><td>Implémentation</td><td>126</td></tr><tr><td>V.A.2</td><td>Test du modèle</td><td>127</td></tr><tr><td>V.B</td><td>Identification des paramètres</td><td>131</td></tr><tr><td>V.C</td><td>Possibilités du modèle élasto-(visco-)plastique</td><td>134</td></tr><tr><td>V.C.1</td><td>Viscosité à temps caractéristique court</td><td>134</td></tr><tr><td>V.C.1.a</td><td>Relaxation</td><td>134</td></tr><tr><td>V.C.1.b</td><td>Influence de la vitesse de sollicitation</td><td>135</td></tr><tr><td>V.C.2</td><td>Viscosité à temps caractéristique long</td><td>136</td></tr><tr><td colspan="2">Conclusions du chapitre</td><td>142</td></tr></table>
V.A Implémentation d'un modèle élasto-(visco-)plastique

Le chapitre précédent a permis d'identifier les paramètres pour un comportement simplifié de type hyperélastique compressible incluant l'effet Mullins. Du point de vue de l'application, cela permet dans un premier temps de faire des simulations en connaissant le comportement au premier chargement (et à la première décharge si on néglige la déformation permanente) des joints compressibles. Dans ce chapitre, l'objectif est, en premier lieu, de rendre compte également de la déformation résiduelle et de l'adoucissement de la contrainte sur une répétition de cycles de même amplitude.

Le modèle élasto-(visco-)plastique est décrit par les équations suivantes (cf. Chapitre 3 partie III.C.1.b) :

$$
\left\{
\begin{aligned}
\bar{\sigma} &= 2J^{-1}\bar{\eta} \left[ (C_{10} + 3C_{30}(I_{\bar{B}} - 3)^2) \bar{B} - C_{01}\bar{B}^{-1} \right] : \mathbb{P}_1 + 2J^{-1} \left[ C_p \bar{B}_e - 3 \right] : \mathbb{P}_1 \\
\hat{\sigma} &= \hat{\eta}\kappa \left[ (1-q)a \tan\left(\frac{J-1}{a}\right) + qb \tanh\left(\frac{J-1}{b}\right) \right] \mathbf{1} + J_a^{-1}\kappa_e \left[ a \tan\left(\frac{J_e-1}{a}\right) \right] \mathbf{1} \\
\dot{\bar{B}}_e &= -\frac{2}{3} \left[ \mathbf{1} : \mathbf{D} \right] \bar{B}_e + \frac{2}{3} \left[ \mathbf{1} : \mathbf{D}_a^\circ \right] \bar{B}_e + \mathbf{L} \cdot \bar{B}_e + \bar{B}_e \cdot \mathbf{L}^T - 2\bar{B}_e \cdot \left[ \frac{\langle\!\langle|\bar{\sigma}_{neq}^D|| - \bar{\chi}\rangle\!\rangle}{\bar{\nu}} \frac{\bar{\sigma}_{neq}^D}{|\bar{\sigma}_{neq}^D|} \right] \\
J_a^{-1} \dot{J}_a &= \frac{\langle\!\langle|\hat{\sigma}_{neq}|| - \hat{\chi}\rangle\!\rangle (\hat{\sigma}_{neq} : \mathbf{1})}{\hat{\nu} |\hat{\sigma}_{neq}||} \quad (V.1)
\end{aligned}
\right.
$$

L'effet Mullins est intégré avec le modèle d'OGDEN & ROXBURGH 1999, utilisé à la fois sur la
partie isochore et la partie volumique de l'énergie libre ψ<sub>eq</sub>.

La (visco-)plasticité est introduite également sur chacune des deux parties pour garantir un comportement dissipatif et une certaine déformation résiduelle, quelle que soit la sollicitation (donc y compris sphérique). Les 18 paramètres du modèle (V.1) devant être identifiés sont :
$C_{10}$, $C_{01}$, $C_{30}$, $C_p$, $\bar{\chi}$, $\bar{\nu}$, $\kappa$, $s$, $q$, $\kappa_e$, $\hat{\chi}$, $\hat{\nu}$, $\bar{r}$, $\bar{m}$, $\bar{\beta}$, $\hat{r}$, $\hat{m}$ et $\hat{\beta}$.

<figure><img src="image_76.png" /><figcaption>FIGURE V.1 – Simulations de traction uniaxiale cyclée. La viscosité $\bar{\nu}$ varie entre $10^{-4}$ et $10^{-1}$ MPa.s.</figcaption></figure>
Le terme « visco » est employé entre parenthèses, car notre but premier est de rendre compte du caractère plastique du comportement. La viscoplasticité est utilisée dans ce but pour des raisons de simplicité d'intégration et de manipulations (potentiel de dissipation dérivable notamment). Avoir une partie viscoplastique est aussi numériquement plus facile à formuler. Dans la suite, les paramètres de dissipation visqueuse $\bar{\nu}$ et $\hat{\nu}$ seront donc fixés à de faibles valeurs afin de faire numériquement tendre la visco-plasticité vers la plasticité. La valeur a été déterminée par simulations de traction simple en faisant varier $\bar{\nu}$. La Figure V.1 montre qu'en dessus de $10^{-3}$, la valeur n'a plus d'influence sur le comportement.

## V.A.1 Implémentation

L'implémentation dans Abaqus de ce modèle est réalisée selon la logique présentée dans l'Algorithme 3. Les grandeurs usuelles permettant le calcul de la contrainte de Cauchy sont d'abord définies, notamment le tenseur de Cauchy-Green gauche isochore à la fin de l'incrément et sa partie déviatorique.

Ensuite, on calcule les énergies de déformation et on crée une variable STATEV qui va permettre de stocker leur maximum et de l'actualiser éventuellement à chaque incrément. Cela permet de définir les variables $\bar{\eta}$ et $\hat{\eta}$ permettant la gestion de l'adoucissement lié à l'effet Mullins, respectivement sur les parties déviatorique et sphérique de la contrainte.

La résolution temporelle des lois d'évolution sur $\bar{\mathbf{B}}_e$ et $J_a$ nécessite d'avoir des variables contenant leurs valeurs en début et en fin d'incrément. On se sert à nouveau de variables d'état STATEV qui sont initialisées au premier incrément de sorte que $\bar{\mathbf{B}}_{e0} = \bar{\mathbf{B}}_e = \mathbf{1}$ et $J_{a0} = J_a = 1$ ($\bar{\mathbf{B}}_{e0}$ et $J_{a0}$ grandeurs en début d'incrément, $\bar{\mathbf{B}}_e$ et $J_a$ à la fin). Une condition fait en sorte qu'en début de chaque nouvel incrément la valeur de $\bar{\mathbf{B}}_e$ soit initialisée à $\bar{\mathbf{B}}_{e0}$ et $J_a$ à $J_{a0}$ (i.e. les valeurs convergées de l'incrément précédent).

La première équation différentielle est résolue par un schéma d'Euler explicite :

$$ 
\left| \begin{array}{l} 
\bar{\mathbf{B}}_{e|0} = \bar{\mathbf{B}}_{e0} \quad \text{pour } i = 0 \\ 
\bar{\mathbf{B}}_{e|_{i+1}} = \bar{\mathbf{B}}_{e|_i} + h \left[ -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{B}}_{e|_i} + \mathbf{L} \cdot \bar{\mathbf{B}}_{e|_i} + \bar{\mathbf{B}}_{e|_i} \cdot \mathbf{L}^T \right. \\ 
\left. \qquad -2\bar{\mathbf{B}}_{e|_i} \cdot \left[ \frac{\langle \|\bar{\sigma}_{neq}^D(\bar{\mathbf{B}}_{e|_i})\| - \hat{\chi}\rangle}{\bar{\nu}} \frac{\bar{\sigma}_{neq}^D(\bar{\mathbf{B}}_{e|_i})}{\|\bar{\sigma}_{neq}^D(\bar{\mathbf{B}}_{e|_i})\|} \right] \right] 
\end{array} \right. 
\quad (V.2) 
$$

avec $h$ le pas choisi pour les $n$ incréments constants de l'algorithme d'Euler subdivisant l'itération de Newton-Raphson en cours. L'équation différentielle de la deuxième loi prend une forme intégrable naturellement, s'écrivant en Euler explicite :

$$ 
\left| \begin{array}{l} 
J_{a|0} = J_{a0} \quad \text{pour } i = 0 \\ 
\ln(J_{a|_{i+1}}) = \ln(J_{a|_i}) + h \left[ \frac{\langle \|\hat{\sigma}_{neq}(J_{a|_i})\| - \hat{\chi}\rangle}{\hat{\nu}} \frac{(\hat{\sigma}_{neq}(J_{a|_i}) : \mathbf{1})}{\|\hat{\sigma}_{neq}(J_{a|_i})\|} \right] 
\end{array} \right. 
\quad (V.3) 
$$
**Algorithmme 3** Implémentation d'un opérateur tangent numérique dans une subroutine UMAT d'un comportement hyperélastique dissipatif

À l'incrément de temps $t_{n+1}$, Abaqus retourne la valeur du tenseur gradient de déformation $\mathbf{F}_{n+1}$. Pour simplifier la notation, on le notera simplement $\mathbf{F}$ et les grandeurs associées $\bar{\mathbf{B}}$, $\bar{\mathbf{B}}^D$, $\bar{\mathbf{B}}_e$, etc.

1. Calcul des tenseurs nécessaires au calcul de la contrainte à l'instant $t_{n+1}$ :
   $\bar{\mathbf{F}}, \bar{\mathbf{B}}$ (et ses invariants $I_{\bar{\mathbf{B}}} \text{ et } \mathcal{I}_{\bar{\mathbf{B}}}$) et $\bar{\mathbf{B}}^D$

2. Calcul des énergies libres isochore $\rho_0 \bar{\psi}$ et volumique $\rho_0 \hat{\psi}$

3. Calcul des variables $\bar{\eta}$ et $\hat{\eta}$ de l'effet Mullins

$$ \bar{\eta} = f(\bar{\psi}^D, \bar{\psi}_D^m) \quad \hat{\eta} = f(\hat{\psi}^D, \hat{\psi}_D^m) $$

4. Initialisation des variables utilisées pour la résolution des lois d'évolution :

Si l'incrément KINC est égal à 1 :

$$ \text{alors } \bar{\mathbf{B}}_{e0} = \bar{\mathbf{B}}_e = \mathbf{1} \text{ et } J_{a0} = J_a = 1 $$

Si passage à un incrément supérieur $n$ :

$$ \text{alors } \bar{\mathbf{B}}_{e0}^n = \bar{\mathbf{B}}_e^{n-1} \text{ et } J_{a0}^n = J_a^{n-1} $$

5. Résolution des équations différentielles des lois d'évolution déviatorique et sphérique (Euler explicite) :

Résolution de $\dot{J}_a=f(J_a)$ ⇒ Obtention de $J_a$

Résolution de $\dot{\bar{\mathbf{B}}}_e=f(\bar{\mathbf{B}}_e, J_a)$ ⇒ Obtention de $\bar{\mathbf{B}}_e$

6. Calcul du tenseur des contraintes de Cauchy STRESS :

$$ \boldsymbol{\sigma} = f(\bar{\mathbf{B}}^D, J, \bar{\mathbf{B}}_e, J_a, \bar{\eta}, \hat{\eta}) $$

7. Calcul des 6 gradients de déformation perturbés $\mathbf{F}^{(kl)}$, leur déterminant $J^{(kl)}$ et les 6 tenseurs de Cauchy-Green gauche isochores associés $\bar{\mathbf{B}}^{(kl)}$ et leurs invariants :

$$ \mathbf{F}^{(kl)} = \mathbf{F} + \frac{\varepsilon}{2} [\mathbf{e}_k \otimes \mathbf{e}_l \mathbf{F} + \mathbf{e}_l \otimes \mathbf{e}_k \mathbf{F}] $$

8. Calcul des énergies libres $\rho_0 \bar{\psi}^{(kl)}$ et $\rho_0 \hat{\psi}^{(kl)}$

9. Calcul des variables d'endommagement $\bar{\eta}^{(kl)}$ et $\hat{\eta}^{(kl)}$ de l'effet Mullins

10. Résolution de l'équation différentielle des lois d'évolution déviatorique et sphérique pour les 6 perturbations :

Résolution de $\dot{\bar{\mathbf{B}}}_e^{(kl)}=f(\bar{\mathbf{B}}_e^{(kl)})$ ⇒ Obtention de $\bar{\mathbf{B}}_e^{(kl)}$

Résolution de $\dot{J}_a^{(kl)}=f(J_a^{(kl)})$ ⇒ Obtention de $J_a^{(kl)}$

11. Calcul de la contrainte $\boldsymbol{\sigma}^{(kl)}$ :

$$ \boldsymbol{\sigma}^{(kl)} = f(\bar{\mathbf{B}}^{(kl)}, J^{(kl)}, \bar{\mathbf{B}}_e^{(kl)}, J_a^{(kl)}, \bar{\eta}^{(kl)}, \hat{\eta}^{(kl)}) $$

12. Calcul de l'opérateur tangent $\mathcal{C}_{ijkl}$ DDSDDE :

$$ \mathcal{C}_{ijkl} = \frac{1}{J_{\varepsilon}} [J^{(kl)} \sigma_{ij}^{(kl)}(\mathbf{F}^{(kl)}) - J \sigma_{ij}(\mathbf{F})] $$

13. Actualisation en fin d'incrément des variables STATEV

Une fois les équations différentielles résolues, la contrainte de Cauchy peut être calculée et retournée dans la variable STRESS.

Reste à fournir à Abaqus, l'opérateur tangent qui est calculé par la méthode numérique présentée dans la partie III.D.3. Toutes les grandeurs précitées sont calculées pour chacune des 6 perturbations du gradient de la transformation.

## V.A.2 Test du modèle

Le modèle a été testé avec des paramètres choisis arbitrairement pour dépasser le seuil de plasticité (voir Tableau V.1). Les simulations sont réalisées sur un élément de type C3D8
pour plusieurs types de sollicitations homogènes qui font travailler les parties déviatoriques et/ou sphériques de la contrainte de Cauchy. Les courbes contrainte-déformation obtenues sont tracées sur la Figure V.2. Plusieurs cyclages ont été réalisés :

- 5 cycles à amplitude de déformation croissante allant de 10 à 50%,

- 5 cycles à 50% de déformation.

Ils permettent de mettre en évidence le comportement attendu lié à l'effet Mullins sur des cycles croissants et aussi l'adoucissement de la contrainte menant à un cycle stabilisé lors de sollicitations multiples à la même amplitude.

<table><caption>TABLEAU V.1 – Paramètres choisis arbitrairement pour tester le modèle (visco-)élastoplastique (unités : SI). Les paramètres de viscosité $\bar{\nu}$ et $\hat{\nu}$ sont fixés à 0,001 (valeur à partir de laquelle ils n'ont plus d'influence et le modèle est élastoplastique).</caption><thead><tr><th>C<sub>10</sub></th><th>C<sub>01</sub></th><th>C<sub>30</sub></th><th>C<sub>p</sub></th><th>χ̄</th><th>J<sub>2</sub></th><th>s</th><th>q</th><th>κ</th><th>κ<sub>e</sub></th><th>χ̂</th><th>r̄</th><th>m̄</th><th>β̄</th><th>ρ̂</th><th>m̂</th><th>β̂</th></tr></thead><tbody><tr><td>-0,3</td><td>0,8</td><td>0,03</td><td>0,05</td><td>0,03</td><td>0,77</td><td>0,47</td><td>0,93</td><td>65</td><td>10</td><td>0,1</td><td>1,58</td><td>0,06</td><td>0,4</td><td>1,4</td><td>0,0008</td><td>0,7</td></tr></tbody></table>

Sur toutes les sollicitations considérées, la déformation résiduelle est bien présente grâce au choix des paramètres et en particulier des seuils de plasticité $\bar{\chi}$ et $\hat{\chi}$. La Figure V.2a est issue de simulations de traction uniaxiale. L'utilisation d'un modèle élasto-plastique couplé à l'effet Mullins a l'avantage de présenter de la dissipation sur les cycles après le premier chargement, comme le montre la courbe rouge. On remarque également que dans le cas de cycles à amplitude croissante, la plasticité permet de décaler légèrement les niveaux de contrainte entre une courbe de décharge et la charge suivante. C'est un effet que l'on observe expérimentalement et qu'il sera donc possible de modéliser.

En compressions hydrostatique et confinée, les mêmes observations peuvent être faites. Notons que pour obtenir un tel comportement en compression confinée, il est primordial de tenir compte de l'effet Mullins, mais également d'avoir de la (visco-)plasticité sur la partie sphérique de $\sigma$. En effet, la même simulation avec un comportement dissipatif portant uniquement sur la partie isochore (et donc une partie volumique hyperélastique) donnerait une courbe avec des cycles presque totalement superposés.

La présence de microsphères dans nos matériaux influe fortement sur le comportement hydrostatique et c'est sûrement pour cette même raison que le matériau disipe sous ce type de sollicitation (cette constatation a été faite lors de la caractérisation expérimentale). Les essais de compression hydrostatique au mercure que nous avons ne permettent pas de réaliser de « vrais chargements cyclés ». Nous disposons seulement de la variation de volume de l'éprouvette lors de l'intrusion et de l'extrusion de mercure. La courbe pression-variation de volume présente bien une hystérèse (voir la courbe rouge sur la Figure V.3), mais leur précision ne permet qu'une analyse qualitative et difficilement quantitative. Ainsi, la déformation résiduelle présente n'est
pas forcément représentative. Pour plus de précision, il semble indispensable de pouvoir réaliser un pilotage plus précis de l'essai (sur la décharge en particulier) et complexe (cyclage).

<figure><img src="image_77.png" /><figcaption>FIGURE V.2 – Exemples de comportements obtenus pour diverses sollicitations grâce au modèle élasto-(visco-)plastique implémenté.</figcaption></figure>

Sur la Figure V.3, on superpose à l'essai la courbe issue de la simulation de compression hydrostatique réalisée avec les paramètres du Tableau V.1. La plupart des paramètres sont conservés par rapport à ceux identifiés en hyperélasticité, donc la courbe bleue est assez proche du comportement observé expérimentalement. Cela montre bien que les prises en compte de la (visco-)plasticité et de l'effet Mullins semblent nécessaires pour modéliser au mieux le comportement de ces matériaux. La Figure V.2c présente d'autres résultats illustrant les capacités du modèle sous compression sphérique.
<figure><img src="image_78.png" /><figcaption>FIGURE V.3 – Courbes pression-variation de volume en compression hydrostatique. La courbe expérimentale en rouge représente l'intrusion et l'extrusion de mercure.</figcaption></figure>

Après avoir vérifié les capacités du modèle sur les chargements d'intérêts appliqués à un élément cubique linéaire, on s'intéresse à des sollicitations plus complexes comme le cas de compression/cisaillement ou traction/torsion. On maille pour cela notre cube avec 1000 éléments C3D8. Les déformées obtenues sont représentée sur la Figure V.4 selon l'évolution de l'élongation principale.

<figure><img src="image_79.png" /><figcaption>FIGURE V.4 – Exemples de simulations avec un maillage ayant plus d'un élément et sous chargements complexes. (a) Compression-Cisaillement. (b) Traction-torsion.</figcaption></figure>
## V.B Identification des paramètres

Il n'est pas possible d'identifier seulement les paramètres de plasticité en gardant ceux identifiés sur le modèle hyperélastique, car les modules des énergies libres $C_p$ et $\Lambda_e$ devraient rester petits pour que l'enveloppe de courbe soit identique. Ce choix conduirait à un non dépassement du seuil plastique en raison de cette faible valeur. Si nous voulons identifier des paramètres cohérents, il faut donc (ré)identifier tous les paramètres. Notre objectif est néanmoins de les identifier par « paquets » phénoménologiques, car il ne serait pas aisé de le faire pour 16 paramètres simultanément.

L'hypothèse qui a été faite dans le chapitre 4, de négliger la partie sphérique de la contrainte en traction uniaxiale est à nouveau adoptée. Ce faisant, il reste 8 paramètres à identifier : $C_{10}$, $C_{01}$, $C_{30}$, $C_p$, $\bar{\chi}$, $\bar{r}$, $\bar{m}$, $\bar{\beta}$.

Un plan d'expériences a été créé par la méthode de LHS. Comme précédemment, on garde les points dont la somme de $C_{10}$ et $C_{01}$ est positive. Notre plan d'expériences est constitué de 1558 points. On réalise des simulations de traction uniaxiale cyclée à amplitude de déformation croissante jusqu'à 50%. Les snapshots obtenus sont tracés sur la Figure V.5a avec le vecteur expérimental en croix noires.

Il n'est maintenant plus possible de visualiser simplement la variété de formes afin d'en déterminer sa dimension pour voir si tous les paramètres sont identifiables. Pour cela, nous allons nous servir de l'algorithme de Fukunaga et Olsen pour estimer la dimensionnalité intrinsèque locale des points issus de la POD. Il est déjà implémenté dans le package Python « scikit-dimension » (BAC et al. 2021). L'utilisation de cette fonction retourne pour chacun des 1558 points leur dimensionnalité intrinsèque, calculée en fonction de leur voisinage. La moyenne de cette dimensionnalité est d'environ 6. Cela veut dire que nous pourrons, a priori, identifier 6 paramètres sur les 8 du modèle.

L'interpolation par krigeage a quand même été menée afin d'obtenir notre modèle réduit et le snapshot expérimental projeté sur la variété de formes pour voir les résultats donnés par l'identification. Souvent, quand un paramètre n'est pas identifiable, ou quand il ne varie pas dans une plage pertinente, l'algorithme d'optimisation retourne des valeurs égales aux valeurs initiales ou aux bornes de recherche définies en amont. C'est effectivement le cas dans notre identification, où 2 des paramètres semblent bien non identifiables : $\bar{\chi}$ a une valeur trop grande, rendant ainsi le modèle hyperélastique, et $\bar{r}$ est égal à la borne inférieure de l'intervalle de recherche. Les autres paramètres ont quant à eux des valeurs auxquelles on pouvait s'attendre.

<table><caption>TABLEAU V.2 – Paramètres identifiés en traction uniaxiale cyclée.</caption><tr><td>$C_{10}$</td><td>$C_{01}$</td><td>$C_{30}$</td><td>$C_p$</td><td>$\bar{\chi}$</td><td>$\bar{r}$</td><td>$\bar{m}$</td><td>$\bar{\beta}$</td></tr><tr><td>-0,33</td><td>0,8</td><td>0,029</td><td>0,04</td><td>0,02</td><td>1,58</td><td>0,06</td><td>0,4</td></tr></table>

Pour avoir des paramètres représentant un peu mieux le comportement en traction du R23, nous avons, utilisant notre connaissance du modèle et par retour d'expérience, baissé la valeur
<figure><img src="image_80.png" /></figure>

<figure><img src="image_81.png" /></figure>

<figure><img src="image_82.png" /><figcaption>FIGURE V.5 – Identification des paramètres déviatoriques sur des essais de traction uniaxiale.</figcaption></figure>

de $\bar{\chi}$ impliquant des ajustements des paramètres $C_{ij}$ dans un voisinage des valeurs identifiées. Quant à $\bar{r}$, nous l'avons choisi égal à celui identifié sur le modèle hyperélastique. Les paramètres retenus sont listés dans la Tableau V.2. Le comportement obtenu avec ces paramètres est tracé sur les Figures V.5b et V.5c. La contrainte des deux premiers cycles est bien reproduite, puis pour des niveaux de déformation supérieurs, les contraintes restent très proches et la forme des cycles plutôt bien respectée. Le point négatif est le manque de déformation résiduelle en pied de courbe. Un modèle élasto-plastique seul (ou en tout cas notre formulation) n'est sûrement pas capable d'introduire autant de plasticité, il faudrait également prendre en compte la viscosité. Expérimentalement, elle doit également influer sur le comportement, et une part de la déformation à contrainte nulle, considérée jusqu'à présent comme permanente, pourrait finalement présenter une part de déformation rémanente liée à un ou plusieurs phénomènes visqueux de temps caractéristiques plus ou moins importants.
## REMARQUE V.1

La même démarche a été effectuée avec le R33, cependant, dans le cas de ce matériau, seulement 5 des 8 paramètres étaient identifiables et avaient des valeurs cohérentes. Il aurait été trop difficile de trouver un comportement contrainte-elongation adéquat en ajustant les paramètres.

Suite à l'identification proposée du R23 (cf. Tableau V.2), seule la partie déviatorique est déterminée. Afin d'identifier le modèle dans son intégralité, il nous faudrait identifier simultanément les paramètres sur la compression confinée et la traction simple. Seulement les inconvénients sont nombreux :

* Les paramètres à identifier auraient été au nombre de 16 :

- le plan d'expériences nécessaire devrait contenir un très grand nombre de points,

- et il y aurait de fortes chances qu'un grand nombre d'entre eux ne soient pas identifiables.

* Les essais de compression confinée dont nous disposons pourraient ne pas uniquement refléter le comportement des matériaux. En effet, la forme des cycles et la présence d'une importante déformation résiduelle laissent penser qu'en plus de la plasticité et d'une forte viscosité inhérente au matériau, l'influence du montage d'essai est prépondérante. Le frottement répété entre les parois de la cavité et l'éprouvette peut être à l'origine de ce comportement (voir Figure V.6).

<figure><img src="image_83.png" /><figcaption>FIGURE V.6 – Courbe contrainte-variation de volume en compression confinée. La courbe bleue est obtenue avec le modèle (visco-)élastoplastique pour des paramètres de la partie sphérique non optimisés.</figcaption></figure>

Pour ces raisons, il ne semble pas judicieux d'essayer d'identifier des paramètres sur cet essai. D'autant plus que le modèle ne sera sûrement pas capable de reproduire un tel comportement.
La courbe bleue sur la Figure V.6 est obtenue avec le modèle pour des paramètres proches de ceux identifiés en hyperélasticité. Le seuil de plasticité, fixé assez bas, permet effectivement de faire apparaître un peu de déformation résiduelle, mais qui reste assez modeste en comparaison de l'attendu expérimental. Pour s'affranchir de ces problèmes, il faudrait mettre un accent particulier sur le développement d'essais qui décrivent uniquement le comportement des matériaux et dont le type de chargement soit assuré. Un essai de compression sphérique cyclé serait l'idéal pour identifier les paramètres liés à la partie volumique du comportement. Cela permettrait même d'identifier les paramètres de chaque partie de la contrainte par étapes et d'identifier moins de paramètres simultanément. Nous pourrions également coupler des essais purement déviatoriques (cisaillement simple) ou sphériques (compression hydrostatique) avec des essais comme la traction uniaxiale. Le modèle de krigeage serait bien plus riche en informations pour permettre une meilleure identification des paramètres.

## V.C Possibilités du modèle élasto-(visco-)plastique

La viscosité que nous avions placée dans les phénomènes secondaires à modéliser en conclusion du chapitre 2, semble, comme évoqué dans la section précédente, avoir un rôle à jouer dans la modélisation du comportement des matériaux à plus ou moins long temps caractéristique. De plus, dans leur application, les joints peuvent subir des sollicitations longues comme le fluage ou la relaxation, mettant en jeu la réponse temporelle du matériau.

### REMARQUE V.2

Dans les simulations proposées dans cette section et la suivante (i.e. correspondant aux résultats donnés Figures V.8 et V.9), il a été choisi de ne plus considérer les deux paramètres de viscosité comme de « simples potentiels » de régularisation du comportement élasto-plastique, mais de leur donner des significations physiques.

### V.C.1 Viscosité à temps caractéristique court

#### V.C.1.a Relaxation

La Figure V.7 présente quelques essais de relaxation, en traction (durée 1 heure) et en compression confinée (durée 20 min) sur le R23. La vitesse de déformation avant le palier est de 100 mm/min. On observe bien une relaxation de la contrainte. En faisant varier un peu $C_p$ et en augmentant le paramètre de viscosité $\bar{\nu}$ avec le reste des paramètres identifiés dans la section précédente, on obtient une simulation de cet essai de relaxation présenté sur la Figure V.8.

La valeur de ces deux paramètres est déterminante pour modéliser la relaxation de la contrainte. L'ajout de cet essai aux snapshots pourrait être avantageux pour donner plus d'importance à $C_p$ et améliorer son identifiabilité.

La relaxation en compression confinée sur la Figure V.8b présente le même type de résultats, les paramètres $\kappa_e$ et $\bar{\nu}$ influencent beaucoup le comportement.
<figure><img src="image_84.png" /><figcaption>FIGURE V.7 – Essais de relaxation sur le R23. (a) Relaxation en traction pour deux niveaux de déformation. (b) Relaxation en compression confinée.</figcaption></figure>

<figure><img src="image_85.png" /><figcaption>FIGURE V.8 – Simulations de relaxation sur le R23. (a) Relaxation en traction. (b) Relaxation en compression confinée.</figcaption></figure>

### V.C.1.b Influence de la vitesse de sollicitation

Dans la partie II.B.4.c, nous avions tracé l'évolution de la contrainte de Piola-Kirchhoff I, en traction, compression confinée et cisaillement simple, pour des vitesses de déformation initiale différentes. L'influence sur le comportement n'étant pas très marquée, nous avions d'abord négligé cet aspect dans la modélisation. Grâce au modèle élasto-(visco-)plastique, en fixant $\bar{\eta}$ à 14 MPa.s comme dans le paragraphe précédent, l'influence de la vitesse de sollicitation peut être modélisée. La rigidification du matériau à vitesse de déformation croissante est mise en évidence. On remarque que l'influence apparaît bien à partir du moment où le matériau plastifie. Avant cela les courbes sont superposées. Expérimentalement, les essais de traction ont aussi ce comportement.
<figure><img src="image_86.png" /><figcaption>FIGURE V.9 – Simulations de traction uniaxiale cyclée. La vitesse de déformation initiale varie entre $10^{-3}$ et $10^{-1}$ s$^{-1}$.</figcaption></figure>

## V.C.2 Viscosité à temps caractéristique long

Le CEA étudie également le comportement de recouvrement du R23 après du fluage en compression confinée. Une large campagne d'essais a été menée visant à estimer la recouvrance de volume d'éprouvettes ayant subi du fluage pendant des temps variables et à différentes températures et pressions. Le schéma à gauche de la Figure V.10 décrit le montage permettant de réaliser le fluage.

Une éprouvette cylindrique est placée dans une cavité et une pièce de la même taille que la

<figure><img src="image_87.png" /><figcaption>FIGURE V.10 – Essais de relaxation sur le R23. (a) Relaxation en traction pour deux niveaux de déformation. (b) Relaxation en compression confinée.</figcaption></figure>

cavité la comprime avec une force appliquée par un ressort (la rigidité calibrée de celui-ci permet
par application d'un déplacement d'imposer une force constante). Ce dernier permet d'atteindre des pressions de 5 ou 15 bars. Il y a de nombreux systèmes comme celui-ci fixés sur un montage placé dans une enceinte thermique pour faire des essais à 60°C. La durée de fluage varie entre 1 et 11 mois. Aux temps souhaités, une éprouvette est sortie de la cavité, par exemple après 3 mois de fluage, puis stockée à température ambiante. Sa densité est mesurée (par la technique de poussée d'Archimède) sur une plage de temps allant de 1h à 1 an.

La Figure V.10 résume schématiquement le principe de la campagne d'essais.

Quelques résultats de cette campagne sont tracés sur la Figure V.11. La variation de densité est représentée en fonction du temps de recouvrance après fluage. Le constat global est que cette variation diminue avec le temps signe que le volume de l'éprouvette augmente après le relâchement de la contrainte, mais il tend à se stabiliser. La recouvrance n'est pas totale même au bout d'un an. C'est le signe d'une viscosité à plus long terme et possiblement également d'une déformation plastique permanente.

<figure><img src="image_88.png" /></figure>

<figure><img src="image_89.png" /><figcaption>FIGURE V.11 – Évolution de la variation de densité après des essais de fluage en compression confinée sur le R23. (a) Comparaison de l'effet du temps de fluage. (b) Comparaison de l'effet de la pression de fluage.</figcaption></figure>

En l'état, notre modèle ne pourrait sûrement pas reproduire un tel comportement, d'autant plus que pour « activer » le caractère visqueux, il faut avoir atteint le seuil de plasticité. Or, une sollicitation de 15 bars (soit 0,15 MPa) en compression confinée ne permettrait pas à la contrainte $\hat{\sigma}_{neq}$ d'avoir des valeurs suffisamment élevées pour plastifier. De plus, le temps caractéristique de la branche visco-plastique pourrait être trop court pour des phénomènes ayant lieu sur une voire plusieurs années. Dans cette optique, nous proposons une voie d'amélioration du modèle permettant d'apporter de la viscosité pour différents temps caractéristiques. Cela consiste à ajouter autant d'états intermédiaires que nécessaire à la transformation (voir Figure V.12). La décomposition du gradient de déformation s'écrit alors :

$$F = \hat{F} \cdot \bar{F} \quad \text{avec} \quad \bar{F} = \bar{F}^{(i)} \cdot {}^{(i)}\bar{F} \quad \text{et} \quad \hat{F} = \hat{F}^{(i)} \cdot {}^{(i)}\hat{F} \quad \forall i \qquad (V.4)$$
$\hat{\mathbf{F}}^{(i)}$ et $^{(i)}\hat{\mathbf{F}}$ dépendent par définition, respectivement de, $J^{(i)}$ et $^{(i)}J$.

Les tenseurs de Cauchy-Green gauche et droit associés à ces gradients s'écrivent :

$$ \bar{\mathbf{B}}^{(i)} = \bar{\mathbf{F}}^{(i)} \cdot \bar{\mathbf{F}}^{(i)T} \qquad \bar{\mathbf{C}}^{(i)} = \bar{\mathbf{F}}^{(i)T} \cdot \bar{\mathbf{F}}^{(i)} \qquad (V.5) $$

$$ ^{(i)}\bar{\mathbf{B}} = ^{(i)}\bar{\mathbf{F}} \cdot ^{(i)}\bar{\mathbf{F}}^{T} \qquad ^{(i)}\bar{\mathbf{C}} = ^{(i)}\bar{\mathbf{F}}^{T} \cdot ^{(i)}\bar{\mathbf{F}} \qquad (V.6) $$

Les tenseurs taux de déformation sont dans ce cas :

$$ ^{(i)}\bar{\mathbf{D}} = (^{(i)}\dot{\mathbf{F}} \cdot ^{(i)}\bar{\mathbf{F}}^{-1})_{\text{sym}} \qquad \bar{\mathbf{D}}^{(i)} = (\dot{\bar{\mathbf{F}}}^{(i)} \cdot \bar{\mathbf{F}}^{(i)-1}) \qquad (V.7) $$

Le tenseur objectif du taux de déformation analogue de $\bar{\mathbf{D}}_a^o$ s'écrit :

$$ ^{(i)}\bar{\mathbf{D}}^{\circ} = \bar{\mathbf{R}}^{(i)} \cdot ^{(i)}\bar{\mathbf{D}} \cdot \bar{\mathbf{R}}^{(i)-1} \qquad (V.8) $$

<figure><img src="image_90.png" /><figcaption>FIGURE V.12 – Décomposition de la transformation en parties isochore et volumique. Ces dernières peuvent ensuite être décomposées en une infinité d'état intermédiaires permettant la prise en compte de phénomènes dissipatifs.</figcaption></figure>
L'énergie libre se décompose dans ce cas comme suit :

$$ \psi = \psi(J, J^{(i)}, \bar{\mathbf{B}}, \bar{\mathbf{B}}^{(i)}) = \hat{\psi}_{eq}(J) + \bar{\psi}_{eq}(\bar{\mathbf{B}}) + \sum_{i=1}^{n} \hat{\psi}_{neq}^{(i)}(J^{(i)}) + \sum_{i=1}^{n} \bar{\psi}_{neq}^{(i)}(\bar{\mathbf{B}}^{(i)}) \quad (V.9) $$

Après son insertion dans l'inégalité de dissipation Clausius-Duhem, et après avoir appliqué le cadre des matériaux standards généralisés et le principe de normalité de la dissipation pour un potentiel se décomposant de la sorte,

$$ \varphi = \bar{\varphi}_0(\bar{\mathbf{D}}) + \hat{\varphi}_0(\hat{\mathbf{D}}) + \sum_{i=1}^{n} \hat{\varphi}^{(i)}(\bar{\mathbf{D}}^{\circ}) + \sum_{i=1}^{n} \hat{\varphi}^{(i)}(\hat{\mathbf{D}}^{\circ}) \quad (V.10) $$

on obtient une nouvelle écriture de la loi de comportement et des lois d'évolution. On veut que le modèle soit élasto-visco-plastique, on choisit donc d'avoir un état intermédiaire visco-plastique pour les parties isochore et volumique de la déformation. Afin d'introduire plusieurs temps caractéristiques de viscosité, nous pouvons utiliser $(n-1)$ intermédiaires visco-élastiques, comme présenté sur les équations du modèle (V.11).

$$ \left\{
\begin{array}{l}
\sigma = \bar{\sigma}_{eq} : \mathbb{P}_1 + \sum_{i=1}^{n} \bar{\sigma}_{neq}^{(i)} : \mathbb{P}_1 + \hat{\sigma}_{eq} : \mathbb{P}_1 + \sum_{i=1}^{n} \hat{\sigma}_{neq}^{(i)} : \mathbb{P}_1 \\
\\
\quad {}^1\bar{\mathbf{D}}^{\circ} = \frac{\partial^{(1)}\bar{\varphi}}{\partial \bar{\sigma}_{neq}^{(1)}} \\
\\
\quad \hat{\sigma}_{neq}^{(i)} : \mathbb{P}_1 = \frac{\partial^{(i)}\bar{\varphi}}{\partial \hat{\mathbf{D}}^{\circ}} \quad \forall i \in (2; n) \hfill (V.11) \\
\\
\quad {}^1\hat{\mathbf{D}}^{\circ} = \frac{\partial^{(1)}\hat{\varphi}}{\partial \hat{\sigma}_{neq}^{(1)}} \\
\\
\quad \hat{\sigma}_{neq}^{(i)} \mathbb{1} = \frac{\partial^{(i)}\hat{\varphi}}{\partial \hat{\mathbf{D}}^{\circ}} \quad \forall i \in (2; n)
\end{array}
\right.
$$

### REMARQUE V.3

Le temps caractéristique de la viscosité dans le modèle est défini par le rapport entre les paramètres de rigidités et de viscosités tel que :

$$ \bar{\tau}^{(i)} = \frac{\bar{\nu}^{(i)}}{C_p^{(i)}} \quad \text{pour la partie isochore et pour la partie volumique } \hat{\tau}^{(i)} = \frac{\hat{\nu}^{(i)}}{\kappa^{(i)}} \quad (V.12) $$

On pourrait imaginer un modèle où les temps caractéristiques courts seraient pris en charge par la partie visco-plastique du modèle et à plus long terme, une partie visco-élastique serait utilisée. Cela donne, par exemple, les équations suivantes, où les lois d'évolutions sont écrites
sous forme différentielle :

$$
\left\{
\begin{aligned}
\bar{\sigma} &= 2J^{-1}\bar{\eta} \left[ (C_{10} + 3C_{30}(I_{\bar{B}} - 3)^2)\bar{B} - C_{01}\bar{B}^{-1} \right] : \mathbb{P}_1 + 2J^{-1} \left[ C_p^{(1)}(I_{\bar{B}^{(1)}} - 3) \right] : \mathbb{P}_1 \\
&\quad + 2J^{-1} \left[ C_p^{(i)}(I_{\bar{B}^{(i)}} - 3) \right] : \mathbb{P}_1 \quad \forall i \in (2; n) \\[1em]
\hat{\sigma} &= \hat{\eta}\kappa \left[ (1-q)a \tan\left(\frac{J-1}{a}\right) + qb \tanh\left(\frac{J-1}{b}\right) \right] \mathbf{1} + {}^{(1)}J^{-1}\kappa^{(1)} \left[ a \tan\left(\frac{J^{(1)}-1}{a}\right) \right] \mathbf{1} \\
&\quad + {}^{(i)}J^{-1}\kappa^{(i)} \left[ a \tan\left(\frac{J^{(i)}-1}{a}\right) \right] \mathbf{1} \quad \forall i \in (2; n) \\[1em]
\dot{\bar{B}}^{(1)} &= -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{B}^{(1)} + \mathbf{L} \cdot \bar{B}^{(1)} + \bar{B}^{(1)} \cdot \mathbf{L}^T - 2\bar{B}^{(1)} \cdot \left[ \frac{\langle \|\bar{\sigma}_{neq}^{(1)D}\| - \bar{\chi}^{(1)} \rangle}{\bar{\nu}^{(1)}} \frac{\bar{\sigma}_{neq}^{(1)D}}{\|\bar{\sigma}_{neq}^{(1)D}\|} \right] \\[1em]
\dot{\bar{B}}^{(i)} &= -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{B}^{(i)} + \mathbf{L} \cdot \bar{B}^{(i)} + \bar{B}^{(i)} \cdot \mathbf{L}^T - 2\bar{B}^{(i)} \cdot \frac{\bar{\sigma}_{neq}^{(i)D}}{\bar{\nu}^{(i)}} \quad \forall i \in (2; n) \\[1em]
{}^{(1)}J^{-1(1)}\dot{j} &= \frac{\langle \|\hat{\sigma}_{neq}^{(1)}\| - \hat{\chi}^{(1)} \rangle}{\hat{\nu}^{(1)}} \frac{(\hat{\sigma}_{neq}^{(1)}:\mathbf{1})}{\|\hat{\sigma}_{neq}^{(1)}\|} \\[1em]
{}^{(i)}J^{-1(i)}\dot{j} &= \frac{1}{\hat{\nu}^{(i)}} (\hat{\sigma}_{neq}^{(i)}:\mathbf{1}) \quad \forall i \in (2; n)
\end{aligned}
\right.
\qquad (V.13)
$$

L'identification des paramètres de ce modèle nécessite la réalisation d'essais de relaxation pour identifier des temps caractéristiques $\tau_i$ cohérents. Ils sont possiblement différents selon le type de sollicitation, donc il faudrait peut-être introduire plusieurs états intermédiaires. Un essai de compression hydrostatique est également nécessaire pour apporter plus d'informations sur la partie volumique dissipative de la loi de comportement.

Des essais de relaxation par paliers pourraient également s'avérer intéressants, car ils permettent de reconstruire une courbe de charge/décharge en s'affranchissant des phénomènes visqueux. Théoriquement, si la charge et la décharge sont superposées, alors le matériau est visco-élastique. Si on observe une hystérèse, le matériau présente sans doute aussi de la plasticité. Cela permettrait de mieux cerner le comportement de nos matériaux et orienterait le choix des états intermédiaires pour le modèle.

Après l'identification de ses paramètres, ce modèle pourrait sans doute décrire le comportement en fluage et recouvrance de nos matériaux à température ambiante. La Figure V.11 montre que la température à laquelle le fluage est effectué influence assez fortement la variation de densité. L'éprouvette se déformerait davantage lors d'une sollicitation à 60°C qu'à 20°C. En plus d'être influencées par le niveau de charge et la durée de fluage, la déformation post-fluage et la recouvrance de volume de l'éprouvette semblent dépendantes de la température. Étant donné les temps de maintien de la pression à 60°C, le matériau doit également subir un vieillissement et donc un changement dans sa structure moléculaire. Ce sont plusieurs phénomènes qu'un modèle phénoménologique purement mécanique ne pourra pas intégrer.
De plus, ces phénomènes sont liés. La température influence sûrement la viscosité du matériau et son état chimique. Si la structure est modifiée sous l'effet de la température, cela doit égale-
ment impacter les propriétés mécaniques. Tous ces phénomènes seraient alors modélisables en considérant, un couplage entre la thermique et la mécanique, mais également avec la chimie.
# Conclusions du chapitre

Ce chapitre a présenté l'implémentation d'un modèle élasto-(visco-)plastique par une routine UMAT dans le code de calculs Abaqus. La visco-plasticité est numériquement plus facile à gérer, c'est pourquoi dans un premier temps on donne une faible valeur au terme de viscosité. On a alors un modèle avec deux lois d'évolution, chacune permettant d'introduire une composante plastique (i.e. une sur la partie volumique et une sur son homologue isochore). Des exemples de simulations cyclées en traction, compressions confinée et hydrostatique, et cisaillement permettent de vérifier les réponses attendues d'un point de vue phénoménologique. Les résultats montrent une bonne prise en compte du comportement recherché :

* une déformation résiduelle est présente,

* l'effet Mullins est pris en compte sur des cycles à amplitude croissante,

* l'association de la plasticité et de l'effet Mullins permet de modéliser un comportement de cycle stabilisé lors de multiples sollicitations à amplitude de déformation constante,

* après dépassement du seuil de plasticité, il y a un décalage entre les courbes de décharge et de charge suivante.

La deuxième partie du chapitre a montré qu'avec les essais mécaniques à notre disposition, une seule partie des paramètres recherchés pouvait être identifiée. En effet, les paramètres de la partie déviatorique ont pu être calibrés uniquement sur des essais de traction uniaxiale. Le couplage avec la compression confinée comme dans le chapitre 4 n'est pas possible, car nous devrions alors identifier un nombre trop grand de paramètres (16). L'utilisation de la POD sur l'essai de traction et la dimension de la variété de courbe obtenue, inférieure au nombre de paramètres recherchés, indique que l'unicité du résultat d'identification n'est pas garantie. Le jeu de paramètres légèrement ajusté donne néanmoins des résultats corrects en traction. Il n'est actuellement pas possible d'identifier les paramètres liés à l'énergie libre volumique, car nous manquons d'un essai de compression hydrostatique cyclée.

Comme évoqué en conclusion du chapitre 2, un des phénomènes secondaires à modéliser est la viscosité, jusqu'à présent très peu influente (a priori) dans les cas de chargement appréhendés et donc négligée. Néanmoins quelques essais préliminaires montrent que nos matériaux peuvent relaxer à des temps caractéristiques assez courts et fluer dans le cas de sollicitations longs termes proches de celles rencontrées dans le domaine d'application. Nous nous sommes alors intéressés à la prise en compte d'une dépendance temporelle de la réponse du modèle proposé. Dans un premier temps, on donne une signification physique à $\bar{\nu}$ et $\hat{\nu}$ (jusqu'alors uniquement utilisés comme paramètres de régularisation) et on augmente leur valeur. Cela permet :

* d'observer une relaxation des contraintes en traction comme en compression confinée,

* d'apporter une influence de la vitesse de déformation sur le comportement mécanique.

La Tableau V.3 présente le jeu de paramètres de ce modèle, obtenu après identification et ajustage, représentant à ce jour au mieux le comportement du matériau R23 en traction (cf. Figure V.5) et en compression confinée. Pour cette dernière sollicitation les paramètres de la partie sphérique ont été ajustés pour refléter le comportement de première charge, tout en
ayant un peu de viscosité et un seuil de plasticité bas pour induire une déformation résiduelle. Le comportement cyclé ne pouvant pas être identifié sur nos essais faute de représentativité, les paramètres permettant le calcul de $\hat{\eta}$ sont choisis identiques à la partie déviatorique. Les résultats issus de ce choix sont tracés sur la Figure V.13.

Dans le cas du R33, les résultats d'identification ne sont pas assez satisfaisants car il permettent d'obtenir un nombre trop faible de paramètres.

<table><caption>TABLEAU V.3 – Jeu de paramètres du R23 proposé pour le modèle (visco-)élastoplastique (unités : SI)</caption><thead><tr><th>C<sub>10</sub></th><th>C<sub>01</sub></th><th>C<sub>30</sub></th><th>C<sub>p</sub></th><th>χ̄</th><th>ν̄</th><th>J<sub>2</sub></th><th>s</th><th>q</th><th>κ</th><th>κ<sub>e</sub></th><th>χ̂</th><th>ν̂</th><th>r̄ = r̂</th><th>m̄ = m̂</th><th>β̄ = β̂</th></tr></thead><tbody><tr><td>-0,33</td><td>0,8</td><td>0,029</td><td>0,04</td><td>0,02</td><td>0,001</td><td>0,77</td><td>0,47</td><td>0,93</td><td>30</td><td>3</td><td>0,75</td><td>1,2</td><td>1,58</td><td>0,06</td><td>0,4</td></tr></tbody></table>

<figure><img src="image_91.png" /><figcaption>FIGURE V.13 – Courbe contrainte-variation de volume en compression confinée (R23) pour les paramètres du Tableau V.3.</figcaption></figure>

L'utilisation de la partie visco-plastique permet d'apporter seulement un temps visqueux caractéristique. Or les matériaux peuvent également connaître des phénomènes à plus long terme. En effet, lorsqu'ils sont sollicités longuement en fluage, puis que la charge est retirée, les matériaux recouvrent petit à petit leur volume (sans atteindre le volume de départ), sur un temps s'étalant de plusieurs mois à plusieurs années. La durée (aussi potentiellement très longue) de chargement semble aussi influer sur les réponses.

Afin de rendre compte de ces différents temps caractéristiques, nous proposons une amélioration du modèle par l'introduction d'états intermédiaires multiples. Il pourrait y avoir une partie (visco-)plastique (avec viscosité plus ou moins négligeable), puis autant d'états intermédiaires visco-élastiques que de phénomènes visqueux à prendre en compte.

Afin d'identifier les paramètres de ce modèle, il faudrait réaliser des essais supplémentaires,
notamment des essais de relaxation. L'identification pourrait être réalisée en simultané, pour plusieurs types de sollicitations, avec des chargements cyclés et de la relaxation. Ainsi, les paramètres de la partie hyperélastique auront plus d'influence sur les essais charge/décharge et les paramètres visco-plastique sur la relaxation. Globalement, le modèle de krigeage pourrait donc être plus riche et permettre l'identification de plus de paramètres. Il apparaît également nécessaire d'avoir des essais de compression hydrostatique (cyclage et relaxation) afin de pouvoir découpler un minimum l'identification sans avoir un nombre trop important de paramètres en simultané.

Point sur les phénomènes pris en compte après le Chapitre 5

■ Phénomènes de 1<sup>er</sup> ordre :

✔ grandes déformations et non-linéarités,

✔ grande compressibilité,

✔ effet Mullins,

✔ déformation résiduelle,

>D influence de la fraction volumique de microsphères.

■ Phénomènes de 2<sup>nd</sup> ordre :

X effet de la température,

>D effet visqueux (Fluage, relaxation).
Conclusions et perspectives

Le CEA conçoit et développe des structures multi-matériaux nécessitant l'utilisation de liaisons jointives compressibles entre certaines pièces pour éviter la concentration de contraintes et ainsi prévenir la rupture des parties fragiles lors de cycles thermiques. Les matériaux choisis pour ces éléments de jonction sont des élastomères silicones vulcanisés à température ambiante (RTV) auxquels on ajoute des microsphères à coques thermoplastiques expansées à l'isobutane. Cela permet d'induire une compressibilité importante et recherchée à un matériau initialement incompressible (ou quasiment) par un procédé de fabrication simple qui nécessite un mélange des composants puis une réticulation dans un moule sans apport de chaleur. Dans cette thèse, trois grades de ce matériau sont étudiés, le RTV seul, le RTV avec 23% de microsphères (R23) et celui avec 33% de microsphères (R33).

Le silicone RTV présente des caractéristiques classiques d'un élastomère, à savoir, un caractère quasi-incompressible, une capacité à se déformer de manière hyperélastique sur une large plage de déformation, ainsi que la présence d'une faible viscosité. On observe également un léger adoucissement de la contrainte sur les premiers cycles de charge. Ce phénomène, l'effet Mullins, est souvent imputé à l'interaction des macromolécules avec les charges renforçantes ajoutées à la matrice. L'ajout de microsphères à ce matériau confère une véritable originalité du comportement thermo-mécanique des mélanges, en changeant, en plus de la compressibilité initiale du RTV, les phénomènes mis en jeu. Ainsi, le comportement observé en traction et cisaillement simple semble indiquer un effet de charges renforçantes des microsphères, la contrainte étant initialement plus élevée que pour le silicone pur. Lors de chargements triangulaires cyclés, l'adoucissement de la contrainte après le premier cycle est plus marqué que pour le RTV seul et pour des gammes de déformations inférieures. Cette prématurité renforce l'« effet charge » des microsphères évoqué précédemment. Dans l'hypothèse où l'effet Mullins est accentué par le taux de charges, les sphères doivent avoir localement une forte interaction avec le RTV. De plus cet effet est présent non seulement en cisaillement et en traction, mais également en compression confinée qui est quant à elle principalement sphérique.

Déformé en compressions hydrostatique ou, à défaut, confinée, les R23 et R33 ont des réponses fortement non-linéaires, sans aucun doute induites par le comportement des microsphères. Ces dernières vont en effet se déformer d'abord de manière élastique, puis flamber sous la contrainte pour enfin se fermer complètement induisant ainsi la limite de compressibilité du matériau atteinte à une déformation proche du pourcentage de microsphères dans le matériau (on retrouve
alors le comportement du RTV seul).

Il est notable de rappeler que tous ces phénomènes anélastiques (adoucissement, déformation résiduelle, impact de la vitesse de sollicitation) sont présents pour tout type de chargement y compris sphériques (ou quasiment). L'objectif du CEA est de disposer d'une loi de comportement représentative de ces phénomènes afin de simuler convenablement les structures développées dans leurs applications.

Peu de modèles existent dans la littérature tenant compte d'un caractère fortement compres-
sible des élastomères (puisque peu de matériaux exhibent une telle propriété) et encore moins
pour des cas de trajets de chargements complexes et variés. La plupart des lois de comporte-
ment sont développées pour des élastomères supposés quasi-incompressibles et permettent de
rendre compte de leurs grandes déformations et leurs comportements dissipatifs (adoucissement
compris).

Nous allons repartir de cette base et la développer pour l'étendre au comportement original de
nos matériaux. Dans le cadre des grandes déformations, nous adoptons en premier lieu l'hypo-
thèse majeure de FLORY 1961 qui divise la transformation en deux parties, une comprenant les
déformations à volume constant (isochore) et l'autre les changements de volume (volumique). La
prise en compte de phénomènes dissipatifs se fait par l'introduction d'états intermédiaires, no-
tion initialement utilisée par GREEN & NAGHDI 1965 et SIDOROFF 1973 pour l'élasto-plasticité.
Le cadre des matériaux standards généralisés permet d'introduire deux états visco-plastiques,
respectivement sur chaque contribution de la transformation. De même, deux variables per-
mettent de piloter l'adoucissement des contraintes via le modèle d'OGDEN & ROXBURGH 1999.
On obtient alors une loi de comportement et ses deux lois d'évolutions pilotant les deux va-
riables internes définissant les deux transformations visco-plastiques (les deux lois d'évolution
des deux paramètres d'adoucissement étant quant à elles pilotés indépendamment). Cette pre-
mière hypothèse forte a d'ailleurs, pour principal effet, le découplage de la contrainte en parties
sphérique et déviatorique, conséquence certes forte, mais s'avérant très utile dans la suite.

Ainsi, suivant cette logique, nous avons pu implémenter dans le code de calculs Abaqus (via
des routines utilisateurs de type UMAT) :

■ Un modèle hyperélastique compressible présentant de l'effet Mullins sur la transformation isochore. Pour cela, l'opérateur tangent est développé analytiquement.

- Un modèle hyperélastique compressible ajoutant à la précédente proposition une prise en compte de l'effet Mullins sur la partie volumique du gradient de transformation. Dans ce cas, l'opérateur tangent est calculé numériquement (par perturbation du gradient) pour simplifier la démarche. Il est à noter que bien qu'adoptée pour des raisons pragmatiques, cette approche numérique donne une véritable souplesse de développement en évitant le calcul long et fastidieux de l'opérateur tangent de manière (semi-)analytique.

- Le modèle élasto-(visco)-plastique complet (i.e. formalisant un apport visco-plastique et un adoucissement de type Mullins sur les deux parties de la transformation). Pour plus de liberté de formulation, on continue à utiliser la méthode de l'opérateur tangent numérique. Dans un premier temps, les paramètres de viscosité sont fixés petits et servent de paramètres de régularisation des deux contributions élasto-plastiques.
Naturellement, comme souvent dans le cas d'approches phénoménologiques telles que celle adoptée dans le présent document, tous nos modèles présentent un certain nombre de paramètres à identifier. Pour cela, nous employons une approche innovante, différente du cadre traditionnellement mis en œuvre avec ce genre de modèle. Son atout est la prédiction du nombre de dimensions du problème d'identification, donnant à son utilisateur son positionnement quant à l'unicité de la solution (ce qui ne sera pas le cas pour une méthode classique de minimisation des moindres carrés par exemple). Cette méthode est basée sur un méta-modèle obtenu après représentation des principales variations des données d'entrée dans une base de dimension réduite (par POD), puis interpolation (par krigeage) des données discrètes en une hypersurface continue. La dimension de cette dernière donne des informations sur l'influence des paramètres et permet donc d'estimer leur identifiabilité. Les données expérimentales peuvent être projetées dans la base de POD puis les paramètres de la loi de comportement considérée sont identifiés par projection sur le méta-modèle.

L'approche est employée dans un premier temps sur le deuxième des trois modèles présentés ci-avant. Les paramètres recherchés peuvent être identifiés en deux temps : ceux de la partie sphérique sur des essais de compression hydrostatique purs, et ceux de la partie déviatorique sur des essais faisant intervenir à la fois des contributions déviatoriques et sphériques (traction, compression confinée). Les variables liées à l'effet Mullins ont été identifiées pour les deux parties de la transformation (isochore/volumique) respectivement sur des essais cyclés de traction et compression confinée (à amplitude de déformation croissante). Les résultats obtenus permettent de reproduire des niveaux de contrainte adéquats. Les enveloppes charges/décharges sont également plutôt bien prises en compte traduisant convenablement l'effet Mullins, y compris sur la compression confinée, objectif impossible à réaliser sans la prise en compte d'une variable d'adoucissement sur la partie volumique de l'énergie libre. Seule la déformation résiduelle n'est évidemment pas reproduite de par le caractère hyperélastique adopté. Cette dernière est un phénomène du premier ordre pour l'application visée par le CEA. Il faut donc identifier les paramètres du modèle élasto-(visco-)plastique complet. Le méta-modèle obtenu présente un souci de dimensionnalité ne permettant pas d'identifier tous les paramètres. Ceci n'est pas synonyme d'un défaut de la formulation de comportement proposée, mais traduit plutôt un manque vraisemblable de la base expérimentale, et particulièrement d'essais donnant plus d'importance à certains paramètres du modèle (ce point sera abordé dans la partie perspectives de la présente conclusion).

La suite présente les premières perspectives de ce travail. Elles sont expérimentales dans un premier temps. En vue d'améliorer l'identification et de la simplifier également, il faudrait idéalement réaliser les campagnes d'essais suivantes :

<figure><img src="image_92.png" /></figure>

essais purement déviatoriques comme du glissement simple, dont le pilotage devra permettre du cyclage et de la relaxation. Les essais présentés dans le deuxième chapitre de ce manuscrit présentent des niveaux de contrainte très faibles en raison du montage d'essais. Il faudrait sûrement améliorer la fixation des éprouvettes pour que leurs bords soient mieux intégrés dans les mors. Après cela, il faudrait vérifier que le chargement est déviatorique par corrélation d'images.
essais purement sphériques permettant encore une fois du cyclage et un contrôle de la vitesse de sollicitation, cela permettrait d'identifier aussi les paramètres de la partie volumique visco-plastique du modèle. Des essais de relaxation seraient un plus non négligeable.

Sur ces deux types de chargement, un pilotage de type relaxation par palier permettrait d'exhiber le comportement à l'équilibre du matériau (i.e. enveloppe après les paliers de relaxation) et de s'assurer de l'identification du comportement élasto-plastique en s'affranchissant de tout artefact visqueux.

Quelques essais préliminaires de relaxation ont été réalisés et ils mettent en évidence l'influence du temps sur la réponse de nos matériaux sous de tels chargements. Pour prendre en compte ces phénomènes dans la modélisation et en première approche, on commence par donner un sens physique aux deux paramètres de viscosité, en augmentant leurs valeurs. Cela permet de modéliser la relaxation à court terme et une certaine influence de la vitesse de sollicitation sur la réponse mécanique.

L'importance de cette dépendance visqueuse a été renforcée par la mise à notre disposition d'une campagne d'essais riche et originale menée par le CEA. Contrairement à nos essais classiques de relaxation caractérisant les temps visco-élastiques caractéristiques courts, celle-ci propose de s'intéresser aux homologues longs termes. Cela est permis grâce à des essais de fluage de plusieurs mois à diverses températures puis de suivi de recouvrance sur des périodes tout aussi longues. La prise en compte de tels phénomènes dans notre modèle nécessite *a minima* l'intégration d'autres temps caractéristiques. La première ouverture donnée en ce sens est la définition d'autres variables internes de viscosité, facilement introduites dans notre formulation (avantage de la formulation eulérienne) par le postulat de nouveaux états intermédiaires. Cette modélisation phénoménologique semble être une piste intéressante, mais ne serait qu'un point de départ étant donnée l'influence des phénomènes présents. Il faudrait sans doute recourir à l'adjonction à la mécanique d'autres physiques telles que la thermique et la chimie. La prise en compte de ces phénomènes amène à revoir aussi totalement l'implémentation éléments finis pour intégrer des degrés de liberté thermiques et chimiques aux noeuds. Il faudrait également avoir une base expérimentale représentative avec des essais de vieillissement sous contrainte et sous différentes températures (à l'image de celle proposée par le CEA).

Sur une perspective à plus court terme, sur le plan expérimental, d'autres investigations pourraient être menées. Notamment la vérification de l'hypothèse forte sur laquelle nous basons notre modèle ($\mathbf{F} = \hat{\mathbf{F}} \cdot \bar{\mathbf{F}}$). Cela permettrait peut-être, à défaut de la valider, d'établir un domaine de validité de notre hypothèse. Si ce domaine existe, il serait alors intéressant de le positionner par rapport à celui considéré pour les applications du CEA. En outre, il semble indispensable de s'intéresser à la microstructure du matériau sous chargement, via par exemple des essais mécaniques in-situ en microtomographie afin de :

- comprendre les mécanismes de déformations des microsphères (type de flambement, réversibilité des déformations...),

- s'assurer de l'intégrité des microsphères après chargement,

- cerner l'interaction entre RTV et microsphères à l'état initial et/ou après vieillissement
thermo-mécanique.

Intrinsèquement intéressantes ces informations pourraient aussi être la base pour bâtir des modèles d'homogénéisation numérique qui pourraient remplacer/compléter une base expérimentale non évidente à construire.

Il y a d'autres perspectives du côté numérique, notamment explorer la possibilité de rendre l'opérateur tangent « plus analytique ». Une voie simple (initialement prévue dans ce travail, mais non menée par manque de temps) serait de construire analytiquement la partie inhérente au comportement hyperélastique (avec Mullins) et une contribution totalement numérique liée aux phénomènes élasto-(visco-)plastiques.

Enfin, toujours à court terme, il serait intéressant de se servir de l'approche modèle réduit pour prédéfinir les essais nécessaires à l'identification de tous les paramètres du modèle (il faudrait obtenir une dimension égale au nombre de paramètres recherchés avec les simulations des essais). La logique serait inversée par rapport à l'approche classique, car c'est le modèle qui guiderait la construction de la base expérimentale nécessaire à sa bonne identification.
# Annexe A

## Grandes déformations et cadre thermodynamique

### I.A Grandes déformations

#### I.A.1 Description du mouvement

##### I.A.1.a Position

Contrairement au cas classique des petites perturbations, où la déformation d'un solide est faible, quand on sollicite un élastomère, comme dans cette thèse, les déformations engendrées sont considérables, il est alors d'usage de faire appel à plusieurs configurations, à plusieurs temps $t$, pour décrire le mouvement.

La fonction $\chi$ décrit le mouvement d'un solide entre sa position dans la configuration initiale $C_0$ à l'instant $t_0$ et sa position à $t$ dans la configuration actuelle $C(t)$. La description du mouvement est dans ce cas, lagrangienne, dépendante du temps et des positions initiales du point matériel.

$$ \mathbf{x} = \chi(\mathbf{X}, t) \qquad (A.1) $$

Cette fonction du mouvement est bijective de $C_0$ dans $C(t)$ et permet ainsi de retrouver la position initiale du point X dans $C_0$ en fonction de sa position dans $C(t)$. Il s'agit d'une description eulérienne du mouvement SEIBERT 2017 où l'on repère les points par leurs coordonnées actuelles dépendant du temps.

$$ \mathbf{X} = \chi^{-1}(\mathbf{x}, t) \qquad (A.2) $$

Comme l'indique la figure A.1, les coordonnées d'un point X sont notées $\mathbf{X}$ dans $C_0$ et $\mathbf{x}$ dans $C(t)$.

La position peut également être indiquée par l'introduction d'un vecteur déplacement $\mathbf{u}$, sous la forme SIDOROFF 1982 :

$$ \mathbf{x}(\mathbf{X}, t) = \mathbf{X} + \mathbf{u}(\mathbf{X}, t) \qquad (A.3) $$
<figure><img src="image_93.png" /><figcaption>FIGURE A.1 – Configurations initiales et actuelles</figcaption></figure>

### I.A.1.b Vitesse et accélération

En variables lagrangiennes de $C_0$, la vitesse du point dépend du temps et de ses coordonnées à $t_0$.

$$ \dot{\mathbf{x}} = \mathbf{v}(t, \mathbf{X}) $$

La description eulérienne fait intervenir pour la vitesse, le temps et les coordonnées actuelles du point lors de son mouvement.

$$ \dot{\mathbf{x}} = \mathbf{v}(t, \mathbf{x}) $$

L'accélération de ce point va alors être la dérivée totale de la vitesse. Les coordonnées étant indépendantes du temps en description lagrangienne, on a :

$$ \mathbf{a} = \ddot{\mathbf{x}} = \frac{d\mathbf{v}}{dt}(\mathbf{X}, t) = \frac{\partial \mathbf{v}}{\partial t} + \frac{\partial \mathbf{v}}{\partial \mathbf{X}} \frac{\partial \mathbf{X}}{\partial t} = \frac{\partial \mathbf{v}}{\partial t} \quad (A.4) $$

La dérivée totale est donc égale à la dérivée partielle dans le cas lagrangien. Dans la description eulérienne cependant, la position dépendant du temps, on a alors :

$$ \begin{aligned} \mathbf{a} = \ddot{\mathbf{x}} &= \frac{d\mathbf{v}}{dt}(\mathbf{x}(t), t) = \frac{\partial \mathbf{v}}{\partial t} + \frac{\partial \mathbf{v}}{\partial \mathbf{x}} \frac{\partial \mathbf{x}}{\partial t} \\ &= \frac{\partial \mathbf{v}}{\partial t}(\mathbf{x}, t) + \operatorname{grad}_x \mathbf{v} \cdot \mathbf{v} \end{aligned} \quad (A.5) $$

Dans cette expression le premier terme est appelé local, tandis que le second désigne une partie convective LEMAÎTRE et al. 2009.

Cette forme de dérivée est aussi appelée dérivée particulière et peut s'étendre à toute fonction $\mathbf{f}(\mathbf{x}, t)$ d'ordre 0 ou 1 dépendant du mouvement et du temps de sorte que :

$$ \dot{\mathbf{f}}(\mathbf{x}, t) = \frac{d\mathbf{f}}{dt}(\mathbf{x}, t) = \frac{\partial \mathbf{f}}{\partial t} + \operatorname{grad}_x \mathbf{f} \cdot \mathbf{v} \quad (A.6) $$
## I.A.2 Transport de vecteurs, éléments de surface et de volume

On introduit un vecteur qui relie deux points infinitésimalement proches. Il s'agit de $d\mathbf{X}$ pour $C_0$ et $d\mathbf{x}$ pour $C(t)$, tels que :

$$d\mathbf{X} = \mathbf{Y} - \mathbf{X} \quad \text{et} \quad d\mathbf{x} = \mathbf{y} - \mathbf{x} \qquad (A.7)$$

Pour décrire la déformation locale entre un point $X$ et un point $Y$, on peut écrire le dévlopement limité de la fonction de mouvement en $Y$ au voisinage de $X$ :

$$\chi(\mathbf{Y}, t) = \chi(\mathbf{X}, t) + \frac{\partial\chi(\mathbf{X}, t)}{\partial\mathbf{X}}(\mathbf{Y} - \mathbf{X}) + o(\mathbf{Y} - \mathbf{X})^2$$

Or d'après (A.1) et (A.7) on peut écrire :

$$d\mathbf{x} = \mathbf{y} - \mathbf{x} = \chi(\mathbf{Y}, t) - \chi(\mathbf{X}, t) \qquad (A.8)$$

L'expression (I.A.2) devient donc :

$$d\mathbf{x} = \frac{\partial\mathbf{x}}{\partial\mathbf{X}}d\mathbf{X} \qquad (A.9)$$

L'équation (A.9) met en évidence le gradient de déformation qui relie un vecteur de la configuration de référence à la configuration actuelle. Ce tenseur est noté $\mathbf{F}$ et a donc pour expression :

$$\mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \frac{\partial \mathbf{x}_i}{\partial \mathbf{X}_j} \mathbf{e}_i \otimes \mathbf{e}_j \qquad (A.10)$$

De même des éléments de surface et de volume peuvent être transportés dans $C(t)$ comme suit :

$$dv = \det \mathbf{F} dV \qquad (A.11)$$

$$\mathbf{n} \cdot \mathbf{dS} = (\det \mathbf{F}) \mathbf{F}^{-T} \cdot \mathbf{N} dS \qquad (A.12)$$

Dans ces expressions, le déterminant du gradient de la transformation, aussi dénommé le Jacobi, est noté $J$. La relation (A.11) met en évidence que lorsque la transformation est incompressible, i.e. le volume du point matériel reste inchangé, alors :

$$J = \det \mathbf{F} = 1 \qquad (A.13)$$

## I.A.3 Mesure des déformations

### I.A.3.a Décomposition polaire

La seule connaissance de $\mathbf{F}$ n'est pas adaptée pour décrire la déformation d'un corps puisque ce tenseur contient à la fois des informations sur la translation mais aussi sur sa rotation. Pour pallier ce problème, il faut faire appel à la décomposition polaire de $\mathbf{F}$. Elle repose sur la multiplication d'un tenseur décrivant une translation et d'un autre décrivant une rotation.

$$\mathbf{F} = \mathbf{R} \cdot \mathbf{U} = \mathbf{V} \cdot \mathbf{R} \qquad (A.14)$$
Dans l'équation (A.14), **R** est le tenseur des rotations, orthogonal. **U** et **V** sont respectivement les tenseurs de déformation pure droit et gauche. Ils sont symétriques et définis positifs.

<figure><img src="image_94.png" /><figcaption>FIGURE A.2 – Décomposition translation et rotation</figcaption></figure>

Ces deux tenseurs peuvent être transportés l'un et l'autre entre les configurations initiale et actuelle par les relations en (A.15).

$$ 
\begin{aligned}
\mathbf{V} &= \mathbf{R} \cdot \mathbf{U} \cdot \mathbf{R}^{\mathrm{T}} \\
\mathbf{U} &= \mathbf{R}^{\mathrm{T}} \cdot \mathbf{V} \cdot \mathbf{R}
\end{aligned}
\qquad (\mathrm{A.15})
$$

Dans le cas du déplacement d'un corps rigide, non déformable, ces tenseurs valent l'identité 1. À partir de là, on peut définir des mesures généralisées des déformations :

— en configuration lagrangienne

$$ 
\left\{ 
\begin{array}{l}
\frac{1}{n}(\mathbf{U}^n - \mathbf{1}) \text{ si } n \neq 0 \\
\log \mathbf{U} \text{ sinon}
\end{array}
\right. 
\qquad (\mathrm{A.16})
$$

— en configuration eulérienne

$$ 
\left\{ 
\begin{array}{l}
\frac{1}{n}(\mathbf{V}^n - \mathbf{1}) \text{ si } n \neq 0 \\
\log \mathbf{V} \text{ sinon}
\end{array}
\right. 
\qquad (\mathrm{A.17})
$$

### I.A.3.b Définition directe

Les déformations peuvent également être définies à partir du produit scalaire de vecteurs des deux configurations. Si on considère $d\mathbf{X}$ et $d\mathbf{Y}$ de $C_0$ et $dx$ et $dy$ de $C(t)$, on peut écrire :

$$ 
\begin{aligned}
d\mathbf{x} \cdot d\mathbf{y} &= (\mathbf{F} \cdot d\mathbf{X}) \cdot (\mathbf{F} \cdot d\mathbf{Y}) \\
&= d\mathbf{X} \cdot \mathbf{F}^T \cdot \mathbf{F} \cdot d\mathbf{Y} \\
&= d\mathbf{X} \cdot \mathbf{C} \cdot d\mathbf{Y}
\end{aligned}
\qquad (\mathrm{A.18})
$$
Par identification on obtient le tenseur des dilatations de CAUCHY-GREEN droit :

$$ C = F^T \cdot F \qquad (A.19) $$

En combinant les relations (A.14) et (A.19), on a grâce à la symétrie de $U$ :

$$ C = F^T \cdot F = (R \cdot U)^T \cdot (R \cdot U) = U^T \cdot U = U^2 \qquad (A.20) $$

De même dans la configuration actuelle, il vient :

$$ \begin{aligned} d\mathbf{X} \cdot d\mathbf{Y} &= (\mathbf{F}^{-1} \cdot d\mathbf{x}) \cdot (\mathbf{F}^{-1} \cdot d\mathbf{y}) \\ &=d\mathbf{x} \cdot \mathbf{F}^{-T} \cdot \mathbf{F}^{-1} \cdot d\mathbf{y} \qquad (A.21) \\ &=d\mathbf{x} \cdot \mathbf{B}^{-1} \cdot d\mathbf{y} \end{aligned} $$

Avec $\mathbf{B}$ le tenseur de dilatations de CAUCHY-GREEN gauche :

$$ \mathbf{B} = \mathbf{F} \cdot \mathbf{F}^T \qquad (A.22) $$

La combinaison de (A.14) et (A.22) livre ici :

$$ \mathbf{B} = \mathbf{F} \cdot \mathbf{F}^T = (\mathbf{V} \cdot \mathbf{R}) \cdot (\mathbf{V} \cdot \mathbf{R})^T = \mathbf{V} \cdot \mathbf{V}^T = \mathbf{V}^2 \qquad (A.23) $$

Dans un état non déformé, le gradient de la transformation vaut l'identité, c'est donc aussi le cas des tenseurs de dilatation $\mathbf{B}$ et $\mathbf{C}$. On peut alors utiliser une autre définition de déformation, qui elle sera plus conventionnelle, puisqu'à un état non déformé, la valeur de la déformation sera nulle. Pour cela on considère la déformation comme une différence de produits scalaires des vecteurs des configurations $C_0$ et $C(t)$ :

$$ \begin{aligned} d\mathbf{x} \cdot d\mathbf{y} - d\mathbf{X} \cdot d\mathbf{Y} &= (\mathbf{F} \cdot d\mathbf{X}) \cdot (\mathbf{F} \cdot d\mathbf{Y}) - d\mathbf{X} \cdot d\mathbf{Y} \\ &=d\mathbf{X} \cdot \mathbf{F}^T \cdot \mathbf{F} \cdot d\mathbf{Y} - d\mathbf{X} \cdot d\mathbf{Y} \qquad (A.24) \\ &=d\mathbf{X} \cdot (\mathbf{C} - \mathbf{I}) \cdot d\mathbf{Y} \\ &=d\mathbf{X} \cdot 2\mathbf{E} \cdot d\mathbf{Y} \end{aligned} $$

Il apparaît le tenseur des déformations dans la configuration lagrangienne, dit de GREEN-LAGRANGE :

$$ \mathbf{E} = \frac{1}{2} (\mathbf{C} - \mathbf{I}) \qquad (A.25) $$

En reproduisant l'opération pour la configuration actuelle, on a :

$$ \begin{aligned} d\mathbf{x} \cdot d\mathbf{y} - d\mathbf{X} \cdot d\mathbf{Y} &=d\mathbf{x} \cdot d\mathbf{y} - (\mathbf{F}^{-1} \cdot d\mathbf{x}) \cdot (\mathbf{F}^{-1} \cdot d\mathbf{y}) \\ &=d\mathbf{x} \cdot d\mathbf{y} - d\mathbf{x} \cdot \mathbf{F}^{-T} \cdot \mathbf{F}^{-1} \cdot d\mathbf{y} \qquad (A.26) \\ &=d\mathbf{x} \cdot (\mathbf{I} - \mathbf{B}^{-1}) \cdot d\mathbf{y} \\ &=d\mathbf{x} \cdot 2\mathbf{A} \cdot d\mathbf{y} \end{aligned} $$
On obtient le tenseur des déformations dans la configuration eulérienne, dit de EULER-ALMANSI :

$$A = \frac{1}{2} (I - B^{-1}) \qquad (A.27)$$

### I.A.4 Vitesse de déformation

La variation de vitesse d'un vecteur s'écrit :

$$ \begin{aligned} (dx)^{\bullet} &= \frac{d}{dt} (F \cdot dX) \\ &= \dot{F} \cdot dX \end{aligned} \qquad (A.28) $$

Le vecteur $dX$ de $C_0$ étant indépendant du temps, il ne restait plus qu'à dériver le gradient de déformation. Or $\dot{F}$ peut s'écrire :

$$ \dot{F} = \frac{dF}{dt} = \frac{d}{dt} \left( \frac{\partial x}{\partial X} \right) = \frac{\partial}{\partial X} \left( \frac{\partial x}{\partial t} \right) = \mathrm{grad}_X v \qquad (A.29) $$

C'est le gradient de vitesse dans la configuration initiale $C_0$.

En exprimant $dX$ dans la configuration actuelle, la relation (A.28) devient :

$$ \begin{aligned} (dx)^{\bullet} &= \dot{F} \cdot F^{-1} \cdot dx \\ &= L \cdot dx \end{aligned} \qquad (A.30) $$

$L$ peut être écrit à la manière de l'équation (A.29), on obtient :

$$ L = \dot{F} \cdot F^{-1} = \frac{\partial v}{\partial X} \frac{\partial X}{\partial x} = \frac{\partial v}{\partial x} = \mathrm{grad}_x v \qquad (A.31) $$

$L$ est donc le gradient de vitesse dans la configuration actuelle $C(t)$, il s'écrit également :

$$ L = D + W \qquad (A.32) $$

$D$ est le symétrique de $L$ représentant la vitesse de déformation et $W$ est la partie anti-symétrique décrivant le rotationnel du champ de vitesse. Ces tenseurs s'écrivent aussi :

$$ D = \frac{1}{2} (L + L^T) \qquad W = \frac{1}{2} (L - L^T) \qquad (A.33) $$

Une interprétation de la vitesse de déformation peut également être vue comme la variation de vitesse du scalaire de deux vecteurs. Dans la configuration initiale, on a :

$$ \begin{aligned} (dx \cdot dx)^{\bullet} &= d\dot{x} \cdot dx + dx \cdot d\dot{x} = (\dot{F} \cdot F^{-1} \cdot dx) \cdot dx + dx \cdot (\dot{F} \cdot F^{-1} \cdot dx) \\ &= (\dot{F} \cdot dX) \cdot F \cdot dX + F \cdot dX \cdot (\dot{F} \cdot dX) \\ &= dX \cdot (\dot{F}^T \cdot F + F^T \cdot \dot{F}) dX \\ &= dX \cdot 2\dot{E} \cdot dX \end{aligned} \qquad (A.34) $$
Pour la configuration eulérienne :

$$ 
\begin{align} 
(dx \cdot dx)^* &= \dot{dx} \cdot dx + dx \cdot \dot{dx} = (\mathbf{L} \cdot dx) \cdot dx + dx \cdot (\mathbf{L} \cdot dx) \nonumber \\
&= dx \cdot (\mathbf{L} + \mathbf{L}^T) dx \tag{A.35} \\
&= dx \cdot 2\mathbf{D} \cdot dx \nonumber
\end{align} 
$$

Or comme $\dot{\mathbf{F}} = \mathbf{L} \cdot \mathbf{F}$ et $\dot{\mathbf{F}}^T = \mathbf{F}^T \cdot \mathbf{L}^T$, la ligne 3 de l'équation (A.34), révélant l'expression de $\dot{\mathbf{E}}$ peut se réécrire :

$$ 
\begin{align} 
2\dot{\mathbf{E}} &= \dot{\mathbf{F}}^T \cdot \mathbf{F} + \mathbf{F}^T \cdot \dot{\mathbf{F}} \nonumber \\
&= \mathbf{F}^T \cdot \mathbf{L}^T \cdot \mathbf{F} + \mathbf{F}^T \cdot \mathbf{L} \cdot \mathbf{F} \nonumber \\
&= \mathbf{F}^T \cdot (\mathbf{L} + \mathbf{L}^T) \cdot \mathbf{F} \tag{A.36} \\
&= \mathbf{F}^T \cdot 2\mathbf{D} \cdot \mathbf{F} \nonumber
\end{align} 
$$

Cela vérifie bien le transport des deux tenseurs d'une configuration à l'autre.

## I.A.5 Description des contraintes

Il y a plusieurs possibilités afin de décrire les efforts intérieurs subis par un corps matériel, selon la configuration considérée. Les grandeurs utilisées sont présentées dans la figure A.3.

<figure><img src="image_95.png" /><figcaption>FIGURE A.3 – Configurations initiales et actuelles, état de contrainte</figcaption></figure>

### I.A.5.a Description eulérienne

En se plaçant dans la configuration actuelle $C(t)$, eulérienne, la contrainte peut être vue comme le rapport des efforts intérieurs df sur un élément de surface déformé ds de normale n
et par unité de surface déformée. On écrit alors la relation issue du théorème des contraintes de CAUCHY HOLZAPFEL 2000,

$$df = \sigma \cdot nds \qquad (A.37)$$

$\sigma$ est le tenseur des contraintes de CAUCHY, dont la définition est également utilisée dans le cadre des petites déformations. Il est symétrique donc $\sigma = \sigma^T$.

### I.A.5.b Description mixte

Le tenseur de Cauchy peut être transporté dans la configuration non déformée, afin d'obtenir une description de la contrainte équivalente à un effort de cohésion par élément de surface déformée mais par unité de surface non déformée. Cela définit donc un nouveau tenseur des contraintes. On emploie, pour cela, la formule de Nanson HOLZAPFEL 2000 (A.38), équivalente à la formule (A.12) transportant les éléments de surface.

$$nds = JF^{-T} \cdot N \cdot dS \qquad (A.38)$$

(A.37) peut alors être modifiée :

$$df = \sigma \cdot nds = J\sigma \cdot F^{-T} \cdot NdS \qquad (A.39)$$

Il en ressort :

$$\pi = J\sigma \cdot F^{-T} \qquad (A.40)$$

Avec $\pi$ le tenseur des contraintes de PIOLA-KIRCHHOFF I, étant ainsi défini comme une application de $C_0$ dans $C(t)$ SIDOROFF 1982. C'est donc une description mixte qui n'est ni lagrangienne, ni eulérienne. Son utilisation est très utile dans la définition des conditions aux limites, puisque les efforts appliqués sont exprimés par unité de surface non déformée, beaucoup plus facile à obtenir expérimentalement SIDOROFF 1982. Il n'est pas symétrique.

### I.A.5.c Description lagrangienne

Une autre possibilité d'exprimer la contrainte consiste à transporter l'effort intérieur de cohésion df dans la configuration de référence.

$$ \begin{aligned} df_0 &= F^{-1} \cdot df \\ &= JF^{-1} \cdot \sigma \cdot F^{-T} \cdot NdS \qquad (A.41) \\ &= F^{-1} \cdot \pi \cdot NdS \\ &= S \cdot NdS \end{aligned} $$

Le tenseur S de l'équation (A.41) est le tenseur des contraintes de PIOLA-KIRCHHOFF II. Il est totalement lagrangien mais sans signification physique. Comme il est défini entièrement dans cette configuration initiale et qu'il est symétrique, il permet une implémentation numérique simplifiée SEIBERT 2017.
# I.B Cadre thermodynamique

Les notions de cinématique en grandes déformations introduites à la section précédente décrivent le mouvement et les déformations d'un corps matériel. Les relations sont valables indépendamment des propriétés matérielles et des efforts opérant sur le corps. Dans la configuration actuelle, l'état du corps matériel $\omega$ dépend de sa masse, sa quantité de mouvement, son moment cinétique, dans une description mécanique. Thermodynamiquement parlant, sa température et son énergie entrent également en compte. Ses interactions avec l'environnement sont, par exemple, la puissance mécanique, la chaleur et l'entropie apportées par l'extérieur au corps, au cours de son mouvement. Cette section va présenter les équilibres physiques faisant le lien entre les sources du mouvement et l'influence de l'environnement.

## I.B.1 Équations d'équilibre de la mécanique

### I.B.1.a Équilibre de masse

L'équilibre de masse provient de la conservation de la masse du corps pendant son mouvement. Mathématiquement, cet axiome s'écrit HAUPT 2002 :

— en configuration lagrangienne :

$$m = \int_{\Omega} dm = \int_{\Omega} \rho_{0}(\mathbf{X}, t) \, dV = c_{ste} \qquad (A.42)$$

— en configuration eulerienne :

$$m = \int_{\omega} dm = \int_{\omega} \rho(\mathbf{x}, t) \, dv = c_{ste} \qquad (A.43)$$

La dérivée temporelle de cette intégrale s'écrit donc :

— en configuration lagrangienne :

$$\frac{d}{dt} \int_{\Omega} \rho_{0}(\mathbf{X}, t) \, dV = 0 \qquad (A.44)$$

Il s'agit du bilan global de masse en configuration de Lagrange. L'intégrale portant sur un volume constant, les opérateurs dérivé et intégrale peuvent être intervertis :

$$\int_{\Omega} \frac{\partial}{\partial t} \rho_{0}(\mathbf{X}, t) \, dV = 0 \qquad (A.45)$$

En faisant l'hypothèse que la fonction intégrée est continue par rapport à $\mathbf{X}$, l'intégrale (A.45) est valable sur $\Omega$ mais également sur des sous-ensembles arbitraires de $\Omega$ HAUPT 2002. Localement, il vient alors :

$$\frac{\partial}{\partial t} \rho_{0}(\mathbf{X}, t) = 0 \qquad (A.46)$$
— en configuration eulérienne :

$$ \frac{d}{dt} \int_{\omega} \rho(\mathbf{x}, t) \, dv = 0 \qquad (A.47) $$

Il s'agit de la forme globale du bilan de masse. L'obtention d'une forme locale de cet équilibre peut également être intéressante.

Le théorème de transport de Reynolds s'écrit SCHEFFER 2016,

$$ \frac{d}{dt} \int_{\omega} \rho(\mathbf{x}, t) \, dv = \int_{\omega} (\dot{\rho}(\mathbf{x}, t) + \rho(\mathbf{x}, t) \operatorname{div}_{x} \mathbf{v}) \, dv \qquad (A.48) $$

Localement, il vient donc :

$$ \dot{\rho}(\mathbf{x}, t) + \rho(\mathbf{x}, t) \operatorname{div}_{x} \mathbf{v} = 0 \qquad (A.49) $$

Le premier terme de la relation (A.47) peut aussi être écrit de la sorte :

$$ \begin{aligned} \frac{d}{dt} \int_{\omega} \rho(\mathbf{x}, t) \, dv &= \frac{d}{dt} \int_{\Omega} \rho(\mathbf{X}, t) \det \mathbf{F} \, dV \\ &= \int_{\Omega} \frac{\partial}{\partial t} [\rho(\mathbf{X}, t) \det \mathbf{F}] \, dV \end{aligned} \qquad (A.50) $$

Pour obtenir une autre forme du bilan local de masse, les équations (A.45) et (A.50) sont combinées :

$$ \int_{\Omega} \frac{\partial}{\partial t} \rho_0(\mathbf{X}, t) \, dV = \int_{\Omega} \frac{\partial}{\partial t} [\rho(\mathbf{X}, t) \det \mathbf{F}] \, dV \qquad (A.51) $$

En faisant l'hypothèse que cette intégrale est aussi valable pour d'autres sous-ensembles de $\Omega$, l'expression devient :

$$ \rho(\mathbf{x}) = \rho_0(\mathbf{X}) (\det \mathbf{F})^{-1} \qquad (A.52) $$

Il s'agit de la forme la plus commune permettant le transport de la masse volumique entre $C_0$ et $C(t)$.

### I.B.1.b Équilibre de quantité de mouvement

Pour la quantité de mouvement, l'axiome veut que cette grandeur varie en fonction des forces appliquées au corps. Si on note **I** la quantité de mouvement et **f** les forces agissant sur le corps, on a $\dot{\mathbf{I}} = \mathbf{f}$ BERGSTRÖM 2015.

**I** s'écrit :

$$ \mathbf{I} = \int_{\omega} \dot{\mathbf{x}} \, dm = \int_{\omega} \rho \dot{\mathbf{x}}(\mathbf{x}, t) \, dv \qquad (A.53) $$

Dans **f** sont prises en compte les forces agissant sur le contour $d\omega$ de $\omega$ et les forces exercées sur le volume comme la pesanteur ($\rho\mathbf{b}$ est la densité massique de force).

$$ \mathbf{f} = \int_{\partial \omega} \mathbf{t} \, ds + \int_{\omega} \rho \mathbf{b} \, dv \qquad (A.54) $$
La forme globale de l'équilibre de quantité de mouvement s'écrit ((A.37) est utilisée dans (A.54)) :

$$ \frac{d}{dt} \int_{\omega} \rho \dot{\mathbf{x}}(\mathbf{x}, t) \, dv = \int_{\partial\omega} \boldsymbol{\sigma} \cdot \mathbf{n} \, ds + \int_{\omega} \rho \mathbf{b} \, dv \quad (A.55) $$

Le théorème de transport de Reynolds, avec le changement de variable $\Psi = \rho \mathbf{v}$, devient :

$$ \begin{aligned} \frac{d}{dt} \int_{\omega} \Psi(\mathbf{x}, t) \, dv &= \int_{\omega} (\dot{\Psi}(\mathbf{x}, t) + \Psi(\mathbf{x}, t) \operatorname{div}_{x} \mathbf{v}) \, dv \\ &= \int_{\omega} (\rho \dot{\mathbf{v}} + \dot{\rho} \mathbf{v} + \rho \mathbf{v} \operatorname{div}_{x} \mathbf{v}) \, dv \quad (A.56) \\ &= \int_{\omega} \rho \dot{\mathbf{v}} \, dv \end{aligned} $$

Le théorème de la divergence impose :

$$ \int_{\partial\omega} \boldsymbol{\sigma} \cdot \mathbf{n} \, ds = \int_{\omega} \operatorname{div}_{x} \boldsymbol{\sigma} \, dv \quad (A.57) $$

Après utilisation du théorème de transport de Reynolds et du bilan de masse (A.49), sur le membre de gauche de l'équation (A.55), puis de l'équation (A.57), la forme globale devient :

$$ \int_{\omega} \rho \dot{\mathbf{v}}(\mathbf{x}, t) \, dv = \int_{\omega} \operatorname{div}_{x} \boldsymbol{\sigma} \, dv + \int_{\omega} \rho \mathbf{b} \, dv \quad (A.58) $$

La forme locale de l'équilibre de quantité de mouvement dans la configuration actuelle $C(t)$, est exprimée par :

$$ \rho \dot{\mathbf{v}}(\mathbf{x}, t) = \operatorname{div}_{x} \boldsymbol{\sigma} + \rho \mathbf{b} \quad (A.59) $$

Si on reprend l'équation (A.55) et qu'on transporte les intégrales dans la configuration initiale, il vient :

$$ \frac{d}{dt} \int_{\Omega} \rho \dot{\mathbf{x}}(\mathbf{X}, t) \, \mathrm{det} \mathbf{F} \, dV = \int_{\partial \Omega} \boldsymbol{\sigma} \cdot \left( (\mathrm{det} \mathbf{F}) \mathbf{F}^{-T} \cdot \mathbf{N} \, dS \right) + \int_{\Omega} \rho \mathbf{b} \, \mathrm{det} \mathbf{F} \, dV \quad (A.60) $$

Le transport fait apparaître l'expression du premier tenseur de PIOLA-KIRCHHOFF. On peut également faire usage du transport de masse volumique (A.52). On obtient donc :

$$ \frac{d}{dt} \int_{\Omega} \rho_0 \dot{\mathbf{x}}(\mathbf{X}, t) \, dV = \int_{\partial\Omega} \boldsymbol{\pi} \cdot \mathbf{N} \, dS + \int_{\Omega} \rho_0 \mathbf{b} \, dV \quad (A.61) $$

En suivant la même méthode que dans la configuration $C(t)$, avec l'emploi des théorèmes de Reynolds et de la divergence on obtient le bilan local de la quantité de mouvement dans la configuration de référence $C_0$ :

$$ \rho_0 \dot{\mathbf{v}}(\mathbf{X}, t) = \operatorname{div}_{x} \boldsymbol{\pi} + \rho_0 \mathbf{b} \quad (A.62) $$
I.B.2 Premier Principe de la thermodynamique

Le premier principe de la thermodynamique stipule que l'énergie d'un corps matériel varie en fonction des puissances mécanique et thermique qu'il reçoit. Mathématiquement cela équivaut à :

$$
\dot{\mathcal{K}} + \dot{\mathcal{E}} = \mathcal{P}_{\text{ext}} + \mathcal{Q} \qquad (\text{A.63})
$$

$\mathcal{K}$ est l'énergie cinétique du corps et $\mathcal{E}$ son énergie interne. On note $e$ son énergie interne spécifique.

$$
\dot{\mathcal{K}} + \dot{\mathcal{E}} = \frac{d}{dt} \int_{\omega} \left( \frac{1}{2} \rho \dot{\mathbf{x}} \cdot \dot{\mathbf{x}} + \rho e \right) \mathrm{d}v \quad (\mathrm{A}.64)
$$

La puissance $P_{ext}$ dans la relation (A.63) se décompose en une puissance exercée à la surface du corps et dans son volume par effet de la gravité,

$$
\mathcal{P}_{\text{ext}} = \int_{\partial\omega} \dot{\mathbf{x}} \cdot \mathbf{t} \, \mathrm{d}s + \int_{\omega} \dot{\mathbf{x}} \cdot \rho \, \mathbf{b} \, \mathrm{d}v \qquad (\text{A.65})
$$

et $\mathcal{Q}$ est un flux de chaleur $\mathbf{q}$ (eulérien) à travers la frontière $\partial\omega$ de $\omega$ et un rayonnement $r$
libéré à l'intérieur du corps.

$$
\mathcal{Q} = - \int_{\partial\omega} \mathbf{q} \cdot \mathbf{n} \, ds + \int_{\omega} \rho \, r \, dv \qquad (\text{A.66})
$$

Le bilan global d'énergie dans la configuration actuelle a alors pour expression :

$$
\frac{d}{dt} \int_{\omega} \left( \rho e + \frac{1}{2} \rho \dot{\mathbf{x}} \cdot \dot{\mathbf{x}} \right) \mathrm{d}v = \int_{\partial\omega} \left[ \left( \boldsymbol{\sigma}^T \dot{\mathbf{x}} \right) \cdot \mathbf{n} - \mathbf{q} \cdot \mathbf{n} \right] \mathrm{d}s + \int_{\omega} \left( \rho \dot{\mathbf{x}} \cdot \mathbf{b} + \rho r \right) \mathrm{d}v \quad (\text{A.67})
$$

Pour le terme de gauche, on utilise le théorème de transport de Reynolds introduit en (A.48)
et pour l’intégrale sur le contour, le théorème de la divergence :

$$
\int_{\omega} (\rho \dot{e} + \rho \ddot{\mathbf{x}} \cdot \dot{\mathbf{x}}) \, dv = \int_{\omega} \left( \mathrm{div}_x (\boldsymbol{\sigma}^T \dot{\mathbf{x}}) - \mathrm{div}_x \mathbf{q} + \rho \dot{\mathbf{x}} \cdot \mathbf{b} + \rho r \right) \, dv \quad (\text{A.68})
$$

Le terme en divergence, $\mathrm{div}_x(\boldsymbol{\sigma}^T \dot{\mathbf{x}})$ peut être transformé par la relation :

$$
\mathrm{div} (\boldsymbol{\sigma}^T \dot{\mathbf{x}}) = (\mathrm{div} \boldsymbol{\sigma}) \cdot \dot{\mathbf{x}} + \boldsymbol{\sigma} : \mathrm{grad} \dot{\mathbf{x}} = (\mathrm{div} \boldsymbol{\sigma}) \cdot \dot{\mathbf{x}} + \boldsymbol{\sigma} : \mathbf{D} \quad (\text{A.69})
$$

L'équation (A.68) s'écrit donc :

$$
\int_{\omega} [(\rho \ddot{\mathbf{x}} - \mathrm{div}_x \boldsymbol{\sigma} - \rho \mathbf{b}) \cdot \dot{\mathbf{x}} + \rho \dot{e} + \mathrm{div}_x \mathbf{q} - \rho r - \boldsymbol{\sigma} : \mathbf{D}] \, \mathrm{d}v = 0 \quad (\text{A.70})
$$

On reconnaît l'équilibre local de quantité de mouvement (A.59) multiplié par la vitesse, ce terme est donc nul et l'équilibre d'énergie en configuration eulérienne s'écrit :

$$
\int_{\omega} (\rho \dot{e} + \mathrm{div}_x \mathbf{q} - \rho r - \boldsymbol{\sigma} : \mathbf{D}) \, \mathrm{d}v = 0 \quad (\text{A.71})
$$
Pour transporter cet équilibre dans les configurations lagrangienne et mixte, on a besoin du transport du flux de chaleur, il est issu de la relation,

$$ \mathbf{Q} \cdot \mathrm{NdS} = \mathbf{q} \cdot \mathrm{nds} = \mathbf{q} \cdot \left[ (\det \mathbf{F}) \mathbf{F}^{-T} \mathrm{NdS} \right] = \left[ (\det \mathbf{F}) \mathbf{F}^{-1} \mathbf{q} \right] \cdot \mathrm{NdS} \quad (\text{A.72}) $$

avec **q** le flux de chaleur de CAUCHY dans $C(t)$ et **Q** le flux de chaleur de PIOLA-KIRCHHOFF dans $C_0$ HAUPT 2002.

Ensuite il faut modifier le terme représentant la puissance des efforts intérieurs respectivement dans les configurations lagrangienne puis mixte HOLZAPFEL 2000 :

$$ \begin{align} \int_{\omega} \boldsymbol{\sigma} : \mathbf{D} \, \mathrm{d}V &= \int_{\Omega} J \boldsymbol{\sigma} : \mathbf{D} \, \mathrm{d}V = \int_{\Omega} J \boldsymbol{\sigma} : \mathbf{F}^{-T} \dot{\mathbf{E}} \mathbf{F}^{-1} \mathrm{d}V \nonumber \\ &= \int_{\Omega} J \mathbf{F}^{-1} \boldsymbol{\sigma} \mathbf{F}^{-T} : \dot{\mathbf{E}} \, \mathrm{d}V \tag{A.73a} \\ &= \int_{\Omega} \mathbf{S} : \dot{\mathbf{E}} \, \mathrm{d}V \nonumber \end{align} $$

$$ \begin{align} \int_{\omega} \boldsymbol{\sigma} : \mathbf{D} \, \mathrm{d}V &= \int_{\Omega} J \boldsymbol{\sigma} : \mathbf{D} \, \mathrm{d}V = \int_{\Omega} J \boldsymbol{\sigma} : \dot{\mathbf{F}} \mathbf{F}^{-1} \mathrm{d}V \nonumber \\ &= \int_{\Omega} J \boldsymbol{\sigma} \mathbf{F}^{-T} : \dot{\mathbf{F}} \, \mathrm{d}V \tag{A.73b} \\ &= \int_{\Omega} \boldsymbol{\pi} : \dot{\mathbf{F}} \, \mathrm{d}V \nonumber \end{align} $$

En utilisant les équations (A.72), (A.73) et le bilan de masse (A.49) dans l'équilibre (A.71), les équilibres d'énergie peuvent être écrits pour les configurations mixte et lagrangienne.

Les expressions finales des équilibres locaux de conservation de l'énergie sont :

— en configuration eulérienne

$$ \rho \dot{\mathbf{e}} = \boldsymbol{\sigma} : \mathbf{D} - \operatorname{div}_x \mathbf{q} + \rho r \quad (\text{A.74}) $$

— en configuration mixte

$$ \rho_0 \dot{\mathbf{e}} = \boldsymbol{\pi} : \dot{\mathbf{F}} - \operatorname{div}_x \mathbf{Q} + \rho_0 r \quad (\text{A.75}) $$

— en configuration lagrangienne

$$ \rho_0 \dot{\mathbf{e}} = \mathbf{S} : \dot{\mathbf{E}} - \operatorname{div}_x \mathbf{Q} + \rho_0 r \quad (\text{A.76}) $$
I.B.3 Second principe de la thermodynamique

D'après la deuxième loi fondamentale de la thermodynamique, un processus physique doit toujours être accompagné d'une augmentation d'entropie, signifiant que la production d'entropie $S_p$ doit donc toujours être positive.

Mathématiquement, l'équilibre d'entropie décrit l'évolution temporelle de l'entropie S en fonction des apports surfacique ou volumique et de la production d'entropie BERGSTRÖM 2015.

$$ \frac{d}{dt} \int_{\omega} \rho S \, dv = \int_{\partial \omega} -\frac{\mathbf{q}}{T} \cdot \mathbf{n} \, ds + \int_{\omega} \frac{\rho r}{T} \, dv + \int_{\omega} S_p \, dv \quad (A.77) $$

Avec le théorème de transport de Reynolds, le premier terme peut être dérivé. Le théorème de la divergence est employé sur l'intégrale en surface. On obtient, la formulation locale de l'inéquation dans la configuration actuelle $C(t)$ :

$$ \rho \dot{S} + \operatorname{div}_x \left( \frac{\mathbf{q}}{T} \right) - \frac{\rho r}{T} = S_p \geq 0 \quad (A.78) $$

De manière analogue dans $C_0$, il vient :

$$ \rho_0 \dot{S} + \operatorname{div}_X \left( \frac{\mathbf{Q}}{T} \right) - \frac{\rho_0 r}{T} = S_p \geq 0 \quad (A.79) $$

Le terme en divergence peut être transformé par :

$$ \operatorname{div}_x \left( \frac{\mathbf{q}}{T} \right) = \frac{1}{T} \operatorname{div}_x \mathbf{q} - \frac{\mathbf{q}}{T^2} \cdot \operatorname{grad}_x (T) \quad (A.80) $$

En introduisant les bilans d'énergie locaux des différentes configurations resp. (A.74), (A.75) et (A.76) dans les équations précédentes (A.78) et (A.79), on obtient :

— en configuration eulérienne

$$ \rho (\dot{e} + T \dot{S}) + \frac{\mathbf{q}}{T} \cdot \operatorname{grad}_x T - \boldsymbol{\sigma} : \mathbf{D} \geq 0 \quad (A.81) $$

— en configuration mixte

$$ \rho_0 (\dot{e} + T \dot{S}) + \frac{\mathbf{Q}}{T} \cdot \operatorname{grad}_X T - \boldsymbol{\pi} : \dot{\mathbf{F}} \geq 0 \quad (A.82) $$

— en configuration lagrangienne

$$ \rho_0 (\dot{e} + T \dot{S}) + \frac{\mathbf{Q}}{T} \cdot \operatorname{grad}_X T - \mathbf{S} : \dot{\mathbf{E}} \geq 0 \quad (A.83) $$
# Annexe B

## Détails des calculs : modèle hyperélastique
compressible

### II.A Développement de la loi de comportement

L'énergie libre dans ce cas hyperélastique compressible dépend de $J$ et $\bar{\mathbf{B}}$. Sa dérivée s'exprime alors par :

$$ \dot{\psi}(J, \bar{\mathbf{B}}) = \frac{\partial \hat{\psi}}{\partial J} \dot{J} + \frac{\partial \bar{\psi}}{\partial \bar{\mathbf{B}}} : \dot{\bar{\mathbf{B}}} \quad (B.1) $$

On commence par le terme en $\dot{J}$ :

$$ \dot{J} = \frac{1}{2} \left[ (\det \mathbf{B})^{-\frac{1}{2}} \right] \det \bar{\mathbf{B}} \quad (B.2) $$

D'après (ITSKOV 2015), $\det \bar{\mathbf{B}}$ s'écrit :

$$ \det \bar{\mathbf{B}} = \frac{\partial \det \mathbf{B}}{\partial \bar{\mathbf{B}}} : \dot{\mathbf{B}} = \det \mathbf{B} \mathbf{B}^{-T} : \dot{\mathbf{B}} = \det \mathbf{B} \mathbf{B}^{-1} : \dot{\mathbf{B}} $$

Donc,

$$ \dot{J} = \frac{1}{2} \left[ (\det \mathbf{B})^{-\frac{1}{2}} \right] \det \mathbf{B} \mathbf{B}^{-1} : \dot{\mathbf{B}} = \frac{1}{2} (\det \mathbf{B})^{\frac{1}{2}} \mathbf{B}^{-1} : \dot{\mathbf{B}} = \frac{1}{2} J \mathbf{B}^{-1} : \dot{\mathbf{B}} \quad (B.3) $$

On peut aussi écrire $\dot{J}$ en fonction de $\mathbf{F}$ :

$$
\begin{align*}
\dot{J} &= \frac{1}{2} J \mathbf{B}^{-1} : \dot{\mathbf{B}} = \frac{1}{2} J (\mathbf{F}^{-T} \cdot \mathbf{F}^{-1}) : (\dot{\mathbf{F}} \cdot \mathbf{F}^T + \mathbf{F} \cdot \dot{\mathbf{F}}^T) \\
&= \frac{1}{2} J (\mathbf{F}^{-T} \cdot \mathbf{F}^{-1}) : (\dot{\mathbf{F}} \cdot \mathbf{F}^T) + \frac{1}{2} J (\mathbf{F}^{-T} \cdot \mathbf{F}^{-1}) : (\mathbf{F} \cdot \dot{\mathbf{F}}^T) \\
&= \frac{1}{2} J (\mathbf{F}^{-T} \cdot \mathbf{F}^{-1} \cdot \mathbf{F}) : \dot{\mathbf{F}} + \frac{1}{2} J (\mathbf{F}^{-T} \cdot \mathbf{F}^T \cdot \mathbf{F}^{-1}) : \dot{\mathbf{F}}^T \\
&= \frac{1}{2} J \mathbf{F}^{-T} : \dot{\mathbf{F}} + \frac{1}{2} J \mathbf{F}^{-1} : \dot{\mathbf{F}}^T \tag{B.4} \\
&= \frac{1}{2} J \mathbf{F}^{-T} : (\mathbf{L} \cdot \mathbf{F}) + \frac{1}{2} J \mathbf{F}^{-1} : (\mathbf{F}^T \cdot \mathbf{L}^T) \\
&= \frac{1}{2} J \mathbf{1} : \mathbf{L} + \frac{1}{2} J \mathbf{1} : \mathbf{L}^T \\
&= J \mathbf{1} : \mathbf{D}
\end{align*}
$$
Donc,

$$ \frac{\partial \hat{\psi}}{\partial J} \dot{j} = J \frac{\partial \hat{\psi}}{\partial J} [\mathbf{1} : \mathbf{D}] \qquad (B.5) $$

Ensuite, on développe le terme en $\dot{\mathbf{B}}$ :

$$ \dot{\mathbf{B}} = \frac{\partial \bar{\mathbf{B}}}{\partial \mathbf{B}} : \dot{\mathbf{B}} \qquad (B.6) $$

Or

$$ \frac{\partial \bar{\mathbf{B}}}{\partial \mathbf{B}} = J^{-\frac{2}{3}} \left( \mathbf{1}^{4} - \frac{1}{3} \mathbf{B} \otimes \mathbf{B}^{-1} \right) \qquad (B.7) $$

et

$$ \begin{aligned} \dot{\mathbf{B}} &= \dot{\mathbf{F}} \cdot \mathbf{F}^T + \mathbf{F} \cdot \dot{\mathbf{F}}^T \\ &= \mathbf{L} \cdot \mathbf{F} \cdot \mathbf{F}^T + \mathbf{F} \cdot \mathbf{F}^T \cdot \mathbf{L}^T \qquad (B.8) \\ &= \mathbf{L} \cdot \mathbf{B} + \mathbf{B} \cdot \mathbf{L}^T \end{aligned} $$

On utilise la relation générale (HOLZAPFEL 2000),

$$ \mathbf{A} \otimes \mathbf{B} : \mathbf{C} = \mathbf{A} (\mathbf{B} : \mathbf{C}) \qquad (B.9) $$

L'expression de $\dot{\mathbf{B}}$ devient alors :

$$ \begin{aligned} \dot{\mathbf{B}} &= J^{-\frac{2}{3}} \left( \mathbf{1}^{4} - \frac{1}{3}\mathbf{B} \otimes \mathbf{B}^{-1} \right) : \left( \mathbf{L} \cdot \mathbf{B} + \mathbf{B} \cdot \mathbf{L}^T \right) \\ &= J^{-\frac{2}{3}} \mathbf{1}^{4} : \left( \mathbf{L} \cdot \mathbf{B} + \mathbf{B} \cdot \mathbf{L}^T \right) - \frac{1}{3} \mathbf{B} \otimes \mathbf{B}^{-1} : \left( \mathbf{L} \cdot \mathbf{B} + \mathbf{B} \cdot \mathbf{L}^T \right) \\ &= J^{-\frac{2}{3}} \left( \mathbf{L} \cdot \mathbf{B} + \mathbf{B} \cdot \mathbf{L}^T \right) - \frac{1}{3} \left[ \mathbf{B}(\mathbf{B}^{-1} : \mathbf{L} \cdot \mathbf{B}) + \mathbf{B}(\mathbf{B}^{-1} : \mathbf{B} \cdot \mathbf{L}^T) \right] \qquad (B.10) \\ &= J^{-\frac{2}{3}} \left( \mathbf{L} \cdot \mathbf{B} + \mathbf{B} \cdot \mathbf{L}^T \right) - \frac{1}{3} \left[ \mathbf{B}(\mathbf{1} : \mathbf{L}) + \mathbf{B}(\mathbf{1} : \mathbf{L}^T) \right] \\ &= J^{-\frac{2}{3}} \left[ \mathbf{L} \cdot \mathbf{B} + \mathbf{B} \cdot \mathbf{L}^T - \frac{2}{3} \mathbf{B}(\mathbf{1} : \mathbf{D}) \right] \\ &= \mathbf{L} \cdot \bar{\mathbf{B}} + \bar{\mathbf{B}} \cdot \mathbf{L}^T - \frac{2}{3} \bar{\mathbf{B}}(\mathbf{1} : \mathbf{D}) \end{aligned} $$

Pour la suite on a besoin de l'expression suivante (MEO 2012) :

$$ [\mathbf{A} : \mathbf{B}] [\mathbf{C} : \mathbf{D}] = \mathbf{A} : [\mathbf{B} \otimes \mathbf{C}] : \mathbf{D} \qquad (B.11) $$
En écrivant le premier terme de $\dot{\psi}$ on a :

$$ 
\begin{aligned}
\frac{\partial \psi}{\partial \bar{\mathbf{B}}} : \dot{\bar{\mathbf{B}}} &= \frac{\partial \psi}{\partial \bar{\mathbf{B}}} : \left[ \mathbf{L} \cdot \bar{\mathbf{B}} + \bar{\mathbf{B}} \cdot \mathbf{L}^T - \frac{2}{3} \bar{\mathbf{B}} (\mathbf{1} : \mathbf{D}) \right] \\
&= \left[ 2 \bar{\mathbf{B}} \cdot \frac{\partial \psi}{\partial \bar{\mathbf{B}}} \right] : \mathbf{D} - \left[ \left[ 2 \bar{\mathbf{B}} \cdot \frac{\partial \psi}{\partial \bar{\mathbf{B}}} \right] : \mathbf{1} \right] \left[ \frac{1}{3} \mathbf{1} : \mathbf{D} \right] \tag{B.12} \\
&= \left[ 2 \bar{\mathbf{B}} \cdot \frac{\partial \psi}{\partial \bar{\mathbf{B}}} \right] : \mathbf{D} - \left[ 2 \bar{\mathbf{B}} \cdot \frac{\partial \psi}{\partial \bar{\mathbf{B}}} \right] : \frac{1}{3} [\mathbf{1} \otimes \mathbf{1}] : \mathbf{D} \\
&= \left[ 2 \bar{\mathbf{B}} \cdot \frac{\partial \psi}{\partial \bar{\mathbf{B}}} \right] : \mathbb{P}_1 : \mathbf{D}
\end{aligned}
$$

Avec

$$ \mathbb{P}_1 = \mathbf{1} - \frac{1}{3}\mathbf{1} \otimes \mathbf{1} \qquad (B.13) $$

L'égalité de dissipation s'écrit donc :

$$ \left( \sigma - 2\rho_0 J^{-1} \bar{\mathbf{B}} \cdot \frac{\partial \bar{\psi}}{\partial \bar{\mathbf{B}}} : \mathbb{P}_1 - \rho_0 \frac{\partial \hat{\psi}}{\partial J} \mathbf{1} \right) : \mathbf{D} = 0 \qquad (B.14) $$

Comme $\mathbf{D}$ ne peut être nul, l'autre partie du produit doublement contracté doit être nulle et on peut alors exprimer $\sigma$ :

$$ \sigma = 2\rho_0 J^{-1} \bar{\mathbf{B}} \cdot \frac{\partial \bar{\psi}}{\partial \bar{\mathbf{B}}} : \mathbb{P}_1 + \rho_0 \frac{\partial \hat{\psi}}{\partial J} \mathbf{1} \qquad (B.15) $$
## III.A Développement de la loi de comportement

Dans le cadre dissipatif compressible, l'énergie libre est exprimée de la sorte :

$$ \psi = \bar{\psi}_{eq}(\bar{\mathbf{B}}) + \bar{\psi}_{neq}(\bar{\mathbf{B}}_e) + \hat{\psi}_{eq}(J) + \hat{\psi}_{neq}(J_e) \quad (C.1) $$

Sa dérivée est alors de la forme :

$$ \dot{\psi} = \frac{\partial \bar{\psi}_{eq}}{\partial \bar{\mathbf{B}}} : \dot{\bar{\mathbf{B}}} + \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} : \dot{\bar{\mathbf{B}}} + \frac{\partial \hat{\psi}_{eq}}{\partial J} j + \frac{\partial \hat{\psi}_{neq}}{\partial J_e} j_e \quad (C.2) $$

Le premier terme de (C.2), comme démontré dans l'annexe B, s'écrit :

$$ \frac{\partial \psi_{eq}}{\partial \bar{\mathbf{B}}} : \dot{\bar{\mathbf{B}}} = \left[ 2\bar{\mathbf{B}} \cdot \frac{\partial \psi_{eq}}{\partial \bar{\mathbf{B}}} \right] : \mathbb{P}_1 : \mathbf{D} \quad (C.3) $$

Le troisième terme quant à lui devient :

$$ \frac{\partial \hat{\psi}_{eq}}{\partial J} j = J \frac{\partial \hat{\psi}_{eq}}{\partial J} [\mathbf{1} : \mathbf{D}] \quad (C.4) $$

Pour développer le quatrième terme, $\dot{J}_e$ devient :

$$ \begin{align} \dot{J}_e &= (JJ_a^{-1}) \cdot JJ_a^{-1} [\mathbf{1} : \mathbf{D}] - JJ_a^{-2} \dot{j}_a \nonumber \\ &= J_e [\mathbf{1} : \mathbf{D}] - JJ_a^{-1} [\mathbf{1} : \mathbf{D}_a] \tag{C.5} \\ &= J_e [\mathbf{1} : \mathbf{D}] - J_e [\mathbf{1} : \mathbf{D}_a^o] \nonumber \end{align} $$

donc,

$$ \frac{\partial \hat{\psi}_{neq}}{\partial J_e} \dot{J}_e = J_e \frac{\partial \hat{\psi}_{neq}}{\partial J_e} [\mathbf{1} : \mathbf{D}] - J_e \frac{\partial \hat{\psi}_{neq}}{\partial J_e} [\mathbf{1} : \mathbf{D}_a^o] \quad (C.6) $$

$\mathbf{D}_a^o$ est une mesure objective de la vitesse de déformation anélastique.

Ensuite pour écrire $\dot{\bar{\mathbf{B}}}_e$, on exprime d'abord $\dot{\bar{\mathbf{F}}}_e$. Pour cela on considère également que $\bar{\mathbf{F}}_a =$
$J_a^{-1/3} \mathbf{F}_a$.

On exprime $\dot{\mathbf{F}}_e$ :

$$
\begin{align*}
\dot{\mathbf{F}}_e &= [J^{-\frac{1}{3}} \mathbf{F} \cdot \bar{\mathbf{F}}_a^{-1}] \cdot \\
&= -\frac{1}{3} J^{-\frac{1}{3}} [\mathbf{1} : \mathbf{D}] \mathbf{F} \cdot \bar{\mathbf{F}}_a^{-1} + J^{-\frac{1}{3}} \dot{\mathbf{F}} \cdot \bar{\mathbf{F}}_a^{-1} + J^{-\frac{1}{3}} \cdot \mathbf{F} \cdot \dot{\bar{\mathbf{F}}}_a^{-1} \tag{C.7} \\
&= -\frac{1}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{F}}_e + \mathbf{L} \cdot \bar{\mathbf{F}}_e - \bar{\mathbf{F}}_e \cdot \bar{\mathbf{F}}_a \cdot \bar{\mathbf{F}}_a^{-1} \cdot \bar{\mathbf{L}}_a \\
&= -\frac{1}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{F}}_e + \mathbf{L} \cdot \bar{\mathbf{F}}_e - \bar{\mathbf{F}}_e \cdot \bar{\mathbf{L}}_a
\end{align*}
$$

Or

$$
\begin{align*}
\bar{\mathbf{L}}_a &= \dot{\bar{\mathbf{F}}}_a \cdot \bar{\mathbf{F}}_a^{-1} = (J_a^{-1/3} \mathbf{F}_a + J_a^{-1/3} \dot{\mathbf{F}}_a) \cdot J_a^{1/3} \mathbf{F}_a^{-1} \\
&= -\frac{1}{3} J_a^{1/3} J_a^{-4/3} \dot{\mathbf{J}}_a \mathbf{F}_a^{-1} + J_a^{-1/3} J_a^{1/3} \dot{\mathbf{F}}_a \cdot \mathbf{F}_a^{-1} \\
&= -\frac{1}{3} J_a^{-1} \dot{\mathbf{J}}_a \mathbf{1} + \mathbf{L}_a \cdot \mathbf{F}_a \cdot \mathbf{F}_a^{-1} \tag{C.8} \\
&= -\frac{1}{3} J_a^{-1} \dot{\mathbf{J}}_a \mathbf{1} + \mathbf{L}_a \\
&= -\frac{1}{3} [\mathbf{1} : \mathbf{D}_a^\circ] \mathbf{1} + \mathbf{L}_a
\end{align*}
$$

$$
\operatorname{car} \dot{J}_a = J_a[\mathbf{1} : \mathbf{D}_a] = J_a[\mathbf{1} : \mathbf{D}_a^\circ].
$$

$$
\dot{\mathbf{F}}_e = -\frac{1}{3}[\mathbf{1} : \mathbf{D}] \bar{\mathbf{F}}_e + \frac{1}{3}[\mathbf{1} : \mathbf{D}_a^\circ] \bar{\mathbf{F}}_e + \mathbf{L} \cdot \bar{\mathbf{F}}_e - \bar{\mathbf{F}}_e \cdot \mathbf{L}_a \quad (C.9)
$$

L'expression de $\dot{\bar{\mathbf{B}}}_e$ recherchée s'écrit :

$$
\begin{align}
\dot{\bar{\mathbf{B}}}_e &= \dot{\bar{\mathbf{F}}}_e \cdot \bar{\mathbf{F}}_e^T + \bar{\mathbf{F}}_e \cdot \dot{\bar{\mathbf{F}}}_e^T \nonumber \\
&= -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{B}}_e + \frac{2}{3} [\mathbf{1} : \mathbf{D}_a^\circ] \bar{\mathbf{B}}_e + \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2 \bar{\mathbf{F}}_e \cdot \mathbf{D}_a \cdot \bar{\mathbf{F}}_e^T \tag{C.10} \\
&= -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{B}}_e + \frac{2}{3} [\mathbf{1} : \mathbf{D}_a^\circ] \bar{\mathbf{B}}_e + \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2 \bar{\mathbf{V}}_e \cdot \bar{\mathbf{R}}_e \cdot \mathbf{D}_a \cdot \bar{\mathbf{R}}_e^T \cdot \bar{\mathbf{V}}_e \nonumber
\end{align}
$$

La décomposition polaire de $\mathbf{F}_e$ s'écrit aussi $\mathbf{F}_e = \mathbf{V}_e \cdot \mathbf{R}_e = J_e^{\frac{1}{3}} \bar{\mathbf{F}}_e = J_e^{\frac{1}{3}} \bar{\mathbf{V}}_e \cdot \bar{\mathbf{R}}_e$. Et comme la décomposition est unique, on peut dire que $\mathbf{V}_e = J_e^{\frac{1}{3}} \bar{\mathbf{V}}_e$ et $\mathbf{R}_e = \bar{\mathbf{R}}_e$. Donc comme $\mathbf{D}_a^\circ = \mathbf{R}_e \cdot \mathbf{D}_a \cdot \mathbf{R}_e^T$ :

$$
\begin{align}
\dot{\bar{\mathbf{B}}}_e &= -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{B}}_e + \frac{2}{3} [\mathbf{1} : \mathbf{D}_a^\circ] \bar{\mathbf{B}}_e + \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2 \bar{\mathbf{V}}_e \cdot \mathbf{R}_e \cdot \mathbf{D}_a \cdot \mathbf{R}_e^T \cdot \bar{\mathbf{V}}_e \tag{C.11} \\
&= -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{B}}_e + \frac{2}{3} [\mathbf{1} : \mathbf{D}_a^\circ] \bar{\mathbf{B}}_e + \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2 \bar{\mathbf{V}}_e \cdot \mathbf{D}_a^\circ \cdot \bar{\mathbf{V}}_e \nonumber
\end{align}
$$
Le deuxième terme de $\dot{\psi}$ devient alors :

$$
\begin{align*}
\frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} : \dot{\mathbf{B}}_e &= \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} : \left[ -\frac{2}{3} [\mathbf{1} : \mathbf{D}] \bar{\mathbf{B}}_e + \frac{2}{3} [\mathbf{1} : \mathbf{D}_a^o] \bar{\mathbf{B}}_e \right. \\
&\qquad \left. + \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2 \bar{\mathbf{V}}_e \cdot \mathbf{D}_a^o \cdot \bar{\mathbf{V}}_e \right] \\
&= \left[ 2 \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \right] : \mathbf{D} - \left[ \left[ 2 \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \right] : \mathbf{1} \right] \left[ \frac{1}{3} \mathbf{1} : \mathbf{D} \right] \\
&\quad - \left[ 2 \bar{\mathbf{V}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \cdot \bar{\mathbf{V}}_e \right] : \mathbf{D}_a^o + \left[ \left[ 2 \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \right] : \mathbf{1} \right] \left[ \frac{1}{3} \mathbf{1} : \mathbf{D}_a^o \right] \\
&= \left[ 2 \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \right] : \mathbf{D} - \left[ 2 \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \right] : \frac{1}{3} [\mathbf{1} \otimes \mathbf{1}] : \mathbf{D} \\
&\quad - \left[ 2 \bar{\mathbf{V}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \cdot \bar{\mathbf{V}}_e \right] : \mathbf{D}_a^o + \left[ 2 \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \right] : \frac{1}{3} [\mathbf{1} \otimes \mathbf{1}] : \mathbf{D}_a^o \\
&= \left[ 2 \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \right] : \mathbb{P}_1 : \mathbf{D} - \left[ 2 \bar{\mathbf{V}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \cdot \bar{\mathbf{V}}_e \right] : \mathbb{P}_1 : \mathbf{D}_a^o
\end{align*}
\tag{C.12}
$$

L'inégalité de dissipation s'écrit donc :

$$
\begin{align}
& - \rho \dot{\psi} + \sigma : \mathbf{D} \ge 0 \nonumber \\
& \left[ -2\rho_0 J^{-1} \bar{\mathbf{B}} \cdot \frac{\partial \bar{\psi}_{eq}}{\partial \bar{\mathbf{B}}} : \mathbb{P}_1 - 2\rho_0 J^{-1} \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} : \mathbb{P}_1 - \rho_0 \frac{\partial \hat{\psi}_{eq}}{\partial J} \mathbf{1} - \rho_0 J^{-1} J_e \frac{\partial \hat{\psi}_{neq}}{\partial J_e} \mathbf{1} + \sigma \right] : \mathbf{D} \nonumber \\
& \quad + \left[ \rho_0 J^{-1} J_e \frac{\partial \hat{\psi}_{neq}}{\partial J_e} \mathbf{1} + 2\rho_0 J^{-1} \bar{\mathbf{V}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \cdot \bar{\mathbf{V}}_e : \mathbb{P}_1 \right] : \mathbf{D}_a^o \ge 0 \tag{C.13}
\end{align}
$$

À partir de l'équation (C.8) et de la définition même de **D**, il est possible d'écrire :

$$
\begin{align}
\mathbf{D} &= \bar{\mathbf{D}} + \frac{1}{3} [\mathbf{1} : \mathbf{D}] \mathbf{1} \tag{C.14} \\
\mathbf{D}_a^o &= \bar{\mathbf{D}}_a^o + \frac{1}{3} [\mathbf{1} : \mathbf{D}_a^o] \mathbf{1} \nonumber
\end{align}
$$

Si on exprime $\hat{\mathbf{D}}$ et $\hat{\mathbf{D}}_a^o$ en fonction de $\hat{\mathbf{F}}_a$, on obtient les expressions suivantes :

$$
\hat{\mathbf{D}} = \hat{\mathbf{F}} \hat{\mathbf{F}}^{-1} = j^{1/3} J^{-1/3} \mathbf{1} = \frac{1}{3} J^{-2/3} j J^{-1/3} \mathbf{1} = \frac{1}{3} [\mathbf{1} : \mathbf{D}] \mathbf{1} \quad (C.15)
$$

De même :

$$
\hat{\mathbf{D}}_a^o = \frac{1}{3} [\mathbf{1} : \mathbf{D}_a^o] \mathbf{1} \qquad (C.16)
$$

Grâce à la décomposition multiplicative de **F** en $\bar{\mathbf{F}}$ et $\hat{\mathbf{F}}$, il vient donc naturellement une dé-
composition de **D** et $\mathbf{D}_a^o$ en tenseurs déviatorique et sphérique.
En introduisant cette expression dans l'inégalité de dissipation (C.13), on obtient :

$$
\begin{equation}
\begin{aligned}
& \left[ -2\rho_0 J^{-1} \bar{\mathbf{B}} \cdot \frac{\partial \bar{\psi}_{eq}}{\partial \bar{\mathbf{B}}} : \mathbb{P}_1 - 2\rho_0 J^{-1} \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} : \mathbb{P}_1 + \bar{\boldsymbol{\sigma}} \right] : \bar{\mathbf{D}} \\
& \quad + \left[ -\rho_0 \frac{\partial \hat{\psi}_{eq}}{\partial J} \mathbf{1} - \rho_0 J^{-1} J_e \frac{\partial \hat{\psi}_{neq}}{\partial J_e} \mathbf{1} + \hat{\boldsymbol{\sigma}} \right] : \hat{\mathbf{D}} \\
& \quad + \left[ 2\rho_0 J^{-1} \bar{\mathbf{V}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \cdot \bar{\mathbf{V}}_e : \mathbb{P}_1 \right] : \bar{\mathbf{D}}_a^o \\
& \quad + \left[ \rho_0 J^{-1} J_e \frac{\partial \hat{\psi}_{neq}}{\partial J_e} \mathbf{1} \right] : \hat{\mathbf{D}}_a^o \geq 0
\end{aligned}
\tag{C.17}
\end{equation}
$$

Après avoir introduit le potentiel de dissipation dans le cadre des matériaux standards
généralisés, comme expliqué dans la partie III.C.1.b, on arrive aux lois de comportements

$$
\bar{\sigma} = 2\rho_0 J^{-1} \bar{\mathbf{B}} \cdot \frac{\partial \bar{\psi}_{eq}}{\partial \bar{\mathbf{B}}} : \mathbb{P}_1 + 2\rho_0 J^{-1} \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} : \mathbb{P}_1 + \frac{\partial \bar{\varphi}_0}{\partial \bar{\mathbf{D}}} \quad (C.18)
$$

$$
\hat{\boldsymbol{\sigma}} = \rho_0 \frac{\partial \hat{\psi}_{eq}}{\partial J} \mathbf{1} + \rho_0 J^{-1} J_e \frac{\partial \hat{\psi}_{neq}}{\partial J_e} \mathbf{1} + \frac{\partial \hat{\varphi}_0}{\partial \hat{\mathbf{D}}} \quad (C.19)
$$

avec $\boldsymbol{\sigma} = \bar{\boldsymbol{\sigma}} + \hat{\boldsymbol{\sigma}}$,

et aux lois d'évolution portant sur les parties déviatorique et sphérique :

$$
2\rho_0 J^{-1} \bar{\mathbf{V}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} \cdot \bar{\mathbf{V}}_e : \mathbb{P}_1 = \frac{\partial \bar{\varphi}_a}{\partial \bar{\mathbf{D}}_a^o} \quad (C.20)
$$

ou si $\psi_{neq}$ est isotrope de $\bar{\mathbf{B}}_e$ :

$$
2\rho_0 J^{-1} \bar{\mathbf{B}}_e \cdot \frac{\partial \bar{\psi}_{neq}}{\partial \bar{\mathbf{B}}_e} : \mathbb{P}_1 = \frac{\partial \bar{\varphi}_a}{\partial \bar{\mathbf{D}}_a^o} \quad (C.21)
$$

$$
\rho_0 J^{-1} J_e \frac{\partial \hat{\psi}_{neq}}{\partial J_e} \mathbf{1} = \frac{\partial \hat{\varphi}_a}{\partial \hat{\mathbf{D}}_a^o} \quad (C.22)
$$

Ces lois d'évolution peuvent être écrites sous forme différentielle.
# Annexe D

## Identification de paramètres à différentes températures

Un des objectifs secondaires de la thèse était de prendre en compte l'influence de la température sur le comportement mécanique. Un couplage thermomécanique aurait été de trop grande envergure au vu de tout ce que est déjà développé (mais reste néanmoins une perspective intéressante et envisageable de ce travail).

À défaut, conscients de l'influence de la température sur la réponse mécanique, nous allons identifier les paramètres du modèle à chaque température testée (qui sera donc supposée temporellement constante et homogène dans le matériau) : -20, 24 et 60°C. Une fois les valeurs identifiées, nous proposerons une relation empirique simple de dépendance à la température des paramètres matériaux comme cela a pu être fait dans MEO 2012, 2000 et DELATTRE 2014.

L'identification ici proposée sur un unique chargement de traction simple dans le cas d'un modèle hyperélastique et d'un chargement monotone. En effet, le méta-modèle obtenu en couplant POD et krigeage est construit avec les mêmes snapshots contrainte-elongation que dans le chapitre 4. Il suffit de projeter sur la variété de courbes les snapshots expérimentaux obtenus à chaque température.

La Figure D.1a présente les trois snapshots à -20°C, 24°C et 60°C, superposés aux snapshots du plan d'expérience et la Figure D.1b, leur représentation dans la base de POD. On voit que les trois snapshots se trouvent dans la manifold qui a 3 dimensions. Les valeurs de paramètres auront des chances d'être identifiées.

La Figure D.2 compare pour chaque température les courbes expérimentales et identifiées (pointillés) obtenues pour le R23. Les résultats sont en bonne adéquation avec les essais. L'évolution de la valeur des paramètres en fonction de la température est tracée sur la Figure D.3a. Des courbes de tendance linéaires ou polynomiales (ordre 2) sont établies et leurs équations figurent sur le graphique. À partir de ces expressions, on propose de vérifier que pour d'autres températures les résultats sont cohérents. On prend deux températures de la plage d'application du matériau, 30 et 70°C. Les prédicitions pour ces deux températures sont tracées sur la Figure D.3b. On voit que les courbes obtenues (en pointillés) ont une forme cohérente avec celles identifiées aux autres températures. De pus, l'évolution selon l'élasticité entropique est bien respectée. Les résultats sont donc concluants pour la traction uniaxiale.

Normalement, les paramètres $C_{ij}$ devraient être identifiés sur des essais de compression confinée en température également. Cependant, comme évoqué dans d'autres parties de ce mémoire,
<figure><img src="image_96.png" /><figcaption>FIGURE D.1 – Snapshots servant à l'identification du comportement en traction uniaxiale et leur représentation dans l'espace réduit de la POD. Correspondance couleurs/température : bleu -20°C, vert 24°C et rouge 60°C.</figcaption></figure>

<figure><img src="image_97.png" /><figcaption>FIGURE D.2 – Courbes contrainte-déformation en traction sur le R23 à différentes températures et leurs courbes identifiées.</figcaption></figure>

ces essais ont été jugés comme potentiellement compliqués à utiliser/analyser : frottement sur la paroi, jue initial entre l'échantillon et l'enceinte..., ces points seront a priori tout aussi (si ce n'est plus) présents à températures variables, la dilation thermique n'étant a priori pas contrôlée. Les variations de contrainte de ce cas seront non maîtrisées.
<figure><img src="image_98.png" /><figcaption>FIGURE D.3 – Prédictions du comportement du R23 en traction uniaxiale à différentes températures.</figcaption></figure>
# Annexe E

## Influence du pourcentage de microsphères sur le comportement mécanique

Un des objectifs de la thèse était d'avoir une loi de comportement dépendante du pourcentage de microsphères. C'est partiellement le cas de par l'utilisation d'énergies libres volumiques fixant la position de l'asymptote de comportement incompressible lors de sollicitations sphériques où toutes les microsphères finiraient par être fermées. On cherche, pour compléter cela, à établir une relation empirique entre les paramètres identifiés au Chapitre 4 et le pourcentage de microsphères.

Cela commence avec les paramètres de la partie volumique du modèle. Nous disposons des matériaux R23 et R33, mais aussi d'un mélange obtenu par une heureuse erreur de formulation qui a donné un grade avec 10% de microsphères (R10). Après avoir identifié leurs paramètres, il est possible d'établir une relation linéaire ou polynomiale avec le pourcentage de microsphères, comme présenté sur la Figure E.1.

<figure><img src="image_99.png" /><figcaption>FIGURE E.1 – Mise en place de relations entre les paramètres de la contrainte hydrostatique et le pourcentage de microsphères</figcaption></figure>

Cela permet de prédire pour d'autres pourcentages de microsphères, les réponses en compression hydrostatique, comme présentées sur la Figure E.2 pour 16% (R16) et 28% (R28).