# Étude théorique et numérique de tissus biologiques en hyperélasticité anisotrope

François Peyraut, Zhi-Qiang Feng, Nadia Labed, Christine Renaud

► To cite this version:

François Peyraut, Zhi-Qiang Feng, Nadia Labed, Christine Renaud. Étude théorique et numérique de tissus biologiques en hyperélasticité anisotrope. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422250

HAL Id: hal-01422250

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

Submitted on 24 Dec 2016

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 CC BY 4.0 - Attribution - International License
# Etude théorique et numérique de tissus biologiques en hyperélasticité anisotrope

F. Peyraut¹, Z.-Q. Feng², N. Labed¹, C. Renaud²

¹ Université de Technologie de Belfort-Montbéliard
Laboratoire M3M, 90010 Belfort cedex, France
{francois.peyraut,nadia.labed}@utbm.fr

² Université d'Evry-Val d'Essonne
Laboratoire LME-Evry, 40 rue du Pelvoux, 91020 Evry, France
{feng,renaud}@iup.univ-evry.fr

**Résumé** – Pour déterminer les déformations et les contraintes au sein de tissus biologiques tels que les ligaments, les tendons ou les parois artérielles, les lois de comportements hyperélastiques anisotropes sont souvent utilisées dans le cadre de la méthode des éléments finis [1]. Dans cet article, on se propose de réaliser une telle étude en parallèle avec une analyse théorique portant sur l'angle entre les fibres de collagène et la direction circonférentielle. En particulier, l'importance d'un angle critique égal à 54.73° est discutée.

**Mots clés** – Biomécanique, hyperélasticité anisotrope, modèle HGO, éléments finis.

## 1 Introduction

Pour déterminer les déformations et les contraintes au sein de tissus biologiques tels que les ligaments, les tendons ou les parois artérielles, les lois de comportements hyperélastiques anisotropes sont souvent utilisées dans le cadre de la méthode des éléments finis [1]. Dans cet article, on se propose de réaliser une telle étude en parallèle avec une analyse analytique. Cette analyse complémentaire permet de comprendre pourquoi la correspondance n'est pas biunivoque entre la déformation principale $\lambda_2$ et le quatrième invariant de la matrice de dilatation pour un modèle usuel tel que celui proposé par Holzapfel, Gasser et Ogden [4-5]. On établit en effet qu'une correspondance non bijective apparaît lorsque l'angle $\beta$ entre les fibres de collagène et la direction circonférentielle dépasse une valeur critique égale à 54.73°. L'importance de cet angle critique a déjà été relevée dans [7]. Pour déterminer cet angle, on considère l'essai de traction uniaxiale non contraint habituellement utilisé pour identifier les propriétés mécaniques de bandes de tissus artériels [2]. Les déformations associées à cet essai sont homogènes. Grâce à un changement d'inconnu approprié, on montre que les équations d'équilibre se réduisent alors à une équation polynômiale de degré 3 que l'on peut résoudre algébriquement par les formules de Cardan. Les formules analytiques obtenues montrent que deux solutions peuvent exister lorsque l'angle $\beta$ dépasse une valeur seuil égale à 54.73°.

Le modèle HGO a été par ailleurs été implémenté dans le code de calcul aux éléments FER à l'aide d'une formulation lagrangienne totale [9]. Les calculs numériques corroborent les résultats obtenus avec l'étude analytique. On note également des résultats conformes aux observations réalisées dans [5] dans le cas où les deux extrémités de l'échantillon sont contraintes.
## 2 Hyperélasticité anisotrope

Les densités d'énergie les plus utilisées en hyperélasticité anisotrope impliquent des lois puissance [2] ou des lois exponentielles [3-4]. Il est généralement admis que l'anisotropie est due aux fibres de collagène [5] alors que la matrice se comporte de manière isotrope. C'est la raison pour laquelle la densité d'énergie W est séparée en une partie isotrope et une partie anisotrope [1] :

$$W = W_{iso} + \sum_{a=1}^{n} W_{ani}^{a} \qquad (1)$$

Chaque densité anisotrope $W_{ani}^a$ fait référence à une direction privilégiée du matériau. Pour modéliser les ligaments ou les tendons, le nombre de famille de fibres $n$ est égal à 1. Il est égal à 2 pour représenter le comportement des parois artérielles. Dans la loi HGO proposée par Holzapfel, Gasser et Ogden [4-5], deux densités anisotropes sont superposées avec deux directions privilégiées $\mathbf{a}^1$ et $\mathbf{a}^2$ correspondant à deux familles de fibres distinctes :

$$\mathbf{a}^1 = \{\cos(\beta), \sin(\beta), 0\}^T \quad ; \quad \mathbf{a}^2 = \{\cos(\beta), -\sin(\beta), 0\}^T \qquad (2)$$

où $\beta$ représente l'angle entre les fibres de collagène et la direction circonférentielle sur des tissus extraits d'aorte abdominale humaine (Figure 1a). Dans la suite, on notera en abrégé $c$, $s$ et $t$ le cosinus, le sinus et la tangente de l'angle $\beta$. Les lois hyperélastiques anisotropes incluent les trois invariants classiques $I_1$, $I_2$ et $I_3$ de la matrice de dilation $C$ ainsi que deux invariants additionnels $J_4^a$ et $J_5^a$ liés à l'anisotropie des tissus et définis en fonction du tenseur structurel $\mathbf{M}^a$ [6] :

$$I_1=\mathrm{tr}(\mathbf{C}), \quad I_2=\mathrm{tr}(\mathrm{cov}(\mathbf{C})), \quad I_3=\mathrm{det}(\mathbf{C}), \quad J_4^a=\mathrm{tr}(\mathbf{C}\mathbf{M}^a), \quad J_5^a=\mathrm{tr}(\mathbf{C}^2\mathbf{M}^a), \quad \mathbf{M}^a=\mathbf{a} \otimes \mathbf{a} \qquad (3)$$

Classiquement, le second tenseur des contraintes de Piola-Kirchhoff $S$ dérive de la densité $W$. En introduisant la matrice des cofacteurs de $C$, $\mathrm{Cof}(C)=\mathrm{det}(C)C^{-T}$, et en dérivant $W$, on obtient [2] :

$$\mathbf{S} = 2 \left[ \frac{\partial W}{\partial I_1} \mathbf{I} + \frac{\partial W}{\partial I_2} (I_1 \mathbf{I} - \mathbf{C}) + \frac{\partial W}{\partial I_3} \mathrm{cov}(\mathbf{C}) + \frac{\partial W}{\partial J_4^a} \mathbf{M}^a + \frac{\partial W}{\partial J_5^a} (\mathbf{C}\mathbf{M}^a + \mathbf{M}^a\mathbf{C}) \right] \qquad (4)$$

Par ailleurs, le comportement anisotrope du modèle HGO est de type exponentiel :

$$J_4^a \ge 1: W_{ani}^a = \frac{k_1}{2k_2} \left[ \exp \left( k_2 (J_4^a - 1)^2 \right) - 1 \right] \quad ; \quad J_4^a < 1: W_{ani}^a = 0 \qquad (5)$$

La densité $W_{ani}^a$ ne dépend que du quatrième invariant $J_4^a$ qui représente l'allongement dans la direction des fibres. Comme elle est nulle lorsque $J_4^a$ est plus petit que 1, le modèle HGO prend en compte le fait que les fibres ne travaillent pas en compression. La partie isotrope de la densité, décrivant le comportement de la matrice, est représentée par le modèle néo-hookéen :
$$ W_{iso} = c_1(I_1 I_3^{-1/3} - 3) \qquad (6) $$

On reprendra dans cet article les valeurs des paramètres matériels $c_1$, $k_1$ et $k_2$ qui ont été choisis dans [8] de manière à ajuster le modèle avec des données expérimentales : $c_1=10.2069$ kPa ; $k_1=0.0017$ kPa ; $k_2=882.847$. On supposera que ces paramètres restent constant quand $a$ change dans l'expression (5) de $W_{ani}^a$. Cette hypothèse est légitime car les propriétés des fibres sont indépendantes de leur orientation. Dans la suite, on ne fera plus apparaître l'indice $a$.

## 3 Etude théorique

On considère le modèle HGO dans le cadre d'un essai de traction uniaxiale non contraint (Figure 1b). Ce test est habituellement utilisé pour identifier les propriétés mécaniques de bandes de tissus artériels.

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

<figure><img src="image_3.png" /><figcaption>Figure 1 – a) Angle $\beta$; b) Essai de traction uniaxiale non contraint</figcaption></figure>

Pour calculer les déformations principales $\lambda_1$, $\lambda_2$ et $\lambda_3$, on se propose de traiter analytiquement puis numériquement le problème. La contrainte d'incompressibilité $J=\det(F)=1$, F étant la matrice gradient des déformations, est prise en compte dans le modèle éléments finis par un terme de pénalité de la forme $k(J-1)^2$. Ce terme est ajouté à la densité d'énergie (1) avec $k$ égal à $10^5$. D'un point de vue analytique, la contrainte d'incompressibilité est traitée à l'aide d'un multiplicateur de Lagrange $p$ représentant une pression hydrostatique. On rappelle que $S$ est relié aux contraintes de Cauchy $\sigma$ par la relation classique :

$$ \sigma = \frac{1}{\det(F)} F S F^T \qquad (7) $$

En utilisant la condition de bord libre $\sigma_{33}=0$ avec l'expression (4) des contraintes et les densités d'énergie (5) et (6), le multiplicateur de Lagrange disparaît et on obtient :

$$ J_{4 \ge 1}: \begin{cases} \sigma_{11} = 2 \left[ c_1 (\lambda_1^2 - \lambda_1^{-2} \lambda_2^{-2}) + 2k_1 \lambda_1^2 c^2 (J_4-1) \exp(k_2(J_4-1)^2) \right] \\ \sigma_{22} = 2 \left[ c_1 (\lambda_2^2 - \lambda_1^{-2} \lambda_2^{-2}) + 2k_1 \lambda_2^2 s^2 (J_4-1) \exp(k_2(J_4-1)^2) \right] \end{cases} \qquad (8) $$

$$ J_{4 < 1}: \sigma_{11} = 2 c_1 (\lambda_1^2 - \lambda_1^{-2} \lambda_2^{-2}); \sigma_{22} = 2 c_1 (\lambda_2^2 - \lambda_1^{-2} \lambda_2^{-2}) \qquad (9) $$
D'après les formules (8) et (9), les contraintes s'expriment différemment selon la valeur de $J_4$ par rapport à 1. Il est naturel de distinguer deux cas :

Cas 1 : $J_4 \ge 1$

La condition de bord libre $\sigma_{22}=0$ appliquée à la formule (8) donne :

$$c_1(\lambda_2^2 - \lambda_1^{-2}\lambda_2^{-2}) + 2k_1\lambda_2^2 s^2 (J_4 - 1) \exp[k_2(J_4 - 1)^2] = 0 \quad (10)$$

Par ailleurs, la définition (3) du quatrième invariant peut se formuler de la manière suivante :

$$(\lambda_1 c J_4^{-1/2})^2 + (\lambda_2 s J_4^{-1/2})^2 = 1 \quad (11)$$

Cela suggère l'utilisation d'un système de coordonnées polaires avec $\theta$ comme angle polaire :

$$\lambda_1 c = J_4^{1/2} \cos(\theta); \quad \lambda_2 s = J_4^{1/2} \sin(\theta) \quad (12)$$

En reportant (12) dans l'équation non linéaire (10), les déformations principales $\lambda_1$ et $\lambda_2$ sont remplacées par deux nouvelles variables $\theta$ and $J_4$ et on obtient une équation polynômiale cubique :

$$\cos(\theta) - \cos^3(\theta) = A, \quad A = s^2 c J_4^{-3/2} \left\{ 1 + \frac{2k_1}{c_1} s^2 (J_4 - 1) \exp \left[ k_2 (J_4 - 1)^2 \right] \right\}^{-1/2} \quad (13)$$

Les trois racines de cette équation peuvent être calculées analytiquement à l'aide des formules de Cardan. Après discussion sur le signe du discriminant $\Delta=A^2/4-1/27$ et élimination des racines indésirables (celles purement imaginaires par exemple), on obtient une solution unique définie par :

$$\lambda_1 = J_4^{1/2} \frac{\cos(\theta)}{c} ; \quad \lambda_2 = J_4^{1/2} \frac{\sin(\theta)}{s} ; \quad \cos(\theta) = \frac{2\sqrt{3}}{3\cos(\alpha/3)} ; \quad \alpha = \arccos\left(-\frac{3\sqrt{3}A}{2}\right) \quad (14)$$

Cas 2 : $J_4 < 1$

Le cas où $J_4$ est inférieur à 1 est plus simple à traiter car la loi de comportement se réduit au modèle néo-hookéen. Il résulte alors de la condition de bord libre $\sigma_{22}=0$ appliquée à (9) que :

$$\lambda_2 = \lambda_1^{-1/2} \quad (15)$$

Cependant, si on souhaite étudier la dépendance de $\lambda_1$ et $\lambda_2$ par rapport à $J_4$, pour avoir une approche homogène à celle réalisée pour le cas 1 précédent, on est amené à résoudre une nouvelle équation polynômiale cubique. Celle-ci se déduit facilement de (11) et (15) :

$$\lambda_1^3 - J_4 c^{-2} \lambda_1 + t^2 = 0 \quad (16)$$

La discussion portant sur le signe du discriminant $\Delta_1=t^4/4-J_4^3/(27c^6)$ montre que trois cas peuvent se présenter. Soit il n'existe aucune solution, soit il en existe une seule, soit il en existe
deux distinctes. Ces différentes situations se présentent en fonction de la valeur prise par l'angle $\beta$ relativement à un angle critique défini par $\beta_c = \arccos(\sqrt{3}/3) \approx 54.73^\circ$. Si $\beta < \beta_c$, ou si $\beta \ge \beta_c$ et $\Delta_1 > 0$, il n'existe aucune solution. Si $\beta \ge \beta_c$ et $\Delta_1 \le 0$, il en existe une ($\lambda_1$) ou deux ($\lambda_{1a}$ et $\lambda_{1b}$) :

$$ \Delta_1=0 : \quad \lambda_1 = \sqrt{J_4 \cos(\beta_c) / \cos(\beta)} \qquad (17) $$

$$ \Delta_1 < 0 : \quad \lambda_{1a} = 2 \sqrt{\frac{J_4}{3c^2}} \cos\left(\frac{\gamma}{3}\right), \quad \lambda_{1b} = 2 \sqrt{\frac{J_4}{3c^2}} \cos\left(\frac{\gamma + 4\pi}{3}\right), \quad \gamma = \arccos\left[-\frac{3\sqrt{3} s^2 c}{2J_4^{3/2}}\right] \qquad (18) $$

Lorsque la solution $\lambda_1$ existe, $\lambda_2$ est calculé en fonction de $\lambda_1$ par la formule (15). La variation de $\lambda_2$ en fonction de $J_4$ est présentée sur la Figure 2a pour deux angles $\beta$ différents (20° et 70°). Lorsque $\beta$ est strictement inférieur à $\beta_c$ (20°), on ne trouve aucune solution correspondante à $J_4 < 1$. Le tissu biologique ne peut donc être comprimé dans la direction des fibres ($J_4 < 1$) que si $\beta$ est supérieur à $\beta_c$. Dans ce cas (70°), le caractère non biunivoque de la dépendance de $\lambda_2$ par rapport à $J_4$ apparaît clairement. Il existe bien deux solutions distinctes conformément à la prédiction de la formule (18). Enfin, on peut montrer que le point de rebroussement qui apparaît sur la Figure 2a est une racine double de l'équation cubique (16) et que cette racine double est donnée par la formule (17).

<figure><img src="image_4.png" /><figcaption>Figure 2 - a) Courbe $\lambda_2 - J_4$; b) Comparaison numérique-analytique ($\beta=70^\circ$)</figcaption></figure>

# 4 Calcul éléments finis

Le modèle HGO a été implémenté dans le code de calcul aux éléments finis FER à l'aide d'une formulation lagrangienne totale [9]. Les calculs numériques corroborent les résultats obtenus avec l'étude analytique. En effet, on constate sur la Figure 2b une parfaite correspondance entre les déformations principales calculées analytiquement et numériquement.

Un calcul numérique a également réalisé dans le cas où les extrémités de l'échantillon sont contraintes (Figure 3). On note un rétrécissement de la section dans la direction 2 et un gonflement dans la direction 3, ce qui est typique d'un comportement non isotrope. Ce résultat est conforme aux observations réalisées dans [5]. Il est également cohérent avec le calcul analytique présenté sur la Figure 2b qui montre, qu'à partir d'une valeur seuil de $\lambda_1$, $\lambda_3$ augmente alors que $\lambda_2$ diminue.
<figure><img src="image_5.png" /><figcaption>Figure 3 - Test de traction uniaxiale contraint - calcul numérique avec FER</figcaption></figure>

# 5 Conclusions

On a présenté dans cet article une étude théorique ainsi qu'une analyse éléments finis du modèle hyperélastique anisotrope HGO. Les calculs portent sur l'essai de traction uniaxiale habituellement utilisé pour identifier les propriétés mécaniques de bandes de tissus artériels. Moyennant un changement d'inconnue, on a montré que les déformations principales peuvent être calculées analytiquement. Ce calcul permet d'expliquer la correspondance non bijective entre la déformation principale $\lambda_2$ et le quatrième invariant de la matrice de dilatation lorsque l'angle entre les fibres de collagène et la direction circonférentielle dépasse une valeur critique égale à 54.73°.

L'analyse éléments finis a été réalisée à l'aide du code de calcul FER dans lequel nous avons implémenté la loi HGO. Dans le cas où les extrémités de l'échantillon sont contraintes, on trouve des résultats en bon accord avec la littérature.

# Références

[1] Weiss, J.A., Maker, B.N., Govindjee, S. «Finite element implementation of incompressible, transversely isotropic hyperelasticity», *Computer methods in applied mechanics and engineering*, 135, 107-128, 1996.

[2] Schröder, J., Neff, P., Balzani D. «A variational approach for materially stable anisotropic hyperelasticity», *International Journal of Solids and Structures*, 42, 4352–4371, 2005.

[3] Fung, Y.C., Fronek, K., Patitucci, P. «Pseudoelasticity of arteries and the choice of its mathematical expression», *Am. J. Physiol.*, 237, H620-H631, 1979.

[4] Holzapfel, G.A., Gasser, T.C., Ogden, R.W. «A new constitutive framework for arterial wall mechanics and a comparative study of material models», *Journal of Elasticity*, 61, 1-48, 2000.

[5] Gasser, T.C., Ogden, R.W., Holzapfel, G.A. «Hyperelastic modelling of arterial layers with distributed collagen fibre orientations», J. R. Soc. Interface, 3, 15-35, 2006.

[6] Spencer, A.J.M. Isotropic polynomial invariants and tensor functions, In: Boehler, J.P. (Ed.), *Applications of Tensor Functions in Solids Mechanics*, CISM Course No. 282. Springer Verlag, 1987.

[7] Guo, Z.Y., Peng, X.Q., Moran, B. «A composites-based hyperelastic constitutive model for soft tissue with application to the human annulus fibrosus» *J. Mech. Phys. Solids*, 54:1952-1971, 2006.

[8] Balzani, D., Neff, P., Schröder, J., Holzapfel G.A. «A polyconvex framework for soft biological tissues. Adjustment to experimental data», *International Journal of Solids and Structures*, 43, 6052–6070, 2006.

[9] http://lmee.univ-evry.fr/~feng/FerSystem.html