# Modèle d'homogénéisation élasto-plastique pour des structures composites

Minh Hoang, Pham Tuong Minh Duong, Boussad Abbès, Ying-Qiao Guo

► To cite this version:

Minh Hoang, Pham Tuong Minh Duong, Boussad Abbès, Ying-Qiao Guo. Modèle d'homogénéisation élasto-plastique pour des structures composites. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717065

HAL Id: hal-01717065

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

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 d'homogénéisation élasto-plastique pour des structures composites

Minh Phuc HOANG<fnref n="1" />, Pham Tuong Minh DUONG<fnref n="2" />, Boussad ABBÈS<fnref n="3" />, Ying-Qiao GUO<fnref n="4" /> *

**Résumé** — Cet article présente un modèle d'homogénéisation élasto-plastique pour des structures composites de type carton ondulé. Une plaque composée de deux couches planes et d'une couche ondulée est modélisée par une plaque orthotrope homogène. L'homogénéisation est réalisée par une intégration locale à travers l'épaisseur dans chaque couche, suivie d'une cumulation des contributions des trois couches. Un modèle élasto-plastique en contraintes planes est adopté pour modéliser le comportement des papiers cartons. Notre H-modèle est implémenté dans le code de calcul Abaqus à l'aide du sous-programme utilisateur Ugens. Les résultats obtenus par le H-modèle sont comparés à ceux obtenus par les simulations 3D du carton ondulé avec des éléments de coque. La comparaison montre l'efficacité et la précision de notre modèle d'homogénéisation élasto-plastique.

**Mots clés** — homogénéisation, élasto-plasticité, structures composites, carton ondulé, modélisation

## 1. Introduction

Les panneaux en carton ondulé sont des structures sandwich constituées de deux feuilles planes en papier carton écartées par une feuille ondulée appelée cannelure. La simulation de ce type de structures nécessite un maillage très fin pour bien décrire la forme ondulée, ce qui conduit à des temps de calcul importants. Dans nos travaux antérieurs, nous avons proposé un modèle d'homogénéisation du carton ondulé en élasticité [1-2] et nous l'avons appliqué à l'étude des comportements statique et dynamique d'emballages [3]. A notre connaissance, il n'existe pas de modèle d'homogénéisation en plasticité dans la littérature. L'objectif de ce travail est de présenter un modèle d'homogénéisation élasto-plastique orthotrope pour la simulation des structures sandwich de type carton ondulé.

Dans la première partie, nous présentons le comportement mécanique du papier. Les courbes de contrainte-déformation sont présentées pour donner un aperçu de leur comportement. Nous décrivons ensuite le modèle de comportement adopté (IPE-Isotropic Plasticity Equivalent model). Ce modèle combine l'élasticité orthotrope et la plasticité orthotrope. Des essais de traction sont effectués sur chacune des feuilles constituant le carton ondulé dans le but d'obtenir les paramètres de ces matériaux. Dans la deuxième partie, nous présentons notre modèle d'homogénéisation élasto-plastique pour les panneaux sandwich en carton ondulé : volume élémentaire représentatif (VER), repère local de la cannelure, intégration suivant les épaisseurs. Les résultats obtenus par notre H-modèle sont comparés à ceux donnés par les simulations 3D de la structure composite à l'aide du logiciel Abaqus.

## 2. Comportement élasto-plastique du papier carton

### 2.1. Modèle du comportement élasto-plastique anisotrope du papier

Le papier carton est principalement constitué de fibres de bois ayant trois composants : la cellulose, l'hémicellulose et la lignine. À cause du processus de fabrication en continu du papier, les fibres de

<footnote n="1">1 Université de Reims Champagne-Ardenne, GRESPI/MPSE, hmphucvn@gmail.com</footnote>

<footnote n="2">2 Université de Reims Champagne-Ardenne, GRESPI/MPSE, tuongminh80@gmail.com</footnote>

<footnote n="3">3 Université de Reims Champagne-Ardenne, GRESPI/MPSE, boussad.abbes@univ-reims.fr</footnote>

<footnote n="4">4 Université de Reims Champagne-Ardenne, GRESPI/MPSE, yq.guo@univ-reims.fr</footnote>

<footnote>* Auteur correspondant</footnote>
bois sont orientées davantage dans le sens machine. Trois principaux axes sont distingués dans un papier : direction machine (MD), direction transversale (CD) et direction d'épaisseur (ZD). Le papier est un matériau anisotrope : il est 2 à 4 fois plus rigide dans la direction MD que CD, tandis qu'il peut subir une plus grande déformation en CD. Le papier est un matériau hygroscopique, les changements de l'humidité relative et de la température ont beaucoup d'influence sur ses propriétés mécaniques.

Les propriétés dans le plan d'un papier sont relativement faciles à déterminer par des essais de traction. A cause de sa faible épaisseur, les propriétés hors plan sont plus difficiles à obtenir. Stenberg [4] a montré que le module d'Young suivant ZD est environ 200 fois plus faible que celui MD. Stenberg et al. [5] ont observé que la déformation dans le plan est négligeable lors de la compression suivant l'épaisseur. Les coefficients de Poisson $v_{xz}$ et $v_{yz}$ par conséquent sont proches de zéro.

Xia et al. [6] ont proposé un modèle élasto-plastique anisotrope pour le papier et le papier carton. Ce modèle 3D est formulé pour des déformations finies avec un critère de plasticité non quadratique multi-surface, un écrouissage anisotrope, un comportement différent en traction et en compression et un comportement élastique non linéaire. Mäkelä et Östlund [7] ont proposé un modèle de comportement élasto-plastique orthotrope en contraintes planes. Ce modèle simple donne une description précise du comportement des papiers en diverses qualités, n'ayant que quelques paramètres qui sont facilement déterminés à partir des essais de traction simples.

Dans cette étude nous utilisons le modèle élasto-plastique en contraintes planes (IPE) [7] pour modéliser le comportement des constituants du carton ondulé. Le comportement d'élasticité orthotrope en contraintes planes est défini par :

$$ \{\sigma\} = \begin{cases} \sigma_x \\ \sigma_y \\ \sigma_{xy} \end{cases} = [H]\{\epsilon^e\} = \frac{1}{(1-\nu_{xy}\nu_{yx})} \begin{bmatrix} E_x & \nu_{yx}E_x & 0 \\ \nu_{xy}E_y & E_y & 0 \\ 0 & 0 & G_{xy}(1-\nu_{xy}\nu_{yx}) \end{bmatrix} \begin{bmatrix} \epsilon_x^e \\ \epsilon_y^e \\ \gamma_{xy}^e \end{bmatrix} \quad (1) $$

Le vecteur des contraintes déviatoriques et le critère de plasticité du modèle IPE sont donnés par:

$$ \{s\} = \begin{cases} s_x \\ s_y \\ s_z \\ s_{xy} \end{cases} = [L]\{\sigma\} = \frac{1}{3} \begin{bmatrix} 2A & C-A-B & 0 \\ C-A-B & 2B & 0 \\ B-C-A & A-B-C & 0 \\ 0 & 0 & 3D \end{bmatrix} \begin{bmatrix} \sigma_x \\ \sigma_y \\ \sigma_{xy} \end{bmatrix} \quad (2) $$

$$ f = \sigma_{eq} - Y = \left( \frac{3}{2} \langle s \rangle \{s\} \right)^{1/2} - E_0 (\varepsilon_0 + \varepsilon_{eq}^p)^{1/n} = 0 \quad (3) $$

où Y est le seuil de plasticité, A, B, C, D, $E_0$, $\varepsilon_0$, $\varepsilon_{eq}^p$, n sont des paramètres du modèle IPE qui peuvent être déterminés à l'aide d'essais expérimentaux. Nous avons implémenté ce modèle dans le logiciel ABAQUS à l'aide de sous-programme utilisateur UMAT.

## 2.2. Identification inverse du comportement élasto-plastique anisotrope du papier

Les géométries et les dimensions des papiers constituant un carton ondulé sont données dans la Fig. 1. Des éprouvettes standards ont été réalisées à l'aide d'une table à découpe (ZÜND M-1600). Afin d'assurer une meilleure prise au niveau des mors lors du serrage des éprouvettes en papier carton, nous avons collé des morceaux de carton compact rigide aux deux extrémités des éprouvettes (Fig. 2).

<figure><img src="image_2.png" /><figcaption>Fig.1. Géométrie et dimensions du carton ondulé.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Géométrie de l'éprouvette en traction.</figcaption></figure>

Les essais ont été réalisés pour les trois directions caractéristiques (MD, CD et 45°), à une vitesse de 10 mm/min et dans les conditions standards de 23°C et 50% d'humidité. Pour chacune des directions, une dizaine d'essais ont été effectués pour assurer la répétitivité des résultats. Nous avons développé une procédure d'identification inverse en couplant le code de calcul Abaqus et le logiciel d'identification SiDoLo pour identifier les paramètres matériaux des papiers. Fig. 3 montre les courbes expérimentales et celles numériques obtenues par identification du modèle IPE. Le tableau 1 résume les paramètres identifiés pour les deux peaux et la cannelure.

<table><caption>Tableau 1. Paramètres du modèle IPE identifiés pour les peaux et la cannelure.</caption><thead><tr><th colspan="11">Peaux (E<sub>x</sub>, E<sub>y</sub>, G<sub>xy</sub> et E<sub>0</sub> en MPa)</th></tr><tr><th>E<sub>x</sub></th><th>E<sub>y</sub></th><th>v<sub>xy</sub></th><th>G<sub>xy</sub></th><th>E<sub>0</sub></th><th>n</th><th>A</th><th>B</th><th>C</th><th>D</th><th>ε<sub>0</sub></th></tr></thead><tbody><tr><td>2433.2</td><td>859.91</td><td>0.0829</td><td>1077.2</td><td>96.45</td><td>4.97</td><td>1.0</td><td>2.498</td><td>2.498</td><td>1.622</td><td>0.486e-3</td></tr></tbody></table><table><thead><tr><th colspan="11">Cannelure (E<sub>x</sub>, E<sub>y</sub>, G<sub>xy</sub> et E<sub>0</sub> en MPa)</th></tr><tr><th>E<sub>x</sub></th><th>E<sub>y</sub></th><th>v<sub>xy</sub></th><th>G<sub>xy</sub></th><th>E<sub>0</sub></th><th>n</th><th>A</th><th>B</th><th>C</th><th>D</th><th>ε<sub>0</sub></th></tr></thead><tbody><tr><td>1130.4</td><td>625.85</td><td>0.0717</td><td>303.05</td><td>87.31</td><td>4.247</td><td>1.0</td><td>2.178</td><td>2.178</td><td>1.871</td><td>0.923e-3</td></tr></tbody></table>

### 3. Modèle d'homogénéisation élasto-plastique du carton ondulé

L'homogénéisation consiste à représenter un panneau sandwich par une plaque homogène. Nous considérons un VER dans Fig. 4. Les déformations planes peuvent être décomposées en parties de membrane et de flexion :

$$ \{\epsilon\} = \{\epsilon_m\} + z\{\kappa\} \qquad (4) $$

où $\{\epsilon_m\}$ est le vecteur de déformations de membrane et $\{\kappa\}$ est le vecteur de courbures.

Le vecteur de contraintes planes est défini par l'équation (1). Les efforts de membrane, les moments de flexion et de torsion sont obtenus par une intégration des contraintes suivant l'épaisseur du carton ondulé :

<figure><img src="image_4.png" /><figcaption>Fig. 3. Courbes de tractions expérimentales et identifiées pour les peaux (gauche) et la cannelure (droite).</figcaption></figure>
$$ \{N\} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \begin{bmatrix} \sigma_x \\ \sigma_y \\ \sigma_{xy} \end{bmatrix} dz \quad ; \quad \{M\} = \int_{-\frac{h}{2}}^{\frac{h}{2}} z \begin{bmatrix} \sigma_x \\ \sigma_y \\ \sigma_{xy} \end{bmatrix} dz \qquad (5) $$

Pour chaque strate du carton ondulé (peaux et cannelure), nous utilisons 3 points d'intégration suivant son épaisseur. L'algorithme de plasticité est appliqué en chaque point d'intégration pour déterminer l'état de contraintes et la matrice d'élasto-plasticité dans le repère local de chaque strate. Il faut noter que le repère local de la cannelure est défini à l'aide de l'angle $\theta(x)$ (Fig. 4) :

$$ \theta(x) = \tan^{-1} \left( \frac{dh(x)}{dx} \right) \quad ; \quad h(x) = \left( \frac{h_c}{2} - \frac{e_2}{2} \right) \sin \left( 2\pi \frac{x}{P} \right) \qquad (6) $$

Le vecteur des contraintes et la matrice d'élasto-plasticité de la cannelure sont définis dans le repère local (x', z'), et puis transformés dans le repère global (MD-x, CD-y, ZD-z) pour effectuer l'intégration dans l'épaisseur de la plaque carton. D'autre part, la cannelure ne peut quasiment pas supporter la contrainte normale suivant x relative à l'effort normal $N_x$ et au moment de flexion $M_x$ sur une section MD, ceci revient à annuler la contrainte normale suivant la direction x ainsi que la première ligne et la première colonne de la matrice d'élasto-plasticité (Eq. 1).

<figure><img src="image_5.png" /><figcaption>Fig. 4. VER du carton ondulé et points d'intégration numérique suivant l'épaisseur.</figcaption></figure>

Finalement, les efforts de membrane, les moments de flexion et de torsion homogénéisés pour le carton ondulé sont donnés par les expressions suivantes :

$$ \left\{ \begin{aligned} N_x &= \frac{e_l}{2} \sum_{i=1}^{3} \sigma_{xi}^{(1)} \omega_i + \frac{e_3}{2} \sum_{i=1}^{3} \sigma_{xi}^{(3)} \omega_i \\ N_y &= \frac{e_l}{2} \sum_{i=1}^{3} \sigma_{yi}^{(1)} \omega_i + \frac{e_2}{2 \cos(\theta)} \sum_{i=1}^{3} \sigma_{yi}^{(2)} \omega_i + \frac{e_3}{2} \sum_{i=1}^{3} \sigma_{yi}^{(3)} \omega_i \\ N_{xy} &= \frac{e_l}{2} \sum_{i=1}^{3} \sigma_{xyi}^{(1)} \omega_i + \frac{e_2}{2 \cos(\theta)} \sum_{i=1}^{3} \sigma_{xyi}^{(2)} \omega_i + \frac{e_3}{2} \sum_{i=1}^{3} \sigma_{xyi}^{(3)} \omega_i \end{aligned} \right. \qquad (7) $$

$$ \left\{ \begin{aligned} M_x &= \frac{e_l}{2} \sum_{i=1}^{3} \sigma_{xi}^{(1)} \omega_i + \frac{e_3}{2} \sum_{i=1}^{3} \sigma_{xi}^{(3)} \omega_i \\ M_y &= \frac{e_l}{2} \sum_{i=1}^{3} \sigma_{yi}^{(1)} z_i \omega_i + \frac{e_2}{2 \cos(\theta)} \sum_{i=1}^{3} \sigma_{yi}^{(2)} z_i \omega_i + \frac{e_3}{2} \sum_{i=1}^{3} \sigma_{yi}^{(3)} z_i \omega_i \\ M_{xy} &= \frac{e_l}{2} \sum_{i=1}^{3} \sigma_{xyi}^{(1)} z_i \omega_i + \frac{e_2}{2 \cos(\theta)} \sum_{i=1}^{3} \sigma_{xyi}^{(2)} z_i \omega_i + \frac{e_3}{2} \sum_{i=1}^{3} \sigma_{xyi}^{(3)} z_i \omega_i \end{aligned} \right. \qquad (8) $$
où $e_1$, $e_2$ et $e_3$ représentent les épaisseurs des 3 strates, $\omega_i$ représente le poids d'intégration numérique correspondant au point d'intégration $i$ de la strate considérée. Le modèle d'homogénéisation élasto-plastique a été implémenté dans le code de calcul Abaqus à l'aide de sous-programme utilisateur UGENS.

# 4. Résultats

## 4.1. Validation numérique du H-modèle en traction suivant directions MD, 45° et CD

Pour valider notre H-modèle élasto-plastique, des tests de traction suivant les directions MD, 45° et CD sont simulés avec la structure 3D du carton et une plaque 2D homogénéisée (Fig. 5). Les paramètres géométriques du carton ondulé sont (Fig. 4): $P=8$ mm, $h_c=2$ mm, $e_1=0.208$ mm, $e_2=0.144$ mm, $e_3=0.208$ mm. Les propriétés élasto-plastiques des 3 strates sont celles identifiées dans Tableau 1. Le carton ondulé 3D et la plaque homogénéisée 2D sont maillées avec des éléments de coques rectangulaires à intégration réduite (S4R) de taille 1 mm.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Orientation des éprouvettes suivant MD, 45° et CD pour carton ondulé et plaque homogénéisée.</figcaption></figure>

Fig. 6 présente les courbes de force-allongement pour la traction suivant MD, 45° et CD sur le carton ondulé 3D et la plaque homogénéisée. L'écart maximum entre les deux modèles est inférieur à 2%. Les temps CPU sont donnés dans le tableau 2 où nous remarquons qu'avec notre modèle d'homogénéisation le temps CPU est réduit d'un facteur 2 à 4.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Courbes de force-allongement suivant MD, 45° et CD pour carton 3D et plaque homogénéisée.</figcaption></figure>

<table><caption>Tableau 2. Comparaison des temps CPU (en second) pour le carton 3D et la plaque homogénéisée.</caption><thead><tr><th>Modèle</th><th>MD-3D</th><th>MD-H</th><th>45°-3D</th><th>45°-H</th><th>CD-3D</th><th>CD-H</th></tr></thead><tbody><tr><th>Traction</th><td>342</td><td>108</td><td>636</td><td>281</td><td>435</td><td>91</td></tr><tr><th>Flexion</th><td>276</td><td>74</td><td>413</td><td>214</td><td>449</td><td>76</td></tr></tbody></table>
## 4.2. Validation numérique du H-modèle en flexion sur sections MD, 45° et CD

Des tests de flexion sur les sections MD, 45° et CD sont également simulés sur le carton 3D et la plaque 2D homogénéisée. Chaque éprouvette est encastrée à une extrémité et un déplacement perpendiculaire à la plaque est imposé à l'autre extrémité. Les mêmes paramètres géométriques et matériels que ceux de §4.1 sont adoptés.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Courbes de force-flèche en flexion MD, 45° et CD pour le carton 3D et la plaque homogénéisée.</figcaption></figure>

Fig. 7 présente les courbes de force-flèche pour la flexion en MD, 45° et CD avec le carton 3D et la plaque homogénéisée. L'écart maximum entre les deux modèles est inférieur à 2.6%. Les temps CPU sont donnés dans Tableau 2. Nous constatons que notre H-modèle permet de réduire le temps CPU d'un facteur 2 à 4. Fig. 8 montre les isovaleurs des contraintes de von Mises dans la cannelure 3D et dans la plaque homogénéisée sous un moment de flexion CD. Nous remarquons que la distribution des contraintes est très proche entre les 2 modèles.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Contraintes de von Mises dans plaque homogénéisée et cannelure pour flexion CD.</figcaption></figure>

## 4.3. Validation du H-modèle sur une éprouvette à double entaille en traction

Une éprouvette à double entaille soumise à une traction implique un état de contrainte plus complexe. La géométrie et les dimensions sont données sur Fig. 9. La Fig. 10 montre les isovaleurs des contraintes de von Mises sur la plaque homogénéisée et sur la cannelure 3D. Nous remarquons que les distributions des contraintes sont très proches entre les 2 modèles, même autour des entailles où il y a une forte concentration des contraintes. Les courbes de charge-allongement sont quasiment confondues (Fig. 11). Le temps CPU est de 808 s pour la simulation 3D, tandis que le temps CPU n'est que de 309 s pour la simulation H-2D. On note que la plaque homogénéisée pouvait être maillée plus grossièrement pour réduire davantage le temps CPU, car elle n'a pas de cannelure courbe qui nécessite un maillage très fin.
4.4. Validation numérique du H-modèle sur une plaque trouée en traction

Un panneau en carton ondulé trouvé soumis à une traction implique un état de contrainte complexe. Les dimensions du panneau est de 96x96 mm² et le diamètre du trou est de 24 mm. Fig. 11 montre les isovaleurs des contraintes de von Mises dans la plaque homogénéisée et dans la cannelure 3D. Nous constatons que les distributions des contraintes sont très proches entre les 2 modèles. Les courbes de charge-allongement sont quasiment confondues (Fig. 12).

<figure><img src="image_10.png" /><figcaption>Fig. 9. Géométrie des éprouvettes à double entaille pour la plaque homogénéisée et le carton ondulée.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Contraintes de von Mises dans la plaque homogénéisée et dans la cannelure à double entaille.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Contraintes de von Mises dans la plaque homogénéisée et dans la cannelure de la plaque trouée.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12. Courbes de force-allongement de la plaque trouée en traction pour carton 3D et plaque homogénéisée.</figcaption></figure>

**5. Conclusion**

Un modèle d'homogénéisation élasto-plastique est présenté pour la simulation des plaques composites de type carton ondulé. Ce genre de structures 3D est homogénéisé par une plaque 2D orthotrope homogène. Des papiers carton sont caractérisés expérimentalement pour identifier les paramètres du modèle IPE en élasto-plasticité. Dans le VER du carton ondulé, l'intégration de la plasticité s'effectue sur chaque point d'intégration numérique dans le repère local, les contraintes et les matrices élasto-plastiques obtenues sont transformées dans le repère global (MD, CD, ZD), et puis une intégration à travers l'épaisseur (ZD) est effectuée pour obtenir les efforts résultants. Des exemples ont montré que le présent H-modèle est très efficace en préparation des géométries et en calcul de grandes structures composites. La bonne concordance des résultats entre les simulations des structures composites 3D et celles des plaques 2D homogénéisées a prouvé la performance du présent H-modèle.

Le présent H-modèle élasto-plastique sera validé pour d'autres sollicitations plus complexes : cisaillements dans le plan et hors du plan du carton, torsion, compression suivant l'épaisseur... Une validation expérimentale avec des essais réalisés à l'aide d'un système de mesure de déformations sans contact (Aramis) est également envisagée.

**6. Références bibliographiques**

[1] N. Talbi, A. Batti, R. Ayad, Y.Q. Guo, An analytical homogenization model for finite element modelling of corrugated cardboard, Composite Structures, 88 (2), pp. 280-289, 2009.

[2] B. Abbès, Y.Q. Guo, Analytic homogenization for torsion of orthotropic sandwich plates: application to corrugated cardboard, Composite Structures, 92, pp. 699-706, 2010.

[3] A.D. Hammou, P.T.M. Duong, B. Abbès, M. Makhlouf, Y.Q. Guo, Finite-element simulation with a homogenization model and experimental study of free drop tests of corrugated cardboard packaging, Mechanics & Industry, 13(3), pp. 175-184, 2012.

[4] N. Stenberg, A model for the through-thickness elastic-plastic behaviour of paper, International Journal of Solids and Structures, 40, pp. 7483-7498, 2003.

[5] N. Stenberg, C. Fellers, S. Östlund, Plasticity in the thickness direction of paperboard under combined shear and normal loading, Journal of Engineering Materials and Technology, 123(4), pp. 184-190, 2003.

[6] Q.S. Xia, M.C. Boyce, D.M. Parks. A constitutive model for the anisotropic elastic-plastic deformation of paper and paper board, International Journal of Solids and Structures, 39, 4053-4071, 2002.

[7] P. Mäkelä, S. Östlund. Orthotropic elastic-plastic material model for paper materials, International journal of Solids and Structures, 40, 5599-5620, 2003.