# Modèle viscoélastique-viscoplastique couplé avec endommagement non-local pour polymères semi-cristallins

Romain Balieu, Franck Lauro, Bruno Bennani, Tsukatada Matsumoto,
Ernesto Mottola

► To cite this version:

Romain Balieu, Franck Lauro, Bruno Bennani, Tsukatada Matsumoto, Ernesto Mottola. Modèle viscoélastique-viscoplastique couplé avec endommagement non-local pour polymères semi-cristallins. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717085

HAL Id: hal-01717085

https://hal.science/hal-01717085v1

Submitted on 25 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Modèle viscoélastique-viscoplastique couplé avec endommagement non-local pour polymères semi-cristallins.

Romain BALIEU <fnref n="1" />*, Franck LAURO <fnref n="1" />, Bruno BENNANI <fnref n="1" />, Tsukatada MATSUMOTO <fnref n="2" />, Ernesto MOTTOLA <fnref n="2" />

**Résumé** — Les polymères semi-cristallins sont très utilisés dans de nombreux secteurs industriels et leurs comportements complexes requièrent des modèles numériques précis pour la simulation éléments finis. Dans ce travail, un modèle de comportement non-associatif viscoélastique-viscoplastique couplé a un modèle d'endommagement non-local est implémenté dans un code éléments finis commercial afin de modéliser le comportement d'un polypropylène chargé de talcs.

## 1 Introduction

Différents modèles constitutifs ont été proposés dans la littérature afin de modéliser le comportement complexe des polymères. Ces modèles peuvent être rangés en deux catégories : les modèles physiques où la micro-structure du matériau est prise en compte pour décrire le comportement macroscopique [1] et les modèles phénoménologiques où les discontinuités du matériau, à l'échelle micro-structurale, sont homogénéisées dans un volume élémentaire représentatif. En suivant cette approche, des modèles de comportement élasto-viscoplastiques basés sur le concept de « overstress » (VBO) [2] utilisant la théorie des variables d'états unifiées ont été étendus pour la modélisation des matériaux polymères [3]. L'addition de charges minérales dans la matrice semi-cristalline accentue le phénomène de cavitation dû à la décohésion de l'interface matrice-particule dans ces types de matériaux. Dans ce cas, la déformation viscoplastique du matériau est accompagnée du processus d'endommagement sous la forme de nucléation, croissance et coalescence de cavités. Différents modèles d'endommagements appliqués aux matériaux polymères ont donc été proposés afin de représenter ce phénomène [4, 5]. L'endommagement présent dans la déformation de ces types de matériaux provoque un adoucissement des contraintes conduisant au phénomène de localisation des déformations entraînant une solution numérique dépendante de la discrétisation de la structure. Les modèles de comportement non-locaux ont été introduits dans la littérature afin de palier à cet effet non-désirable de dépendance au maillage [6].

Dans ce travail, un modèle de comportement viscoélastique-viscoplastique couplé à un potentiel de dissipation non-associatif et une surface de charge non-symétrique est proposé afin de modéliser les phénomènes observés expérimentalement. Les équations constitutives du modèle sont inscrites dans le cadre des grandes transformations en utilisant un formalisme hypoélastique. Les propriétés intéressantes du tenseur des taux de rotation logarithmique sont utilisées dans ce modèle afin de faire le lien entre les contraintes de Cauchy énergiquement conjuguées aux déformations de Henky. Afin d'obtenir une solution indépendante du maillage, un modèle d'endommagement non-local est utilisé dans le modèle développé.

## 2 Modèle constitutif

### 2.1 Formulation corotationnelle du modèle viscoélastique

L'objectif principal en mécanique des milieux continus utilisant les grandes déformations est de respecter l'objectivité au cours du temps. La formulation d'un modèle élastique ou viscoélastique-viscoplastique

<footnote n="1">¹ LAMIH, Université de Valenciennes et du Hainaut-Cambrésis, UMR 8201</footnote>

<footnote n="2">² TOYOTA MOTOR EUROPE, Bruxelles</footnote>

<footnote>* Auteur correspondent: romain.balieu@univ-valenciennes.fr</footnote>
doit satisfaire l'axiome d'objectivité qui déclare que « les équations constitutives doivent être invariantes durant un changement de repère référentiel »[7]. Afin de satisfaire cette indifférence matérielle, de nombreux taux objectifs ont été introduits dans la littérature pour remplacer la dérivée temporelle matérielle. Le taux objectif d'un tenseur Eulérien du second ordre a est défini par

$$ \mathring{\mathbf{a}}^* = \mathbf{a} - \mathbf{\Omega}^* \mathbf{a} + \mathbf{a} \mathbf{\Omega}^*, \qquad (1) $$

où $\mathbf{a}$ est la dérivée temporelle matérielle de $\mathbf{a}$ et $\mathbf{\Omega}^*$ est le tenseur des taux de rotation, i.e., un tenseur du second ordre antisymétrique dépendant du temps qui défini un repère tournant $(\cdot)^*$ relatif à un repère fixe. Le tenseur des taux de rotation est défini par son tenseur orthogonal propre $\mathbf{Q}^*$ tel que

$$ \mathbf{\Omega}^* = \dot{\mathbf{Q}}^{*T} \mathbf{Q}^* = -\mathbf{Q}^{*T} \dot{\mathbf{Q}}^*. \qquad (2) $$

Dans le repère tournant $(\cdot)^*$, un tenseur Eulérien du second ordre $\mathbf{a}$ devient $\mathbf{a}^* = \mathbf{Q}^* \mathbf{a} \mathbf{Q}^{*T}$, le retour au repère spatial du tenseur $\mathbf{a}$ est alors réalisé par l'opération inverse, i.e., $\mathbf{a} = \mathbf{Q}^{*T} \mathbf{a}^* \mathbf{Q}^*$. La forme générale d'une relation hypoélastique est donnée par

$$ \mathring{\mathbf{T}}^* = \mathcal{F}(\mathbf{T}, \mathbf{D}), \qquad (3) $$

où $\mathbf{T}$ est un tenseur Eulérien de contraintes du deuxième ordre symétrique et différentiable. $\mathbf{D}$ est le tenseur des taux de déformation défini par

$$ \mathbf{D} = \frac{1}{2} (\mathbf{L} + \mathbf{L}^T), \qquad (4) $$

où $\mathbf{L}$ est le tenseur des gradients de vitesse calculé à partir du tenseur des gradients de déformation $\mathbf{F}$, i.e.,

$$ \mathbf{L} = \dot{\mathbf{F}} \mathbf{F}^{-1}. \qquad (5) $$

Dans le modèle constitutif, la relation directe entre le taux de déformation et le tenseur de déformation de Henky (logarithmique) proposée par Xiao et al. [8] est utilisée. Dans les travaux effectués par les précédents auteurs, il a été prouvé que le taux de déformation $\mathbf{D}$ est identique au taux corotationnel du tenseur de déformation de Henky<fnref n="2" /> $\mathbf{h}$ en utilisant le tenseur des taux de rotation logarithmique $\mathring{\mathbf{\Omega}}^{\text{log}}$, i.e.,

$$ \mathring{\mathbf{h}}^{\text{log}} = \mathbf{h} - \mathbf{\Omega}^{\text{log}} \mathbf{h} + \mathbf{h} \mathbf{\Omega}^{\text{log}} = \mathbf{D}. \qquad (7) $$

Le tenseur des taux de rotation logarithmique est défini par

$$ \mathbf{\Omega}^{\text{log}} = \mathbf{W} + \sum_{i,j=1, i \neq j}^{m} \left[ \left( \frac{1 + (b_i/b_j)}{1 + (b_i/b_j)} + \frac{2}{\log(b_i/b_j)} \right) \mathbf{P}_i \mathbf{D} \mathbf{P}_j \right], \qquad (8) $$

où $\mathbf{W}$ est la partie antisymétrique du tenseur des gradients de vitesse, i.e. $\mathbf{W} = \frac{1}{2}(\mathbf{L} - \mathbf{L}^T)$, $b_i = \lambda_i^2$ sont les valeurs propres du tenseur Eulérien des déformations de Cauchy-Green $\mathbf{B}$ et $\mathbf{P}_i$ les vecteurs propres associés aux valeurs propres $\lambda_i > 0$ de $\mathbf{V}$. De plus, la notion de mesure Lagrangienne des déformations énergétiquement conjuguée aux contraintes proposée par Hill [9] a été étendue aux mesures Eulériennes [8]. Il en résulte que l'unique repère tournant où les contraintes de Cauchy sont énergétiquement conjuguées aux déformations de Henky est le repère tournant défini par le tenseur des taux de rotation logarithmique $\mathbf{\Omega}^{\text{log}}$, i.e.,

$$ \boldsymbol{\sigma}: \mathring{\mathbf{h}}^{\text{log}} = \boldsymbol{\sigma}: \mathbf{D}. \qquad (9) $$

$$ \mathbf{h} = \frac{1}{2} \log(\mathbf{B}), \qquad (6) $$

où $\mathbf{B}$ est le tenseur Eulérien des déformations de Cauchy-Green, i.e. $\mathbf{B} = \mathbf{F}^T \mathbf{F} = \mathbf{V}^2$, avec $\mathbf{V}$ le tenseur Eulérien symétrique et défini positif des taux de déformation.

<footnote n="1">1. Le symbole $(\cdot)^T$ désigne la transposition.</footnote>

<footnote n="2">2. Le tenseur des déformations de Henky est défini par</footnote>
Dans le modèle proposé, la décomposition additive du taux de déformation $\mathbf{D}$ en une partie viscoélastique $\mathbf{D}^{\text{ve}}$ et viscoplastique $\mathbf{D}^{\text{vp}}$ est assumée, i.e., $\mathbf{D} = \mathbf{D}^{\text{ve}} + \mathbf{D}^{\text{vp}}$. La forme générale hypoélastique du modèle de comportement proposée dans ces travaux est donc définie par

$$ \mathring{\boldsymbol{\sigma}}^{\log} = \mathcal{F}(\boldsymbol{\sigma}, \mathbf{D}^{\text{ve}}) = \mathcal{F}(\boldsymbol{\sigma}, \mathbf{D} - \mathbf{D}^{\text{vp}}) = \mathcal{F}\left(\boldsymbol{\sigma}, \mathring{\mathbf{h}}^{\log} - (\mathring{\mathbf{h}}^{\text{vp}})^{\log}\right). \quad (10) $$

L'intégration du modèle hypoélastique (Eq. (10)) dans l'intervalle de temps $[0, t]$ peut être réalisé aisément en utilisant l'intégration corotationnelle $\int_{corot}$ d'un tenseur Eulérien du second ordre, défini par

$$ \int_{corot} \mathbf{A} \mathrm{d}t = \mathbf{Q}^{*\mathrm{T}} \left[ \int_{0}^{t} \mathbf{Q}^{*} \mathbf{A} \mathbf{Q}^{*\mathrm{T}} \mathrm{d}t \right] \mathbf{Q}^{*}. \quad (11) $$

Appliqué à un taux corotationnel $\mathring{\mathbf{a}}^*$ possédant la condition initial $\mathbf{a}|_{t=0} = \mathbf{0}$, Eq. (11) devient

$$ \int_{corot} \mathring{\mathbf{A}}^{*} \mathrm{d}t = \mathbf{Q}^{*\mathrm{T}} \left[ \int_{0}^{t} \overline{\mathbf{Q}^{*} \mathring{\mathbf{A}} \mathbf{Q}^{*\mathrm{T}}} \mathrm{d}t \right] \mathbf{Q}^{*} = \mathbf{Q}^{*\mathrm{T}} \mathbf{Q}^{*} \mathbf{A} \mathbf{Q}^{*\mathrm{T}} \mathbf{Q}^{*} = \mathbf{A}. \quad (12) $$

Dans le cas particulier de l'intégration corotationnelle logarithmique, les tenseurs des contraintes de Cauchy et des déformations viscoélastiques et viscoplastiques de Henky sont obtenus suivant

$$ \boldsymbol{\sigma} = \int_{corot} (\boldsymbol{\sigma})^{\log}, \quad (13) $$

$$ \mathbf{h}^{\text{ve}} = \int_{corot} (\mathbf{h}^{\text{ve}})^{\log} \mathrm{d}t = \int_{corot} \mathbf{D}^{\text{ve}} \mathrm{d}t, \quad (14) $$

$$ \mathbf{h}^{\text{vp}} = \int_{corot} (\mathbf{h}^{\text{vp}})^{\log} \mathrm{d}t = \int_{corot} \mathbf{D}^{\text{vp}} \mathrm{d}t. \quad (15) $$

La forme générale du modèle de comportement proposée (Eq. (10)) peut donc être écrite sous la forme suivante

$$ \boldsymbol{\sigma} = \mathcal{F}(\boldsymbol{\sigma}, \mathbf{h}^{\text{ve}}) = \mathcal{F}(\boldsymbol{\sigma}, \mathbf{h} - \mathbf{h}^{\text{vp}}). \quad (16) $$

Le modèle viscoélastique linéaire de Wiechert (modèle de Maxwell généralisé) est utilisé dans le modèle constitutif afin de reproduire la sensibilité de la rigidité du matériau à la vitesse de déformation. Le tenseur des contraintes endommagées défini par le modèle viscoélastique à un temps $t$, i.e. $\boldsymbol{\sigma}(t)$, résultant d'un incrément de déformation viscoélastique à un temps $\zeta$, i.e. $\mathbf{h}^{\text{ve}}(\zeta)$, est donné par

$$ \boldsymbol{\sigma}(t) = (1-D) \int_{-\infty}^{t} \left\{ \mathcal{L}_{\mathring{\mathbf{z}}_{\infty}}^{\text{ve}} + \sum_{i=1}^{N} \mathcal{L}_{\mathring{\mathbf{z}}_{i}}^{\text{ve}} \exp\left(-\frac{t-\zeta}{\tau_{i}}\right) \right\} : \frac{d\mathbf{h}_{\mathring{\mathbf{z}}_{\infty}}^{\text{ve}}(\zeta)}{d\zeta} \mathrm{d}\zeta. \quad (17) $$

$D$ est la variable scalaire d'endommagement isotrope, $\mathcal{L}_{\mathring{\mathbf{z}}_{\infty}}^{\text{ve}}$ et $\mathcal{L}_{\mathring{\mathbf{z}}_{i}}^{\text{ve}}$ sont respectivement les tenseurs de rigidité des ressorts compris dans les éléments Maxwell et le tenseur de rigidité à long terme, définis par

$$ \mathcal{L}_{\mathring{\mathbf{z}}_{\infty}}^{\text{ve}} = 2G_{\infty}\mathbf{I}_{\mathring{\mathbf{z}}_{\infty}} + K_{\infty}\mathbf{I}_{\mathring{\mathbf{z}}_{\infty}}\otimes\mathbf{I}_{\mathring{\mathbf{z}}_{\infty}} \quad \text{et} \quad \mathcal{L}_{\mathring{\mathbf{z}}_{i}}^{\text{ve}} = 2G_{i}\mathbf{I}_{\mathring{\mathbf{z}}_{i}} + K_{i}\mathbf{I}_{\mathring{\mathbf{z}}_{i}}\otimes\mathbf{I}_{\mathring{\mathbf{z}}_{i}}, \quad (18) $$

où $G_{\infty}$, $K_{\infty}$, $G_i$ et $K_i$ sont les modules de cisaillement et de compressibilité à long terme et instantanés, respectivement.

### 3 Formulation viscoplastique non-associative

La surface de charge donnée par Raghava [10] qui dépend du deuxième invariant du tenseur déviatorique des contraintes $J_2(\mathbf{S})$ et du premier invariant du tenseur des contraintes $I_1(\boldsymbol{\sigma})$ est utilisée dans le modèle constitutif. La surface de charge est définie par

$$ f^{\text{vp}}(\boldsymbol{\sigma}, R, D) = \frac{(\eta - 1)I_1(\boldsymbol{\sigma}) + \sqrt{(\eta - 1)^2I_1^2(\boldsymbol{\sigma}) + 12\eta J_2(\mathbf{S})}}{2\eta(1 - D)} - \boldsymbol{\sigma}_t - R(\kappa), \quad (19) $$
où η est le ratio entre la limite élastique en compression et en traction σ_c et σ_t, respectivement. D est la variable d'endommagement isotrope et R représente l'écrouissage isotrope défini par

$$R(\kappa) = Q_1 \kappa \exp(-b_1 \kappa) + Q_2(1 - \exp(-b_2 \kappa)) + b_3 \kappa^3 + b_4 \kappa^2 + b_5 \kappa. \quad (20)$$

$Q_1, Q_2, b_1, b_2, b_3, b_4$ et $b_5$ sont des paramètres matériaux. La variable d'écrouissage $\kappa$ utilisée dans Eq. (20) est la déformation viscoplastique équivalente, i.e.,

$$\kappa = \sqrt{\frac{2}{3} \underline{\mathbf{h}}^{\text{vp}} : \underline{\mathbf{h}}^{\text{vp}}}. \quad (21)$$

La déformation viscoplastique des matériaux polymères n'est pas un phénomène isochorique. La viscoplasticité non-associative est alors utilisée dans le modèle constitutif afin de modéliser la variation de volume du matériau. Un second potentiel de dissipation F est alors introduit pour définir les règles de normalité utilisées dans le modèle. Le potentiel viscoplastique utilisé est défini par

$$F(\underline{\boldsymbol{\sigma}}) = \frac{\sqrt{3J_2(\underline{\boldsymbol{S}}) + \alpha^+ \langle p \rangle^2 + \alpha^- \langle -p \rangle^2}}{1 - D}, \quad (22)$$

où $\alpha^+$ et $\alpha^-$ sont des paramètres qui définissent la dilatation et la rétractation du matériau pour des sollicitations définies par des pressions hydrostatiques positives et négatives, respectivement. $\langle \cdot \rangle$ est le symbole de Macaulet, qui est pour un scalaire $x$, donné par $(x + |x|)/2$. Avec la formulation du potentiel viscoplastique proposée, la direction de flux du tenseur des vitesses de déformations peut évoluer indépendamment avec des pressions hydrostatiques positives (dilatation) et négatives (rétractation). En suivant la règle de normalité non-associative, le tenseur des vitesses de déformations viscoplastiques est donné par

$$\underline{\mathbf{h}}^{\text{vp}} = \dot{\lambda} \frac{\partial F}{\partial \underline{\boldsymbol{\sigma}}} = \dot{\lambda} \underline{\mathbf{n}}, \quad (23)$$

où $\underline{\mathbf{n}}$ est la direction du flux viscoplastique exprimée par

$$\underline{\mathbf{n}} = \frac{1}{g(1-D)} \left( \frac{3}{2} \underline{\mathbf{S}} + \frac{1}{3} (\alpha^+ \langle p \rangle - \alpha^- \langle -p \rangle) \underline{\mathbf{I}} \right), \quad (24)$$

avec le scalaire $g$ défini tel que

$$g = \sqrt{3J_2(\underline{\boldsymbol{S}}) + \alpha^+ \langle p \rangle^2 + \alpha^- \langle -p \rangle^2}. \quad (25)$$

En substituant Eq. (21) dans la précédente définition du tenseur des vitesses de déformations viscoplastiques (Eq. (23), la vitesse de déformation équivalente $\dot{\kappa}$ devient

$$\dot{\kappa} = \sqrt{\frac{2}{3} \underline{\mathbf{h}}^{\text{vp}} : \underline{\mathbf{h}}^{\text{vp}}} = \dot{\lambda} \sqrt{\frac{2}{3} \underline{\mathbf{n}} : \underline{\mathbf{n}}} = \frac{\dot{\lambda}}{g(1-D)} \sqrt{3J_2(\underline{\boldsymbol{S}}) + \frac{2}{9}(\alpha^+ \langle p \rangle - \alpha^- \langle -p \rangle)^2}. \quad (26)$$

Les matériaux polymères sont très dépendants de la vitesse de déformation. Une formulation viscoplastique est alors utilisée dans le modèle constitutif afin de prendre en compte la sensibilité à la vitesse de déformation sur la surface de charge. Pour inclure la viscoplasticité, une contrainte visqueuse $\sigma^v$ est définie pour étendre la surface de charge statique $f$ [11]. La contrainte visqueuse représente la différence entre la surface de charge statique $f$ et dynamique $f^d$. La surface de charge dynamique est donc naturellement définie par

$$f^d = \frac{(\eta - 1)I_1(\tilde{\sigma}) + \sqrt{(\eta - 1)^2 I_1^2(\tilde{\sigma}) + 12\eta J_2(\tilde{S})}}{2\eta} - \sigma_t - R(\kappa) - \sigma^v \equiv 0. \quad (27)$$

Classiquement, dans la formulation d'un modèle viscoplastique, le multiplicateur viscoplastique $\dot{\lambda}$ est une fonction de la contrainte visqueuse $\sigma^v$. Dans le modèle constitutif, le multiplicateur viscoplastique est donné par

$$\dot{\lambda} = \frac{\eta^v}{\sqrt{\frac{2}{3} \underline{\mathbf{n}} : \underline{\mathbf{n}}}} \left\langle \frac{\sigma^v}{\sigma_y + R(\kappa)} \right\rangle^{1/n}, \quad (28)$$
où $n$ et $\eta^{vp}$ sont les paramètres de la sensibilité à la vitesse de déformation. En choisissant la contrainte visqueuse tel que $\sigma^v = \langle f \rangle$ et en combinant Eqs. (28) et (26), la contrainte visqueuse peut être réécrite tel que

$$ \tilde{\sigma}^v = (\sigma_y + R(\kappa)) \left(\frac{\dot{\kappa}}{\eta^v}\right)^n. \qquad (29) $$

Finalement en substituant Eq. (29) dans Eq. (27), la surface de charge dynamique du modèle constitutif est définie par

$$ f^d = \frac{(\eta - 1)I_1(\boldsymbol{\sigma}) + \sqrt{(\eta - 1)^2 I_1^2(\boldsymbol{\sigma}) + 12\eta J_2(\boldsymbol{S})}}{(1 - D)2\eta} - (\sigma_t + R(\kappa)) \left(1 + \frac{\dot{\kappa}}{\eta^v}\right)^n \equiv 0. \qquad (30) $$

De par sa définition, la contrainte visqueuse étant positive ou nulle, le multiplicateur viscoplastique est aussi positif ou nul pour tout les états de contraintes, i.e. $\dot{\lambda} \ge 0$. Lorsque l'état de contraintes se place dans le domaine viscoplastique, la surface de charge dynamique $f^d$ étant l'unique chemin possible, i.e., $f^d = 0$, les conditions standards de chargement/déchargement de Kuhn-Tucker sont donc réspectées, i.e.,

$$ f^d \le 0, \quad \dot{\lambda} \ge 0, \quad \dot{\lambda} f^d = 0. \qquad (31) $$

## 4 Modèle d'endommagement isotrope non-locale

Au cours de la déformation des polymères semi-cristallins où des charges minérales sont introduites, le phénomène bien connu de cavitation intervient. La cavitation est due à la décohésion de l'interface matrice-particule. Dans ce cas la déformation viscoplastique du matériau est accompagnée du processus d'endommagement sous la forme de nucléation, croissance et coalescence de cavités [12]. Dans le modèle constitutif, un modèle phénoménologique d'endommagement est introduit afin de représenter le phénomène de cavitation. Dans ces travaux, un endommagement isotrope, représenté par une variable scalaire $D$ [13] est alors introduite dans les équations constitutives. La section résistante totale d'un solide dans sa configuration réelle (i.e. endommagée) notée $S$ et la surface totale de toutes les fissures et cavités notée $S^D$ sont définies. La surface fictive non-endommagée (i.e. effective) notée $\tilde{S}$ est alors obtenue en supprimant toute les cavités et fissures qui caractérisent l'endommagement. Dans le cas d'un endommagement isotrope, la surface effective est alors définie par $\tilde{S} = S - S^D = S(1 - D)$. En utilisant cette surface, la contrainte effective (i.e. non-endommagée) est alors définie par

$$ \tilde{\sigma} = \frac{\sigma}{1 - D}. \qquad (32) $$

Cette contrainte effective sera alors introduites dans la surface de charge pour l'implémentation du modèle constitutif dans un code éléments finis. Dans le modèle proposé, l'évolution de la variable d'endommagement est choisie tel que

$$ \dot{D} = \frac{\dot{\lambda}}{g(1 - D)k_c} \sqrt{3J_2(\boldsymbol{S}) + \frac{2}{9}(\alpha^+ \langle p \rangle - \alpha^- \langle -p \rangle)^2} \exp\left(-\frac{\kappa}{k_c}\right) = \exp\left(-\frac{\kappa}{k_c}\right) \frac{\sqrt{\frac{2}{3}\boldsymbol{n}:\boldsymbol{n}}}{k_c} \dot{\lambda}, \qquad (33) $$

où $k_c$ est un paramètre matériau. La nucléation, croissance et coalescence de cavités induites dans la déformation des polymères semi-cristallins chargés provoquent un adoucissement des contraintes. La décroissance des contraintes pour un accroissement des déformations conduit au problème numérique bien connu de localisation des déformations dans une infime zone de la structure. La réponse numérique provenant d'une simulation éléments finis utilisant la plasticité continue classique d'un matériau poséant de l'adoucissement présente une dépendance au maillage non désirée. Dans le but d'outrepasser cette limitation due à la nature du matériau, une formulation non-locale où la variable d'endommagement est définie par un pondérage spatial, est utilisé dans le modèle constitutif.

Une variable non-locale $\tilde{a}$ d'un point matériel $\tilde{x}$ est mathématiquement définie comme un pondérage spatial des variables locales $a$ dans tout les points matériels d'un corps $\mathcal{B}$ tel que

$$ \tilde{a}(\tilde{x}) = \frac{1}{\Psi(\tilde{x})} \int_{\mathcal{B}} \Psi(\tilde{x}; \tilde{y}) a(\tilde{y}) d\mathcal{B}(\tilde{y}). \qquad (34) $$
$\Psi(\vec{x};\vec{y})a(\vec{y})$ est une fonction de poids de type Gaussienne donnée par

$$ \Psi(\vec{x}; \vec{y}) a(\vec{y}) = \frac{1}{c} \exp\left(-\frac{r^2}{2l^2}\right), \qquad (35) $$

où le scalaire $c$ dépend de la dimension du problème. $c$ est définie par $\sqrt{2\pi l}$ pour les problèmes unidimensionnels, $2\pi l^2$ pour le bidimensionnel et $\sqrt{\frac{\pi^3}{2}l^3}$ pour le cas tridimensionnel. La fonction de poids dépend donc uniquement de la distance $r = ||\vec{x} - \vec{y}||$. Le paramètre $l$ détermine la taille du volume qui influe sur le calcul de la quantité non-locale. $\Psi(\vec{x})$ est un facteur de normalisation défini par

$$ \Psi(\vec{x}) = \int_{\mathcal{B}} \Psi(\vec{y}; \vec{x}) d\mathcal{B}. \qquad (36) $$

L'implémentation numérique d'une formulation non-locale nécessite d'avoir accès aux informations de tous les points de gauss de la structure dans le but de satisfaire les conditions de consistances du modèle constitutif à tous les points matériels en même temps. Malheureusement, les routines utilisateurs provenant des codes éléments finis commerciaux ne permettent pas d'avoir accès à toutes les informations de la structure en même temps. De plus, de par la nature des équations non-locales, le calcul du module tangent consistent devient très complexe. Une formulation alternative proposée par Tvergaard et al. [14] est utilisée dans le modèle constitutif afin de surmonter ces difficultés. Les précédents auteurs ont proposé une version non-locale du modèle de Gurson où la variable représentant l'évolution de la fraction volumique des vides $\dot{f}$ est remplacée par $K\dot{f}$ avec $K = \dot{f}/f$. Le facteur de pénalité $K$ est donc le ratio de la quantité non-locale et locale du taux de la fraction volumique des vides calculé au pas de temps précédent. En suivant ce concept, le modèle non-local utilisé dans le modèle constitutif est défini par

$$ \dot{D} = K \exp\left(-\frac{\kappa}{k_c}\right) \frac{\sqrt{\frac{2}{3}\mathbf{n}:\mathbf{n}}}{k_c} \dot{\lambda}, \qquad (37) $$

où le facteur de pénalité $K$ est défini tel que

$$ K = \frac{\dot{D}_n}{\dot{D}_n}. \qquad (38) $$

$\dot{D}_n$ est le taux d'endommagement non-local calculé avec l'opérateur non-local (Eq. (34)) à partir des taux d'endommagements locaux environnants provenant du pas de temps précédent. Avec cette formulation, seul le scalaire $K$ est introduit dans les équations constitutives et la résolution locale de ces équations reste inchangée.

Le modèle constitutif est implémenté en FORTRAN 90 dans une routine utilisateur pour le code élément finis implicite ABAQUS/Standard®. Un schéma implicite est utilisé pour la mise à jour des contraintes et la résolutions des équations constitutives est accomplie en utilisant la méthode de la prédiction viscoélastique / retour viscoplastique. Afin d'implémenter la procédure non-locale aisément, une régularisation de type Lagrangienne est utilisée dans le modèle constitutif ce qui signifie, d'un point de vue numérique, que l'opérateur non-local (Eq. (34)) est calculé seulement une fois au début du calcul et n'est pas mis à jour au cours du temps. Une fois cet opérateur défini, la variable non-locale nécessaire au calcul de $K$ est évaluée pour le point de gauss courant et le facteur de pénalité non-local est introduit pour la résolution des équations constitutives du modèle. Les détails de l'implémentation numérique du modèle constitutif et de la procédure d'identification des paramètres matériaux sont donnés par [15].

# 5 Résultats

Dans le but de comparer les réponses numériques avec les essais expérimentaux, les essais de tractions réalisés à différentes vitesses de sollicitations pour l'identification des paramètres sont simulés avec le modèle constitutif. La géométrie des éprouvettes est donnée par la norme AM3 510. Le maillage éléments finis et les conditions limites sont montrés dans la figure 1(a). Pour les simulations, l'éprouvette est
maillée avec des éléments coques à intégration réduite avec trois points d'intégrations dans l'épaisseur. Les comparaisons entre les réponses données par le modèle constitutif et les essais expérimentaux sont montrées dans les figures 1(b-c). Une bonne corrélation entre le modèle numérique et les résultats expérimentaux est obtenue. Les réponses en terme d'efforts-déplacements montrent que le modèle constitutif est en accord avec les essais expérimentaux aussi bien sur les réponses viscoplastiques (figure 1(b) ) que viscoélastiques (figure 1(c) ). Dans le but d'évaluer la dépendance au maillage du modèle constitutif,

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Maillage éléments finis et conditions limites pour les simulations des essais de tractions uni-axiales (a), réponses efforts-déplacements obtenues avec le modèle constitutif pour différentes vitesses de sollicitations (b) et (c).</figcaption></figure>

une analyse de la striction d'une barre cylindrique est réalisée avec différents maillages. La barre cylindrique de longueur 53.334 mm et de rayon 12.826 mm est soumise à une traction uniaxiale jusqu'à une elongation de 8 mm (avec une vitesse constante de 1mm/s). Afin de faciliter la striction, une imperfection géométrique de 1.8 % est présente au centre de la barre. Pour des raisons d'économie de temps de calculs, les simulations sont réalisées sur un huitième de la barre cylindrique avec les conditions limites appropriées. Les simulations éléments finis sont réalisées avec trois différents maillages (663, 2100 et 6300 éléments hexaédriques à intégration réduite) avec le modèle d'endommagement local ($l = 0$ mm) et non-local ($l = 1$ mm). Les géométries déformées des six simulations sont montrées dans les figures 2(b-c). La localisation des déformations non-désirée dans les éléments situés au centre de la barre apparaît pour les simulations réalisées avec le modèle d'endommagement local. Par contre, pour les simulations effectuées avec le modèle d'endommagement non-local, les déformations des éléments situés au centre de la barre ne sont pas surestimées par rapport aux éléments environnants. Les six efforts de réactions résultants des simulations sont montrés dans la figure 2(a). Pour les réponses obtenues avec le modèle d'endommagement local, une forte influence au maille est observée. Ce phénomène n'est plus observé avec les simulations réalisées avec le modèle d'endommagement non-local. Ces résultats montrent que le problème de dépendance au maillage qui est du à l'adoucissement des contraintes observé dans les polymères semi-cristallins chargés est résolu en utilisant la formulation non-locale proposée.

# 6 Conclusions

Un modèle phénoménologique viscoélastique-viscoplastique couplé à un modèle d'endommagement non-local en grandes déformations est développé dans ce papier afin de simuler par éléments finis le comportement observé expérimentalement d'un polymère semi-cristallin chargé. Le modèle constitutif est capable de représenter la dépendance à la vitesse de déformation et à la pression hydrostatique. Le développement d'un potentiel viscoplastique particulier permet aussi de représenter la variation de volume qui évolue différemment et traction et compression (dilatation et rétractation). Afin d'être consistent entre les mesures expérimentales et le modèle proposé, la formulation hypoélastique basé sur les taux de rotation logarithmique qui réalise le lien entre les contraintes de Cauchy et les déformations de Henky est utilisée dans le modèle constitutif. De plus, l'introduction d'un modèle non-local d'endommagement
<figure><img src="image_3.png" /><figcaption>Fig. 2 – Barre cylindrique sollicitée en traction avec 3 maillages différents : réponses efforts-déplacements des 6 simulations (a), profils déformés des simulations réalisées avec l'endommagement local (b) et non-local (c).</figcaption></figure>

ment permet résoudre la pathologie bien connue de dépendance au maillage dû à l'adoucissement des contraintes présent dans le matériau. Les réponses numériques issues des simulations réalisées avec le modèle proposé sont en accord avec les mesures expérimentales.

Références

[1] Arruda, E.M., Boyce, M.C. *A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials*. Journal of the Mechanics and Physics of Solids 41, 389-412, 1993.

[2] Krempl, E., McMahon, J., Yao, D. *Viscoplasticity based on overstress with a differential growth law for the equilibrium stress*. Mechanics of Materials 5, 35-48, 1986.

[3] Colak, O.U. *Modeling deformation behavior of polymers with viscoplasticity theory based on overstress*. International Journal of Plasticity 21, 145-160, 2005.

[4] Lazzeri, A., Bucknall, C.B. *Applications of a dilatational yielding model to rubber-toughened polymers*, Po-lymer 36, 2895-2902, 1995

[5] R. Balieu, F. Lauro, B. Bennani. *Modèles de comportement pour matériaux polymères soumis au crash*, 10ème colloque national en calcul des structures, hal.archives-ouvertes.fr, 2011.

[6] Jirásek, M. and Rolshoven, S. *Comparison of integral-type nonlocal plasticity models for strain-softening materials*. International Journal of Engineering Science, 41, 1553–1602, 2003.

[7] Truesdell, C., Noll, W. *The Non-Linear Field Theories of Mechanics*. Springer, 3ème édition, 1965.

[8] Xiao, H., Bruhns, O.T., Meyers, A. *Logarithmic strain, logarithmic spin and logarithmic rate*. Acta Mechanica 124, 89-105, 1997.

[9] Hill, R. *Aspects of invariance in solid mechanics*. Advances in Applied Mechanics 18, 1-75, 1978.

[10] R. Raghava, R.M. Caddell, G.S.Y. Yeh. *The macroscopic yield behaviour of polymers*, Journal of Materials Science, Springerlink, 225-232, 1973.

[11] Perzyna, P. *Fundamental problems in viscoplasticity*, Advances in Applied Mechanics 9, 243-377, 1966.

[12] Croix, P., Lauro, F., Oudin, J., Christlein, J. *Improvement of damage prediction by anisotropy of microvoids*. Journal of Materials Processing Technology, 143-144, 202-208, 2003.

[13] Kachanov, L *Time of the rupture process under creep conditions*, Izv Akad Nauk USSR Otd Tekh Nauk 8, 26-31, 1958.

[14] V. Tvergaard, A. Needleman. *Effects of nonlocal damage in porous plastic solids*, International Journal of Solids and Structures, Elsevier, 1063-1077, 1995.

[15] R. Balieu. *A fully coupled viscoelastic-viscoplastic damage model for semi-cristalline polymers*. Ph.D thesis, Université de Valenciennes, France, 2012.