# Sur le volume élémentaire représentatif de la loi de perte de charge d'Ergun

H. J. Herrmann, V. Komiwes, P. Mege, Y. Meimon

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

## To cite this version:

H. J. Herrmann, V. Komiwes, P. Mege, Y. Meimon. Sur le volume élémentaire représentatif de la loi de perte de charge d'Ergun. Oil & Gas Science and Technology - Revue d'IFP Energies nouvelles, 1999, 54 (5), pp.577-585. 10.2516/ogst:1999049. hal-02075838

HAL Id: hal-02075838

https://ifp.hal.science/hal-02075838v1

Submitted on 21 Mar 2019

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Sur le volume élémentaire représentatif de la loi de perte de charge d'Ergun

V. Komiwes¹, P. Mège¹, Y. Meimon¹ et H.J. Herrmann²

1 Institut français du pétrole, 1 et 4, avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex - France

2 École de physique et de chimie industrielle de Paris, 15, rue Vauquelin, 75005 Paris - France

e-mail: philippe.mege@ifp.fr

**Résumé** — Diverses expressions permettent de calculer la perte de charge occasionnée par le passage d'un fluide visqueux incompressible en écoulement à travers un milieu granulaire. Ces lois ont toutes une origine phénoménologique et sont obtenues à partir d'expériences réalisées à une échelle macroscopique. L'une d'entre elles, la loi d'Ergun, est couramment utilisée sous une forme locale dans des modèles de couplage fluide-solide pour réaliser le transfert de quantité de mouvement entre les deux phases sans réelle justification.

Dans cet article, la simulation numérique directe des équations de Navier-Stokes est utilisée pour calculer les pertes de charge de fluide visqueux incompressible en écoulement à travers un milieu granulaire. Des simulations effectuées pour différents volumes de matériau granulaire permettent de retrouver la loi de perte de charge d'Ergun, de définir un volume élémentaire représentatif et de prouver le caractère local de cette loi.

Mots-clés : loi d'Ergun, simulation numérique directe, milieux poreux, perte de charge, volume élémentaire représentatif, écoulements de fluide.

**Abstract — On the Representative Elementary Volume of the Ergun's Law** — Several expressions can be used to predict the pressure drop caused by a viscous incompressible fluid flow through a granular media. These laws are empirical and are obtained at macroscopic scale.

One of these, Ergun's law, is currently used in a local form to modelling the momentum transfer between fluid and solid. There is not a real justification.

In these study, the direct simulation of Navier-Stokes Equation is used to calculate pressure drop caused by viscous incompressible fluid flow through a granular media. Some simulations realised for different volume of granular media allow to find a representative elementary volume and prove the local character of the Ergun's law.

Keywords: Ergun's law, direct numerical simulation, porous media, pressure drop, representative elementary volume, fluid flow.
## NOMENCLATURE

$d_p$ diamètre des particules (m)

$r$ rayon des particules (m)

$\rho_f$ masse volumique du fluide (kg·m⁻³)

$\mu_f$ viscosité du fluide (kg·m⁻¹·s⁻¹)

$V_{SF}$ vitesse superficielle du fluide (m·s⁻¹)

$\vec{v} = v_i \vec{u}_{xi}$ champ de vitesse du fluide (m·s⁻¹)

$p$ pression du fluide (Pa)

$\Delta p$ perte de charge (Pa)

$h$ taille de la maille rapportée au rayon des particules
(-)

$\Delta t_f$ pas de temps adimensionnel de discrétisation
(-)

$\varepsilon$ fraction volumique de solide du lit (-)

$Re$ nombre de Reynolds particulière (-)

$L_y$ hauteur du milieu granulaire (m).

## INTRODUCTION

Diverses expressions permettent de calculer la perte de charge occasionnée par le passage d'un fluide visqueux incompressible en écoulement à travers un milieu granulaire. Une revue de ces lois est proposée par Dullien [1]. Parmi les plus connues, nous avons la relation de Carman-Kozeny et de Rumpf-Gupte dans le cas d'écoulement laminaire, la relation de Burke-Plummer dans le cas d'écoulement turbulent, et la loi d'Ergun [2] qui couvre les écoulements laminaires, turbulents et intermédiaires. Ces lois ont toutes une origine phénoménologique et sont obtenues à partir d'expériences réalisées à une échelle macroscopique.

L'une d'entre elles, la loi d'Ergun, est couramment utilisée sous une forme locale dans des modèles de couplage fluide-solide pour réaliser le transfert de quantité de mouvement entre les deux phases. Citons en particulier :

- les modèles à deux fluides [3-6] où une approche de type eulérien est utilisée pour décrire les deux phases ;

- les modèles utilisant une approche eulérienne pour la description du fluide et une représentation discrète des particules [7, 8] de type sphères molles ou sphères dures.

À ce sujet, aucune justification n'est faite quant au caractère local de cette loi et au volume élémentaire représentatif. Le présent travail propose d'utiliser la simulation numérique directe des équations de Navier-Stokes pour vérifier le caractère local de la loi d'Ergun et définir un volume élémentaire représentatif. Ce type d'approche a déjà été utilisé dans différents travaux [9-11] dans le cadre d'écoulements de Stokes et a permis de vérifier la loi de Carman-Kozeny.

# 1 FONDEMENT ET EXPRESSION DE LA LOI D'ERGUN

Le domaine de validité de cette loi concerne les écoulements d'un fluide à travers des milieux poreux constitués d'un assemblage de particules sphériques homogènes et fixes allant des régimes laminaires aux régimes turbulents.

Cette loi s'obtient par l'addition des lois de Blake-Kozeny et de Burke-Plummer valables respectivement dans le cas limite des écoulements laminaire et turbulent. Le fondement théorique de ces lois provient de l'analogie avec l'expression de la perte de charge dans un tube cylindrique. Cette analogie permet d'obtenir des expressions théoriques corrigées expérimentalement par un facteur multiplicatif proche de 2 en raison de la tortuosité [12] pour donner en définitive les lois respectives de Blake-Kozeny et de Burke-Plummer :

$$ \frac{\Delta p}{L_y} = \frac{150 \mu_f}{d_p^2} \frac{\varepsilon^2}{(1-\varepsilon)^3} V_{SF} \quad (1) $$

valable en régime laminaire et qui s'exprime par la relation : $Re/\varepsilon \ll 1$,

$$ \frac{\Delta p}{L_y} = \frac{1,75 \rho_f}{d_p} \frac{\varepsilon}{(1-\varepsilon)^3} V_{SF}^2 \quad (2) $$

valable en régime turbulent et qui s'exprime par la relation : $Re/\varepsilon \ll 1$, où Re désigne le nombre de Reynolds particulière défini comme suit :

$$ Re = \frac{d_p V_{SF} \rho_f}{\mu_f} \quad (3) $$

et où $d_p$, $\rho_f$, $\mu_f$, $V_{SF}$, $L_y$ et $\varepsilon$ désignent respectivement le diamètre des particules, la masse volumique du fluide, la viscosité du fluide, la vitesse superficielle du fluide, la hauteur du milieu granulaire et la fraction volumique du solide.

La somme des expressions de ces deux lois permet d'obtenir la loi de perte de charge d'Ergun donnée ci-dessous et vérifiée expérimentalement :

$$ \frac{\Delta p}{L_y} = \frac{150 \mu_f}{d_p^2} \frac{\varepsilon^2}{(1-\varepsilon)^3} V_{SF} + \frac{1,75 \rho_f}{d_p} \frac{\varepsilon}{(1-\varepsilon)^3} V_{SF}^2 \quad (4) $$

Les termes linéaire et quadratique avec la vitesse superficielle du fluide $V_{SF}$ correspondent respectivement aux pertes d'énergie visqueuse et cinétique du fluide.

En introduisant la pression $p^*$ et la longueur $L_y^*$ adimensionnées :

$$ p = V_{SF}^2 \rho_f p^* \quad (5) $$
et :

$$L_y = d_p L_y^* \quad (6)$$

la loi d'Ergun s'écrit sous la forme adimensionnelle suivante :

$$\frac{\Delta p^*}{L_y^*} = \frac{150}{Re} \frac{\epsilon^2}{(1-\epsilon)^3} + 1,75 \frac{\epsilon}{(1-\epsilon)^3} \quad (7)$$

où les seuls paramètres à considérer sont :

- le nombre de Reynolds particulière Re ;

- la fraction volumique de solide du lit ε ;

- la hauteur du milieu granulaire rapportée au diamètre des particules $L_y^*$.

Dans ce qui suit, les différentes grandeurs sont sans dimension.

## 2 MODÈLE NUMÉRIQUE

Le modèle numérique consiste à résoudre l'écoulement du fluide interstitiel à travers un milieu poreux constitué de particules sphériques homogènes par simulation directe des équations de Navier-Stokes.

### 2.1 Constitution du milieu poreux

Le milieu poreux considéré est un empilement de type "cubique centré" de $n_p = n_x n_y n_z$ particules sphériques homogènes. $n_x, n_y$ et $n_z$ désignent respectivement le nombre de particules présentes dans chaque direction d'espace x, y et z (fig. 1).

Pour ce type d'empilement, la fraction de solide est donnée par :

$$\epsilon = \frac{n_p \times \frac{4}{3} \pi r^3}{L_x L_y L_z} = \frac{\pi}{6} \approx 0.524 \quad (8)$$

où r désigne le rayon des particules.

### 2.2 Modélisation de l'écoulement du fluide

Les hypothèses suivantes sont effectuées :

- le fluide est incompressible,

- le fluide est newtonien.

La modélisation de l'écoulement du fluide interstitiel repose sur la résolution directe des équations de Navier-Stokes à l'échelle de la particule, écrites sous forme adimensionnelle :

$$\frac{\partial v_i}{\partial t} = -A(\vec{v}) - \frac{\partial p}{\partial x_i} + \frac{1}{Re} \Delta(\vec{v}) \quad (9)$$

$$\frac{\partial v_i}{\partial x_i} = 0 \quad (10)$$

avec :

$$A(\vec{v}) = \frac{\partial(v_i v_j)}{\partial x_j} \quad (11)$$

et :

$$\Delta(\vec{v}) = \frac{\partial^2 v_i}{\partial^2 x_j} \quad (12)$$

où $i = 1, ..., 3$ et $j = 1, ..., 3$.

<figure><img src="image_3.png" /><figcaption>Figure 1
Milieu poreux de type "cubique centré".
"Cubic centered porous media.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2
Conditions aux limites.
Boundary conditions.</figcaption></figure>
Pour les conditions aux limites à la frontière du domaine, une vitesse d'entrée et une pression de sortie sont respectivement imposées dans les plans y = 0 et y = L<sub>y</sub>. Les autres plans sont soumis à des conditions périodiques (fig. 2).

Le mouvement du fluide sur la surface des particules est soumis aux conditions d'adhérence, à savoir que la vitesse du fluide en ces points est égale à la vitesse des particules, nulle dans le cas présent.

## 2.3 Schéma numérique

### 2.3.1 Discrétisation en temps

La méthode de discrétisation temporelle utilisée est une méthode à pas fractionnaire explicite du premier ordre. Elle a été proposée par Chorin [14]. La description de cette méthode est donnée ci-après :

$$ \frac{v_i^{n+1} - v_i^n}{\Delta t_f} = -\frac{\partial p^{n+1}}{\partial x_i} - A(v_i) + \frac{1}{Re} \Delta(v_i^n) \quad (13) $$

où i = 1, ..., 3.

La quantité $v_i^{n+1/2}$ est introduite de la façon suivante :

$$ \frac{(v_i^{n+1} - v_i^{n+1/2}) - (v_i^n - v_i^{n+1/2})}{\Delta t_f} = -\frac{\partial p^{n+1}}{\partial x_i} - A(\vec{v}^n) + \frac{1}{Re} \Delta(\vec{v}^n) \quad (14) $$

ce qui permet d'écrire :

$$ \frac{v_i^n - v_i^{n+1/2}}{\Delta t_f} = -A(\vec{v}^n) + \frac{1}{Re} \Delta(\vec{v}^n) \quad (15) $$

et :

$$ \frac{v_i^{n+1} - v_i^{n+1/2}}{\Delta t_f} = \frac{\partial p^{n+1}}{\partial x_i} \quad (16) $$

L'application de l'opérateur de divergence à l'équation (16) permet d'obtenir l'équation de Poisson pour la pression en tenant compte de la contrainte d'incompressibilité :

$$ \frac{\partial^2 p^{n+1}}{\partial^2 x_i} = \frac{1}{\Delta t_f} \frac{\partial v_i^{n+1/2}}{\partial x_i} \quad (17) $$

La condition aux limites pour la pression s'obtient en effectuant le produit scalaire de l'équation (16) avec la normale unitaire sortante de la frontière du domaine. Par exemple, dans le cas où la frontière correspond au plan x = 0, la condition aux limites en pression s'écrit :

$$ \left( \frac{\partial p^{n+1}}{\partial x_i} \right)_{x=0} = \frac{1}{\Delta t_f} (v_{1\Gamma}^{n+1} - v_{1\Gamma}^{n+1/2}) \quad (18) $$

Le schéma étant indépendant de la valeur de $v_i^{n+1/2}$ sur la frontière [13], une possibilité est de prendre :

$$ v_{1\Gamma}^{n+1/2} = v_{1\Gamma}^{n+1} \quad (19) $$

En définitive, la condition aux limites pour la pression s'écrit sur la frontière correspondant au plan x = 0 :

$$ \left( \frac{\partial p^{n+1}}{\partial x_i} \right)_{x=0} = 0 \quad (20) $$

L'algorithme correspondant est le suivant :

- les composantes de la vitesse du fluide : $v_i^{n+1/2}$ où i = 1, ..., 3 sont déterminées explicitement par l'équation (15);

- la pression $p^{n+1}$ est déterminée par résolution de l'équation (17) par l'algorithme du gradient conjugué stabilisé et préconditionné [15];

- les composantes de la vitesse du fluide : $v_i^{n+1/2}$ où i = 1, ..., 3 sont déterminées explicitement par l'équation (16).

### 2.3.2 Discrétisation en espace

La discrétisation spatiale utilise un maillage décalé en différences finies de type *Marker and Cell* (MAC), où la pression est située au centre de la maille et les composantes de la vitesse sur les milieux des faces (fig. 3).

La discrétisation des termes convectif et de diffusion est prise au deuxième ordre. Les formules sont données dans [13].

<figure><img src="image_5.png" /><figcaption>Figure 3</figcaption></figure>

Maillage de type MAC, le symbole + correspond à la composante de la vitesse $v_1$ aux nœuds $(k + 1/2, l, m)$; le symbole × correspond à la composante de la vitesse $v_2$ aux nœuds $(k, l + 1/2, m)$; le symbole • correspond à la pression aux nœuds $(k, l, m)$.

MAC mesh, + corresponds to $v_1$ at node $(k + 1/2, l, m)$; × corresponds to $v_2$ at grid node $(k, l + 1/2, m)$; • corresponds to pressure at grid node $(k, l, m)$.
### 2.3.3 Conditions de stabilité

Le pas de temps est déterminé de façon à vérifier les conditions de stabilité suivantes :

$$ \Delta t_f \le \frac{h^2}{6} Re $$

et :

$$ \Delta t_f \le \frac{2}{(\tilde{v}_1^2 + \tilde{v}_2^2 + \tilde{v}_3^2)Re} \quad (21) $$

où $\tilde{v}_1^2, \tilde{v}_2^2$ et $\tilde{v}_3^2$ désignent les composantes maximales de la vitesse du fluide. La démonstration de ces formules est effectuée en annexes.

## 3 RÉSULTATS DE SIMULATIONS

### 3.1 Étude de sensibilité à la taille de la maille

Ces simulations consistent à calculer la perte de charge en faisant varier le nombre de nœuds par rayon de particule $n_h$ dans le but de déterminer une taille de maille :

$$ h = \frac{r}{n_h} \quad (22) $$

donnant une erreur faible avec la loi d'Ergun (fig. 4), l'erreur étant définie par la formule suivante :

$$ \text{Erfour} = \frac{\left| \frac{\Delta p}{L_y} \text{ calculée} - \frac{\Delta p}{L_y} \text{ loi d'Ergun} \right|}{\max \left( \frac{\Delta p}{L_y} \text{ calculée}, \frac{\Delta p}{L_y} \text{ loi d'Ergun} \right)} \times 100 \quad (23) $$

<figure><img src="image_6.png" /><figcaption>Figure 4</figcaption></figure>

Évolution de l'erreur en fonction du nombre de nœuds par rayon de particules. Les symboles □ représentent les points de simulation.

Error as a function of number of grid nodes.

Pour cette étude, le nombre de Reynolds particulière Re est fixé à Re = 1 et le volume de matériau granulaire contient $n_p = 4^3$ particules. La hauteur totale du domaine dans la direction y est de 8 diamètres de particules.

Un nombre de nœuds par rayon de particules égale à 16 correspondant à $n = 0,0625$ r permet d'obtenir un écart de moins de 5 % avec la loi d'Ergun (fig. 4) et sera adopté pour les simulations effectuées aux paragraphes 3.3 et 3.4.

### 3.2 Profils de pression

Les profils de pression axiaux et transversaux à la tête du lit ont été représentés dans le cas où Re = 1 et pour des volumes de matériau granulaire constitués de $n_p = 4^3$ et de $n_p = 1$ particules (fig. 5 et 6). La hauteur totale du domaine dans la direction y est de 8 diamètres de particules.

<figure><img src="image_7.png" /><figcaption>Figure 5</figcaption></figure>

Profils de pression axiaux pour des volumes constitués de $n_p = 4^3$ et $n_p = 1$ particules à $x = L_x^{1/2}$ et $z = L_z^{1/2}$ fixés.

Pressure profil for $n_p = 4^3$ and $n_p = 1$ particles and at $x = L_x^{1/2}$ and $z = L_z^{1/2}$ fixed.
<figure><img src="image_8.png" /><figcaption>Figure 6</figcaption></figure>

Profils de pression transversaux pour des volumes constitués $n_p = 4^3$ et $n_p = 1$ particules à $y = L_y$ et $z = L_z^{1/2}$ fixés.

Pressure profil for $n_p = 4^3$ and $n_p = 1$ particles at $y = L_y$ and $z = L_z^{1/2}$ fixed.

Pour les deux volumes de matériau granulaire, et en conformité avec la loi, la variation de pression axiale est linéaire avec la hauteur $L_y$ tandis que la variation de pression transversale est négligeable devant la variation de pression axiale. Les oscillations des profils transversaux correspondent à la présence des particules.

### 3.3 Volume élémentaire représentatif et validité du caractère local

Des simulations sont ensuite réalisées pour différents volumes de matériau granulaire définis par $n_p = 4^3$, $n_p = 3^3$ et $n_p = 1$ particules et pour différentes valeurs du nombre de

Reynolds particulaire allant de $Re = 0,02$ à $Re = 20$ (fig. 7). Les nombres de Reynolds particulaires plus importants sont difficilement accessibles avec le modèle utilisé où la description du champ hydrodynamique se fait par des équations non moyennées de Navier-Stokes. La limitation de cette approche provient du fait que la résolution directe des équations directes de Navier-Stokes nécessite une taille de maille rapportée au rayon des particules inversement proportionnelle au nombre de Reynolds particulaire, donnée par la relation :

$$ \frac{h}{r} = O\left(Re^{-\frac{9}{4}}\right) \qquad (24) $$

et conduisant à des temps CPU prohibitifs.

La loi de perte de charge d'Ergun a été retrouvée pour les 3 volumes envisagés et pour un nombre de Reynolds particulaire $Re \le 20$ avec une erreur inférieure à 5 %. Le volume élémentaire représentatif est de l'ordre de grandeur du diamètre des particules.

<figure><img src="image_9.png" /><figcaption>Figure 7</figcaption></figure>

Perte de charge en fonction du rapport $Re/\epsilon$ pour des volumes de matériau granulaire constitués de $n_p = 4^3$, $n_p = 3^3$ et $n_p = 1$ particules.

Pressure drop as a function of $Re/\epsilon$ for different volume of granular media.

### 3.4 Temps CPU et limitations

Le temps CPU est de l'ordre de 6 µs par pas de temps et par nœuds sur une Octane SGI à processeur R10 000. Les limitations de la méthode proviennent :

- du nombre de Reynolds particulaire qui, compte tenu de la contrainte imposée par la relation (22), est inférieur à 20 pour une taille de maille de $h = 0,0625 r$ ;
- de la complexité du calcul qui est en $O(k^5)$, où k désigne le facteur de réduction du maillage ($O(k^2)$ pour la condition de stabilité et $O(k^3)$ pour l'augmentation du nombre de nœuds).

À titre indicatif, les simulations correspondant à un volume de matériau granulaire de $n_p = 4^3$ particules et effectuées avec des tailles de maille respectives de $h = 0,25$ r et $h = 0,0625$ r utilisent respectivement 65 536 et 4 194 304 nœuds et nécessitent respectivement 1000 et 16 000 pas de temps pour atteindre un état stationnaire. Le temps CPU respectif est de l'ordre de 7 minutes et 111 heures.

## CONCLUSION

Un modèle numérique, basé sur la résolution directe des équations de Navier-Stokes est proposé pour calculer les pertes de charge d'un fluide visqueux incompressible à travers un empilement de particules sphériques. Des simulations effectuées pour différents volumes de matériau granulaire permettent de retrouver la loi de perte de charge d'Ergun pour des nombres de Reynolds particulaires $Re \le 20$. Les résultats montrent que le volume élémentaire représentatif est de l'ordre de grandeur du diamètre des particules.

En perspective, cette approche doit nous permettre, en incluant le calcul des trajectoires, de justifier également le caractère local d'autres lois phénoménologiques, comme la loi de Richardson et Zaki [16] valable pour les lits fluidisés, les lits mobiles et la sédimentation. Ce travail pourra s'effectuer selon une méthode de couplage fluide-solide proposée aux références [17-19].

## REMERCIEMENTS

V. Komiwes tient à remercier K. Hofler, M. Muller, B. Wachmann et S. Schwarzer du laboratoire ICAI de l'université de Stuttgart pour les diverses discussions, conseils et l'aide apportée, ainsi qu'à S. Guyons et A. Stoessel de l'Institut français du pétrole pour le travail réalisé concernant l'optimisation du code de calcul.

## RÉFÉRENCES

1 Dullien, F.A.L. (1979) *Porous Media: Fluid Transport and Pore Structure*, Academic Press, New York, London.

2 Ergun, S. (1952) *Fluid Flow Through Packed Columns*. *Chemical Engineering Progress*, **48**, 2.

3 Adler, P.M., Jacquin, C.G. and Quibler, J.A. (1990) *Flow in Simulated Porous Media*. *Int. J. Multiphase Flow*, **16**, 4.

4 Martys, N.S., Torquato, S. and Bentz, D.P. (1994) *Universal Scaling of Fluid Permeability*. *Sphere Packings Physical Review E*, **50**, 1.

5 Schwartz, L.M., Martys, N.S., Bentz, D.P., Garboczi, E.J. and Torquato, S. (1993) *Cross-Property Relations and Permeability Estimation*. *Model Porous Media Physical Review E*, **48**, 6.

6 Simonin, O. (1990) *Eulerian Formulation for Particle Dispersion in Turbulent Two-Phase Flows*. *5th Workshop on Two Phase Flow Predictions*, 19-22, Erlangen, March 1990.

7 Boelle, A. (1997) *Validation d'un modèle à deux fluides appliqué à la simulation des lits fluidisés denses*. *Thèse de doctorat*, université de Paris VI.

8 Loo, S. (1967) *Fluid Dynamics of Multiphase Systems*, Baidsdell Publishing Co.

9 Ferschneider, G. and Mège, P. (1996) *Eulerian Simulation of Dense Phase Fluidized Bed*. *Revue de l'IFP*, Éditions Technip, Paris, **51**, 2.

10 Bird, R.B., Stewart, W.E. and Lightfoot, E.N. *Transport Phenomena*. John Wiley Sons, New York, Chichester, Bisbane, Toronto.

11 Tsuji, Y. and Tanaka, T. (1993) *Discrete Particle Simulation of Flow Patterns in Two-Dimensional Gas Fluidized Beds*. *International Journal of Modern Physics B*, **7**, 9 and 10, 1889-1898.

12 Hoomans, B.P.B., Kuipers, J.A.M., Briels, W.J. and Van Swaaij, W.P.M. (1993) *Discrete Particle Simulation of Bubble and Slug Formation in a Two-Dimensional Gas-Fillidized-Bed: a Hard-Sphere Approach*. *Chem. Eng. Sci.*, **51**, 99-118.

13 Peyret, R. and Taylor, T.D. (1983) *Computational Methods for Fluid Flow Springer Series*. *Computational Physics*, Springer, New York, Berlin, Heidelberg.

14 Chorin, A.J. (1967) *A Numerical Method for Solving Incompressible Viscous Flow Problems*. *Journal of Computational Physics*, **2**, 12-26.

15 Van Der Vorst, H.A. (1992) *Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems* SIAMJ. *Sci. Statist. Comput.*, **13**, 631-644.

16 Richardson, J.F. and Zaki, W.N. (1954) *Sedimentation and Fluidisation*. *Trans. Instn. Chem. Eng.*, **32**.

17 Kalthoff, W., Wachmann, B., Shwarzer, S. and Herrmann, H.J. *Simulating Particles in Fluids in Two and Three Dimensions*. *Communication personnelle*.

18 Kalthoff, W., Shwarzer, S. and Herrmann, H.J. *An Algorithm for the Simulation of Particulate Suspensions with Inertia Effects*. *Communication personnelle*.

19 Komiwes, V., Herrmann, H.J., Mège, P. and Meimon, Y. (1998) *Direct Simulation of Granular Flow with Fluid, Session 20, Granular Materials of the Joint Conferences EPS. 17th Meeting of the Consensed Matter Division*; SFP: 6ᵉ Journée de la matière condensée, Grenoble.
# ANNEXES

## Conditions de stabilité

Une possibilité pour étudier la stabilité du schéma est de négliger les termes en pression et de linéariser le terme convectif en considérant que :

$$v_1(k+1/2,l,m) = \tilde{v}_1(\text{const})$$

$$v_2(k,l+1/2,m) = \tilde{v}_2(\text{const})$$

$$v_1(k,l,m+1/2) = \tilde{v}_3(\text{const})$$

soit pour la première composante du terme de convection :

$$\tilde{A}(\tilde{v}_1) = \tilde{v}_1 \frac{v_1(k+3/2,l,m) - v_1(k-1/2,l,m)}{2h} + \tilde{v}_2 \frac{v_1(k+1/2,l+1,m) - v_1(k+1/2,l-1,m)}{2h} + \tilde{v}_3 \frac{v_1(k+1/2,l,m+1) - v_1(k+1/2,l,m-1)}{2h}$$

Le schéma linéarisé est alors donné pour la première composante par :

$$\tilde{v}_1^{n+1}(k+1/2,l,m) = \tilde{v}_1(k+1/2,l,m) + \Delta t_f \times (\tilde{A}(\tilde{v}_1) + \frac{1}{Re} \Delta(\tilde{v}_1))$$

Les autres composantes s'obtiennent de manière identique. L'injection de la perturbation :

$$\tilde{v}_1^{n+1}(k+1/2,l,m) = \alpha^n e^{i[p(k+1/2)h+qlh+rmh]}$$

à l'intérieur du schéma linéarisé permet d'obtenir le facteur d'amplification $\alpha$ pour $p, q$ et $r$ arbitraire :

$$\alpha = 1 - \frac{\Delta t_f}{h} [i\tilde{v}_1 \sin(ph) + i\tilde{v}_2 \sin(qh) + i\tilde{v}_3 \sin(rh)] + \frac{2\Delta t_f}{h^2 Re} [\cos(ph) + \cos(qh) + \cos(rh) - 3]$$

où : $i^2 = -1$

La condition suivante est suffisante pour respecter la condition de stabilité :

$$|\alpha|^2 \le 1$$

quel que soit $p, q$ et $r$ :

$$\frac{6\Delta t_f}{h^2 Re} \le 1 \qquad (25)$$

et :

$$\frac{(\tilde{v}_1^2 + \tilde{v}_2^2 + \tilde{v}_3^2) Re \Delta t_f}{2} \le 1 \qquad (26)$$

La démonstration des conditions suffisantes de stabilité du schéma numérique est effectuée ci-dessous.

## Démonstration

$$|\alpha|^2 = [1 + a(\tilde{p} + \tilde{q} + \tilde{r} - 3)]^2 + b[\tilde{v}_1 \sin(ph) + \tilde{v}_2 \sin(qh) + \tilde{v}_3 \sin(rh)]^2$$

en posant :

$$a = \frac{2\Delta t_f}{h^2 Re}, \quad b = \frac{(\Delta t_f)^2}{h^2}, \quad \tilde{p} = \cos(ph), \quad \tilde{q} = \cos(qh), \quad \tilde{r} = \cos(rh)$$

D'après l'inégalité de Schwarz :

$$[\tilde{v}_1 \sin(ph) + \tilde{v}_2 \sin(qh) + \tilde{v}_3 \sin(rh)]^2 \le (\tilde{v}_1^2 + \tilde{v}_2^2 + \tilde{v}_3^2)[\sin^2(ph) + \sin^2(qh) + \sin^2(rh)]$$

et en utilisant l'inégalité suivante obtenue d'après (25) et (26) :

$$b(\tilde{v}_1^2 + \tilde{v}_2^2 + \tilde{v}_3^2) \le a$$

ceci implique :

$$|\alpha|^2 \le [1 + a(\tilde{p} + \tilde{q} + \tilde{r} - 3)]^2 + a(3 - \tilde{p}^2 - \tilde{q}^2 - \tilde{r}^2) = 1 + a\Phi(\tilde{p}, \tilde{q}, \tilde{r})$$

avec :

$$\Phi(\tilde{p}, \tilde{q}, \tilde{r}) = (a-1)(\tilde{p}^2 + \tilde{q}^2 + \tilde{r}^2) + 2a(\tilde{p}\tilde{q} + \tilde{p}\tilde{r} + \tilde{q}\tilde{r}) + 2(1-3a)(\tilde{p} + \tilde{q} + \tilde{r}) + 3(3a-1) \\ \Phi(\tilde{p}, \tilde{q}, \tilde{r}) = (\tilde{p} \quad \tilde{q} \quad \tilde{r}) \underbrace{\begin{pmatrix} a-1 & a & a \\ a & a-1 & a \\ a & a & a-1 \end{pmatrix}}_{=Q} \begin{pmatrix} \tilde{p} \\ \tilde{q} \\ \tilde{r} \end{pmatrix} + 2(1-3a)(\tilde{p} + \tilde{q} + \tilde{r}) + 3(3a-1)$$

Le polynôme caractéristique de $Q$ est :

$$P_Q(\lambda) = (1+\lambda)^2(3a-1-\lambda)$$

Les valeurs propres de $Q$ sont :

$$\lambda = -1 \text{ et } \lambda = 3a - 1$$
La condition (25) implique que ces valeurs sont négatives.
$\Phi$ étant définie négative, elle admet un maximum qui vérifie :

$$ \frac{\partial \Phi}{\partial \tilde{p}} = 2a(\tilde{p} + \tilde{q} + \tilde{r}) - 2\tilde{p} + 2(1 - 3a) = 0 $$

$$ \frac{\partial \Phi}{\partial \tilde{q}} = 2a(\tilde{p} + \tilde{q} + \tilde{r}) - 2\tilde{q} + 2(1 - 3a) = 0 $$

$$ \frac{\partial \Phi}{\partial \tilde{r}} = 2a(\tilde{p} + \tilde{q} + \tilde{r}) - 2\tilde{r} + 2(1 - 3a) = 0 $$

soit :

$$ \tilde{p} = \tilde{q} = \tilde{r} = 1 $$

En ce point :

$$ \Phi(1, 1, 1) = 0 $$

d'où, quel que soit $\tilde{p}, \tilde{q}, \tilde{r}$ :

$$ \Phi(\tilde{p}, \tilde{q}, \tilde{r}) \le 0 $$

soit finalement :

$$ |\alpha|^2 \le 1 $$