# Analyses de stabilité en milieux multi-phasiques

Marie-Angèle Abellan, René de Borst, Jean-Michel Bergheau

► To cite this version:

Marie-Angèle Abellan, René de Borst, Jean-Michel Bergheau. Analyses de stabilité en milieux multi-phasiques. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812908

HAL Id: hal-01812908

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

Submitted on 11 Jun 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
Analyses de stabilité en
milieux multi-phasiques

Marie-Angèle Abellan* — René de Borst**,*** —
Jean-Michel Bergheau*

*LTDS-ENISE, UMR CNRS 5513, E.N.I.S.E, F-42023 Saint-Etienne
{abellan, bergheau}@enise.fr

** Faculty of Aerospace Engineering, TU Delft, NL-2600 GB Delft
r.deborst@lr.tudelft.nl

*** LaMCoS, UMR CNRS 5514, I.N.S.A. de Lyon, F-69621 Villeurbanne

RÉSUMÉ. Cette contribution propose l'analyse de la stabilité de milieux multi-phasiques. Tout d'abord les équations d'un milieu multi-phasiques sont rappelées dans le cadre tridimensionnel. Elles sont ensuite appliquées au cas d'un milieu multi-phasiques pour lequel les constituants présents sont le solide, le fluide et les ions. Les relations constitutives du solide, du fluide et des ions sont données. Puis l'analyse de la stabilité d'une barre infinie est réalisée en partant d'un état de contrainte homogène ainsi que d'une distribution homogène du fluide et des ions et en appliquant une perturbation harmonique. Enfin la discussion permet de souligner l'importance de telles analyses pour la stabilité des calculs numériques des solides élastiques et inélastiques.

ABSTRACT. A stability analysis is presented of a multi-phase medium. Firstly, the equations of a multi-phase medium are recalled in a general three-dimensional setting. They are specialised for a multi-phase medium with solid, fluid and ions as the constituents. The constitutive assump- tions for the solid, the fluid and the ions are given. Finally, the stability of an infinitely long bar is examined, where, starting from a homogeneous stress state and a homogeneous distribution of fluid and ions, a harmonic perturbation is applied. Conclusions for the stability of numerical calculations are drawn both for elastic and for inelastic solids.

MOTS-CLÉS: Stabilité, endommagement, milieux multi-phasiques, milieux hétérogènes.

KEYWORDS: Stability, damage, multi-phase media, heterogeneous media.
## 1. Introduction

La simulation numérique des milieux continus adoucissants fait apparaître une forte dépendance de la solution vis à vis du maillage [BOR 04]. Cette dépendance est liée au changement de caractère des équations aux dérivées partielles du problème aux limites. Localement ce changement dans les équations provient de la relation constitutive choisie pour décrire le comportement du milieu continu qui dans ce cas rend compte d'un comportement adoucissant. Mathématiquement cette phase est accompagnée de l'apparition d'une discontinuité dans la solution. Les schémas numériques standards essayent de capturer cette discontinuité en discrétisant le milieu continu à l'aide d'éléments les plus petits possibles. Ces résultats numériques ont été discutées dans de nombreux articles dont [BOR 04] pour un solide seul et [SCH 96] pour les milieux poreux biphasiques.

Les analyses de dispersion permettent de déterminer si les ondes se propagent dans des zones où le solide présente un comportement adoucissant et de voir si le problème mathématique est bien posé. De telles analyses ont été réalisées pour un solide seul mais, à la connaissance des auteurs, elles n'existent pas pour des milieux multi-phasiques. La présente contribution propose l'analyse de stabilité pour un milieu multi-phasique (solide, fluide et ions). Ceci est la première étape en vue d'une analyse complète de la propagation des ondes caractéristique d'un milieu adoucissant. Ceci permettrait de répondre à la question : " la dépendance vis à vis du maillage observée dans le cas du solide adoucissant existe-t-elle de manière identique pour les milieux multi-phasiques ?"

## 2. Equations de bilan

Supposons un milieu multi-phasique [JOU 95] sous les hypothèses de petits gradients de déplacements, de petits gradients de concentrations, en l'absence de toutes réactions chimiques entre les divers constituants présents et dans le cas dévolutions isothermes. Sous ces hypothèses le bilan de quantité de mouvement de chacun des constituants s'écrit :

$$ \nabla \cdot \boldsymbol{\sigma}_{\pi} + \hat{\mathbf{p}}_{\pi} + \rho_{\pi} \mathbf{g} = \frac{\partial(\rho_{\pi} \mathbf{v}_{\pi})}{\partial t} + \nabla(\rho_{s} \mathbf{v}_{\pi} \otimes \mathbf{v}_{\pi}) \quad (1) $$

avec $\rho_{\pi}$ la masse volumique apparente du constituant $\pi$, $\mathbf{v}_{\pi}$ sa vitesse absolue et $\boldsymbol{\sigma}_{\pi}$ le tenseur des contraintes de Cauchy du constituant $\pi$. Dans cette contribution $\pi = s, f, i$, avec $s, f$ et $i$ qui font respectivement référence au solide, au fluide et aux ions. De plus, $\hat{\mathbf{p}}_{\pi}$ est la source de quantité de mouvement du constituant $\pi$ en provenance des autres constituants, elle prend en compte de possible interactions entre le solide, le fluide et les ions. Evidemment cette source doit vérifier $\sum_{\pi=s,f,i} \hat{\mathbf{p}}_{\pi} = 0$. La prise en compte de l'équation (1) écrite pour $\pi = s, f, i$ et de la relation sur les sources de quantité de mouvement conduit à :

$$ \sum_{\pi=s,f,i} \left( \nabla \cdot \boldsymbol{\sigma}_{\pi} + \rho_{\pi} \mathbf{g} - \frac{\partial(\rho_{\pi} \mathbf{v}_{\pi})}{\partial t} - \nabla(\rho_{s} \mathbf{v}_{\pi} \otimes \mathbf{v}_{\pi}) \right) = \mathbf{0} \quad (2) $$
Si l'on néglige l'accélération de la pesanteur et les termes convectifs, dont l'ordre de grandeur est très petit en comparaison avec les autres termes dans le cas des écoulements dans les milieux poreux, et sous l'hypothèse que les constituants sont intrinsèquement incompressibles, le bilan de quantité de mouvement du milieu total s'écrit :

$$ \sum_{\pi=s,f,i} \left( \nabla \cdot \boldsymbol{\sigma}_{\pi} - \rho_{\pi} \frac{\partial \mathbf{v}_{\pi}}{\partial t} \right) = \mathbf{0} \qquad (3) $$

Le tenseur des contraintes de Cauchy du fluide en présence des ions est défini par $\boldsymbol{\sigma}_f + \boldsymbol{\sigma}_i = -p \mathbf{I}$ avec $p$ la pression dans le fluide en présence des ions et $\mathbf{I}$ le tenseur identité d'ordre deux. De sorte que le bilan de quantité de mouvement total devient :

$$ \nabla \cdot \boldsymbol{\sigma}_s - \nabla p - \rho_s \frac{\partial \mathbf{v}_s}{\partial t} - \rho_f \frac{\partial \mathbf{v}_f}{\partial t} - \rho_i \frac{\partial \mathbf{v}_i}{\partial t} = \mathbf{0} \qquad (4) $$

Un traitement identique permet d'écrire le bilan de masse de chacun des constituants :

$$ \frac{\partial \rho_{\pi}}{\partial t} + \nabla \cdot (\rho_{\pi} \mathbf{v}_{\pi}) = 0 \qquad (5) $$

Définissons à présent les quotients volumiques $n_\pi$ de chacun des constituants $\pi$, sous la contrainte $\sum_{\pi=s,f,i} n_\pi = 1$. Nous supposerons que le quotient volumique des ions est négligeable par rapport aux quotients volumiques du solide et du fluide. Ainsi il vient $n_s + n_f = 1$. A l'aide de cette identité et sous l'hypothèse de constituants intrinsèquement incompressibles, le bilan de masse du solide et du fluide s'écrit :

$$ \nabla \cdot \mathbf{v}_s + n_f \nabla \cdot (\mathbf{v}_f - \mathbf{v}_s) = 0 \qquad (6) $$

### 3. Cinématique et équations de comportement

Les équations de bilan, i.e. le bilan de quantité de mouvement pour le milieu total eq. (4) et les bilans de masse eq. (5) pour $\pi = i$ et eq. (6), sont complétées par la relation cinématique $\dot{\boldsymbol{\epsilon}}_s = \nabla^s \mathbf{v}_s$ avec $\dot{\boldsymbol{\epsilon}}_s$ le tenseur taux de déformation du solide, l'exposant $s$ indiquant que l'on travaille avec la partie symétrique de l'opérateur gradient. Pour le squelette solide, le comportement sera supposé décrit par une loi contrainte-déformation incrémentale linéaire $\dot{\boldsymbol{\sigma}}_s = \mathbf{D}^{\text{tan}} : \dot{\boldsymbol{\epsilon}}_s$ avec $\mathbf{D}^{\text{tan}}$ le tenseur de rigidité matériel tangent du solide.

Si l'on néglige les couplages possibles entre vitesse et flux de chaleur, l'écoulement du fluide en présence d'ions à travers le milieu multi-phasiques s'exprime à l'aide de la loi de Darcy écrite en fonction du potentiel chimique du fluide $\mu_f$, du potentiel chimique des ions $\mu_i$ et d'un tenseur de perméabilité $\mathbf{K}$ :

$$ n_f(\mathbf{v}_f - \mathbf{v}_s) = -\mathbf{K} \cdot \left( \nabla \mu_f + \frac{n_i}{n_f} \nabla \mu_i \right) \qquad (7) $$

L'équation d'état donnant le potentiel chimique du fluide s'écrit [KEM 98] :

$$ \mu_f = p - \phi + \psi \qquad (8) $$
avec $\psi = \psi(n_s)$ le matric potentiel prenant en compte les interactions fluide-solide (capillarité et adsorption) et $\phi = 2RTc_i$ la pression osmotique où $R$ est la constante universelle des gas parfaits, $T$ est la température absolue et $c_i = \frac{\rho_i}{n_f M_i}$ est la concentration des ions par unité de volume du fluide. Le potentiel chimique des ions est défini par :

$$ \mu_i = \mu_{i0} + \frac{RT}{\bar{V}_i} \ln(c_i) \qquad (9) $$

avec $\mu_{i0}$ le potentiel chimique standard des ions, $\bar{V}_i = \frac{n_i M_i}{\rho_i}$ le volume molaire partiel des ions et $M_i$ la masse molaire des ions. Finalement, la diffusion des ions à travers le fluide présent dans le milieu poreux s'exprime par la loi de Fick au moyen d'un tenseur de diffusion des ions $D_i$

$$ n_i (\mathbf{v}_i - \mathbf{v}_f) = -\mathbf{D}_i \cdot \nabla \mu_i \qquad (10) $$

En prenant en compte eqs (8)–(10) le bilan de quantité de mouvement (4) devient :

$$ \nabla \cdot \boldsymbol{\sigma}_s + n_f \mathbf{K}^{-1} (\mathbf{v}_f - \mathbf{v}_s) + \frac{n_i^2}{n_f} \mathbf{D}_i^{-1} (\mathbf{v}_i - \mathbf{v}_f) + \nabla \psi - \rho_s \frac{\partial \mathbf{v}_s}{\partial t} - \rho_f \frac{\partial \mathbf{v}_f}{\partial t} - \rho_i \frac{\partial \mathbf{v}_i}{\partial t} = \mathbf{0} \qquad (11) $$

et le bilan de masse (5) pour les ions :

$$ -\frac{M_i n_i^2}{RT} \mathbf{D}_i^{-1} (\dot{\mathbf{v}}_i - \dot{\mathbf{v}}_f) + \nabla^2 (\rho_i \dot{\mathbf{u}}_i) = \mathbf{0} \qquad (12) $$

## 4. Analyse de stabilité

Dans la suite, nous étudierons le problème d'une barre sollicitée uniaxiallement à partir d'une distribution initiale homogène des contraintes pour le solide et pour le fluide ainsi que d'une distribution initiale homogène de concentration pour les ions. Pour $v_x \neq 0, v_y = 0, v_z = 0$ et après différenciation par rapport au temps $t$, eq. (11) devient (nous noterons en passant que pour du cisaillement pur, $v_x = 0, v_y \neq 0, v_z = 0$, un problème similaire serait obtenu) :

$$ \frac{\partial \dot{\sigma}_s}{\partial x} + \frac{n_f}{K} (\dot{v}_f - \dot{v}_s) + \frac{n_i^2}{n_f D_i} (\dot{v}_i - \dot{v}_f) - \rho_s \ddot{v}_s - \rho_f \ddot{v}_f - \rho_i \dot{v}_i = 0 \qquad (13) $$

avec $\sigma_s$ la contrainte axiale dans le solide et $K$ le coefficient de perméabilité. Un traitement identique des eqs (6) et (12) conduit à :

$$ \frac{\partial^2 \dot{v}_s}{\partial x^2} + n_f \left( \frac{\partial^2 \dot{v}_f}{\partial x^2} - \frac{\partial^2 \dot{v}_s}{\partial x^2} \right) = 0 \quad \text{et} \quad -\frac{M_i n_i^2}{RT D_i} (\ddot{v}_i - \ddot{v}_f) + \rho_i \frac{\partial^2 \dot{v}_i}{\partial x^2} = 0 \qquad (14) $$

Dans le cas unidimensionnel, la relation contrainte-déformation incrémentale linéaire est remplacée par $\dot{\sigma}_s = E^{\text{tan}} \dot{\epsilon}_s$ avec $E^{\text{tan}}$ le module tangent du solide. Par substitution de la relation cinématique écrite pour une dimension, il vient $\dot{\sigma}_s = E^{\text{tan}} \partial v_s / \partial x$.
Pour analyser la stabilité d'un état d'équilibre, une perturbation harmonique est appliquée respectivement aux vitesses unidimensionnelles du solide, du fluide et des ions. Pour un milieu multi-phasiques unidimensionnel infini et en partant d'un état de déformation homogène, on obtient :

$$ \begin{pmatrix} \delta v_s \\ \delta v_f \\ \delta v_i \end{pmatrix} = \begin{pmatrix} A_s \\ A_f \\ A_i \end{pmatrix} \exp(ikx + \lambda t) \quad (15) $$

où $A_s, A_f, A_i$ sont respectivement les amplitudes des perturbations pour le solide, le fluide et les ions, $k$ est le nombre d'onde et $\lambda$ est la valeur propre. Substitution de eq. (15) dans eqs (13)-(14) en tenant compte de l'expression ci-dessus pour $\dot{\sigma}_s$ et en requérant que le système d'équations homogènes obtenu à partir des amplitudes $A_s, A_f$ et $A_i$ possède une solution non-triviale conduit à l'équation charactéristique :

$$ \lambda^3 + a\lambda^2 + b\lambda + c = 0 \quad (16) $$

avec

$$ a = \frac{\rho_s n_f k^2 \rho_i + K^{-1} n_f \frac{M_i}{RT} D_i^{-1} n_i^2 - n_s k^2 \rho_f \rho_i}{(\rho_s n_f - n_s \rho_f - n_s \rho_i) \frac{M_i}{RT} D_i^{-1} n_i^2} \quad (17a) $$

$$ b = \frac{E^{\tan} n_f \frac{M_i}{RT} n_i^2 D_i^{-1} k^2 + K^{-1} n_f k^2 \rho_i - n_s \frac{n_i^2}{n_f} D_i^{-1} k^2 \rho_i}{(\rho_s n_f - n_s \rho_f - n_s \rho_i) \frac{M_i}{RT} D_i^{-1} n_i^2} \quad (17b) $$

$$ c = \frac{n_f \rho_i E^{\tan} k^4}{(\rho_s n_f - n_s \rho_f - n_s \rho_i) \frac{M_i}{RT} D_i^{-1} n_i^2} \quad (17c) $$

D'après le théorème de stabilité de Routh-Hurwitz, toutes les solutions $\lambda(k)$ auront une partie réelle négative si les conditions suivantes sont remplies simultanément : $a > 0$, $c > 0$ et $ab - c > 0$. Pour les valeurs classiques des paramètres apparaissant dans les relations ci-dessus la condition $a > 0$ est toujours vérifiée. Ceci est également vrai pour $c > 0$ tant que $E^{\tan} > 0$ c'est à dire tant que le comportement du solide n'est pas adoucissant. La troisième condition, $ab - c > 0$, peut se réécrire :

$$ \begin{aligned} & \rho_s n_f^2 \rho_i^2 K^{-1} k^4 + K^{-2} n_f^2 \frac{M_i}{RT} D_i^{-1} n_i^2 \rho_i k^2 + n_s^2 \rho_f \rho_i^2 D_i^{-1} \frac{n_i^2}{n_f} k^4 \\ & + E^{\tan} \frac{M_i}{RT} n_i^2 D_i^{-1} \left( K^{-1} n_f^2 \frac{M_i}{RT} n_i^2 D_i^{-1} + n_s n_f \rho_i^2 k^2 \right) k^2 \\ & - \rho_s \rho_i^2 n_s n_i^2 D_i^{-1} k^4 - K^{-1} \frac{M_i}{RT} n_i^4 D_i^{-2} n_s \rho_i k^2 - n_s \rho_f \rho_i^2 K^{-1} n_f k^4 > 0 \end{aligned} \quad (18) $$

Une étude paramétrique a été effectuée pour étudier le signe de $(ab - c)$ pour les valeurs moyennes suivantes : $0.15[M] \le c_i \le 0.30[M]$, $0.13 \le n_f \le 0.35$, $0 \le n_i \le 0.2 \times 10^{-4}$, $6MPa \le E^{\tan} \le 2500MPa$ et $10^{-22}m^4/Ns \le K \le 5 \times 10^{-21}m^4/Ns$, $M_i = 58.5 \times 10^{-3}kg$, $R = 8.314J/Kmol$, $T = 394K$, $D_i = 3.3 \times 10^{-11}m^2/s$,
$\rho'_s = 1330\text{kg/m}^3$ et $\rho'_f = 1000\text{kg/k}^3$ avec $\rho'_s, \rho'_f$ les densités du solide et du fluide définies par $\rho_\pi = n_\pi \rho'_\pi$ en fonction des masses volumiques apparentes. Pour ce choix de paramètres qui sont des moyennes représentatives pour des tissus biologiques [KEM 98], le signe de $(ab-c)$ est positif.

**5. Discussion**

L'analyse précédente fait apparaître que la condition de stabilité pour un milieu multi-phasiques est la même que pour le solide au moins pour le cas d'une étude unidimensionnelle. En effet, la stabilité est assurée tant que le module tangent du solide reste positif et au contraire une perte de stabilité se produit dès que ce module tangent devient négatif.

Cette étude de stabilité sera complétée par une analyse de dispersion pour le cas d'ondes se propageant dans un milieu multi-phasiques. Ceci permettra de clarifier si le changement de forme des équations du problémes aux limites associé, c'est à dire la perte d'hyperbolicité, se produit également au même moment pour un milieu multi-phasiques et pour le solide associé. Si ceci était le cas, la prise en compte de l'écoulement des fluides en plus du comportement purement mécanique n'aurait aucun effet régularisant. Si au contraire, des différences existes, les milieux multi-phasiques aparaient comme moins sensibles à la perte d'hyperbolicité et ainsi les phénomènes de dépendance de la solution vis à vis du maillage, qui handicapent les analyses numériques des matériaux adoucissants, pourraient être éliminée.

**6. Bibliographie**

[BOR 04] DE BORST R., « Damage, material instabilities, and failure », *Encyclopedia of Computational Mechanics*, Stein E., de Borst R. and Hughes T.J.R. (eds), Wiley, Chichester, 2004, vol. 2, ch. 10.

[SCH 96] SCHREFLER B.A., SANAVIA L., MAJORANA C.E., « A multiphase medium model for localisation and post-localisation simulation in geomaterials », *Mechanics of Cohesive-frictional Materials*, vol. 1, 1996, p. 95-114.

[JOU 95] JOUANNA P., ABELLAN M.-A., « Generalized approach to heterogeneous media », *Modern Issues in Non-Saturated Soils*, Gens A., Jouanna P. and Schrefler B.A. (eds), Springer-Verlag, Wien - New York, 1995, p. 1-128.

[KEM 98] VAN KEMENADE P.M., « Water and ion transport through intact and damaged skin », Dissertation, Eindhoven University of Technology, 1998.