# Imbibition capillaire d'un échantillon de bois : reconstitution 3D de l'anatomie et simulation Lattice Boltzmann

Xavier Frank, Giana Almeida, Patrick Perré

► To cite this version:

Xavier Frank, Giana Almeida, Patrick Perré. Imbibition capillaire d'un échantillon de bois : reconstitution 3D de l'anatomie et simulation Lattice Boltzmann. CFM 2009 - 19ème Congrès Français de Mécanique, Aug 2009, Marseille, France. hal-03390796

HAL Id: hal-03390796

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

Submitted on 21 Oct 2021

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
# Imbibition capillaire d'un échantillon de bois : reconstitution 3D de l'anatomie et simulation Lattice Boltzmann

X. FRANKª, G. ALMEIDAᵇ, P. PERRÉᶜ

a. INRA, UMR 1092 LERFOB, F-54042 Nancy cedex, France

b. Luiz de Queiroz College of Agriculture. University of Sao Paulo (ESALQ, USP), Brazil

c. AgroParisTech, UMR 1092 LERFOB, ENGREF 14, rue Girardet F-54042 Nancy cedex, France

E-mail: xfrank@nancy-engref.inra.fr

## Résumé :

Le bois est un milieu poreux hydrophile, il est donc crucial de tenir compte de ses interactions avec l'eau. Nous nous intéresserons plus particulièrement à la dynamique d'imbibition capillaire d'un échantillon de bois. Celle-ci est étudiée numériquement par une approche Lattice Boltzmann diphasique. L'anatomie 3D de l'échantillon est reconstituée par interpolation au moyen du logiciel Meshpore à partir d'une série d'images en microscopie optique de coupes 2D réalisées au microtome.

## Abstract :

As wood is an hydrophilic material, it is crucial to take it's interactions with water into account. Capillary imbibition dynamics inside a wood sample is adressed here. To do so, numerical experiments are carried out with the help of a multiphase Lattice Boltzmann scheme. An interpolation routine, implemented in the Meshpore software, is used to deduce wood anatomy in three dimensions from micrographs of 2D sections of a real wood sample.

**Mots clefs :** milieu poreux, capillarité, Lattice Boltzmann, bois

## Introduction

Les écoulements multiphasiques en milieux poreux ont pendant longtemps été appréhendés par des paramètres macroscopiques (perméabilités relatives, pression capillaire)[1]. En effet, il y a encore une ou deux décennies, prendre en considération à la fois la morphologie tridimensionnelle du milieu poreux et la répartition et l'écoulement des différentes phases fluides présentes dans ce milieu n'était qu'un rêve. Aujourd'hui, les progrès remarquables effectués sur différents fronts de science (dispositifs d'acquisition, puissance des ordinateurs, méthodes numériques...) ont donné vie à ce rêve [2, 3].

Le travail présenté dans ce texte s'incrit dans ces travaux de modélisation des phénomènes à l'échelle des pores. Nous avons travaillé sur une essence de bois feuillu (Fagus sylvatica). L'anatomie 3D a été définie par une succession de coupes réalisées à l'aide d'un microtome à glissière. Chaque image a ensuite été traitée afin de décrire le système vaculaire par un ensemble de contours fermés. Ces coupes 2D sériées ont ensuite été analysées simultanément afin de reconstituer la morphologie tridimensionnelle du système vasculaire de l'arbre. L'invasion de ce système vasculaire a finalement été simulée par une approche Lattice Boltzmann diphasique.

## 1 Acquisition et traitement des images

Cette étude a été réalisée sur du bois de hêtre (Fagus sylvatica). Cette essence produit, par croissance secondaire, un bois à porosité diffuse qui possède des vaisseaux d'environ 50 µm de diamètre. L'exploration anatomique a été faite sur un barreau de 10 × 5 × 30 mm³ suivant les directions radiale, tangentielle et longitudinale, respectivement. Ce barreau a été découpé en tranches successives grâce à un microtome. Afin d'acquérir une image de l'échantillon après chaque découpe, nous avons utilisé un microtome à glissière. Un montage a été réalisé au laboratoire afin de fixer sur ce microtome un microscope optique modulaire avec éclairage épicentrique (Nachet MS 98) équipé d'une caméra digitale Basler Vision Technologies. Le microtome utilisé (ancien modèle de chez Leitz) a pour particularité d'avoir une avance de l'échantillon parfaitement verticale. Ceci permet au microscope d'observer la même zone de l'échantillon au fur et à mesure des découpes. Des sections microscopiques de 20 µm d'épaisseur (L) ont été réalisées et des images ont été acquises une coupe sur deux, donc avec un intervalle d'exploration de 40 µm. Un total de 58 images a été analysé avec une épaisseur d'exploration de 2,28 mm (L).
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Étapes du traitement d'images : (A) image initiale ; (B) Sélection de la zone d'intérêt ; (C) segmentation du système vasculaire avec le logiciel *MeshPore* ; (D) Contours vectorisés</figcaption></figure>

Le traitement des images a été effectué grâce au logiciel *MeshPore*, une application Windows développée au laboratoire en Fortran 95 avec la bibliothèque Winteracter pour l'interface et les sorties graphiques [4]. Ce logiciel a pour particularité de réaliser tout le traitement d'image sous forme vectorielle, ce qui permet ultérieurement de générer des grilles de points avec une résolution indépendante de la taille initiale des pixels de l'image. La première opération réalisée par *MeshPore* consiste à générer automatiquement, sous la forme de chaînes de points, les contours d'isovaleurs d'un niveau de gris choisi. Ces chaînes peuvent ensuite être modifiées manuellement (problèmes locaux de qualité d'image) et rééchantillonnées (répartition régulière des points du contours, choix de la résolution spatiale, abandon des objets trop petits...). En raison de la grande quantité de vaisseaux présents dans chaque image, un langage macro disponible dans *MeshPore* a permis de choisir la même région d'intérêt et les traitements automatiques pour toutes les images étudiées. La Figure 1 montre les principales étapes du traitement d'images : (1A) image initiale acquise par la camera digitale ; (1B) région délimité pour faire le traitement d'image ; (1C) traitement d'image en utilisant le logiciel *MeshPore* ; (1D) image finale qui montre la zone d'étude (carré jaune).

## 2 Reconstruction de la morphologie 3D

Afin de passer de n plans 2D à une véritable morphologie 3D, deux étapes crucales sont nécessaires : i) trouver les connectivités entre les chaînes de deux plans successifs et ii) choisir une méthode d'interpolation entre les plans afin de générer un réseau cubique. Pour franchir ces deux étapes, une version 3D de *MeshPore* a été développée. Elle permet de charger tous les plans 2D et de faire un traitement de l'ensemble de ces plans. La partie graphique de cette version 3D permet de visualiser à loisir un plan ou un ensemble de plans en choisisant le plan de départ, le nombre et l'increment. Un menu permet aussi de choisir les clusters visualisés. Les projections 3D sont réalisés en utilisant les quaternions, qui permettent une formulation élégante. De nombreuses options sont envisageables pour trouver les connectivités. Plutôt qu'un critère basé sur l'aire commune entre deux contours [5], nous avons préféré ici un critère plus exigeant pour décider si deux contours de plans successifs sont connectés : deux contours sont considérés connectés si le barycentre de l'un est à l'intérieur de l'autre, ou vice versa. La réciprocité du critère est importante, notamment pour les branchements (deux vaisseaux d'un plan connectés à un seul du plan supérieur ou inférieur). Conformément au théorème de Stokes, le barycentre est calculé par une intégrale de contour sur la chaîne (1,2). L'intégrale sur le contour de l'angle de l'arc partant du point permet de savoir si ce point est à l'intérieur du contour. Pour ces intégrales de contours, les chaînes ouvertes sont fermées en reliant dernier et premier point. Les clusters (ensemble des contours connectés) sont construits de proches en proches à partir de ces connections. Des traitements sur les clusters sont ainsi possibles (par exemple, la sélection de tous les clusters traversant l'échantillon). La figure 2 représente les contours de 20 plans (un plan tracé tous les deux plans) de 3 clusters traversant entièrement l'échantillon.
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Exemple de 3 clusters de contours tracés sur 20 plans successifs.</figcaption></figure>

$$x_{bary} = \frac{1}{S} \oint_C x^2 dy \qquad (1)$$

$$y_{bary} = -\frac{1}{S} \oint_C y^2 dx \qquad (2)$$

L'autre difficulté consiste à interpoler la morphologie entre deux plans successifs. Cette opération est en effet nécessaire pour obtenir un réseau cubique car les distances selon l'axe Oz entre plans successifs (40 µm) sont plus grandes que la résolution typique nécessaire dans le plan Oxy pour représenter correctement la morphologie des vaisseaux. Afin de respecter au mieux la structure vasculaire de l'arbre, ce sont les chaînes qui sont interpolées avant de générer le réseau. Pour chaque couple de chaînes connectées (chaînes initiales) de deux plans successifs, une chaîne est construite pour chacun des plans d'interpolation. Cette chaîne est construite point par point. Pour chaque point la première chaîne initiale, le "meilleur" point de l'autre chaîne est sélectionné pour former une génératrice du cheminement. Un point de la chaîne interpolée est alors calculé comme intersection de cette génératrice avec le plan d'interpolation. Le "meilleur" point de la seconde chaîne initiale est choisi comme étant celui qui est le plus proche du plan défini par les deux barycentres et le point de la première chaîne (Fig. 3). Le résultat de cette interpolation est visible sur la figure 4 qui montre deux exemples typiques de connections : deux chaînes reliées légèrement déformées et décalées dans le plan Oxy et un exemple de deux chaînes très proches connectées à une chaîne unique plus grande (cas de deux vaisseaux connectés à un seul).

## 3 Approche Lattice Boltzmann

La simulation numérique des phénomènes de mouillage est une tâche des plus complexes, en particulier dans le cas des milieux poreux réels, dont la géométrie rend d'autant plus difficile l'implémentation des conditions aux limites avec les approches classiques, du type suivi d'interface. L'approche Lattice Boltzmann constitue une alternative crédible aux méthodes classiques de la mécanique des fluides numérique pour répondre à ces difficultés. En particulier, les interfaces émergent naturellement, ce qui évite la très délicate opération de suivi des interfaces. Issue récemment des la physique statistique des systèmes hors équilibre, l'approche Lattice Boltzmann ne repose pas sur une résolution de l'équation de Navier-Stokes, mais s'appuie sur une description statistique du système considéré [6]. On définit les fonctions densité de probabilité de présence (PPDF) $f_i$, où $f_i(\vec{r}, t)$ est la probabilité de trouver au point $\vec{r}$ à la date $t$ une particule se déplaçant à la vitesse $\vec{c}_i$. Les particules se déplacent sur un réseau, aussi les seules vitesses $\vec{c}_i$ possibles sont celles qui permettent à une particule de se déplacer d'un noeud à un noeud voisin en un pas de temps. L'ensemble des vitesses $\vec{c}_i$ définissent donc le réseau. Ici, nous utiliserons le réseau D3Q19, un réseau tridimensionnel comportant 19 vitesses. L'évolution des $f_i$ obéit à l'équation Lattice Boltzmann (3).

$$f_i(\vec{r} + \delta t \vec{c}_i, t + \delta t) - f_i(\vec{r}, t) = -\frac{1}{\tau} (f_i - f_i^{eq}) \qquad (3)$$

Les distributions d'équilibre $f_i^{eq}$ sont déduites des grandeurs macroscopiques (4-7).
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Principe développé pour la génération de chaînes sur les plans d’interpolation. Ce procédé est répété pour toutes les chaînes connectées entre les plans 1 et 2.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Exemple de 9 plans d’interpolation construits entre les contours successifs de deux clusters : réseau de points obtenus pour les 2ème (haut) et 5ème (bas) plans d’interpolation.</figcaption></figure>
$$ \sum_i f_i^{eq} = \rho \qquad (4) $$

$$ \sum_i f_i^{eq} \vec{c}_i = \rho \vec{u} \qquad (5) $$

$$ \sum_i f_i^{eq} c_{i\alpha} c_{i\beta} = \rho u_{\alpha} u_{\beta} + \rho c_s^2 \delta_{\alpha\beta} \qquad (6) $$

$$ \sum_i f_i^{eq} c_{i\alpha} c_{i\beta} c_{i\gamma} = \rho c_s^2 (u_{\gamma}\delta_{\alpha\beta} + u_{\alpha}\delta_{\gamma\beta} + u_{\beta}\delta_{\alpha\gamma}) \qquad (7) $$

La densité $\rho$ et la vitesse locale $\vec{u}$ sont des moments des $f_i$ (8-9).

$$ \rho = \sum_i f_i \qquad (8) $$

$$ \rho \vec{u} = \sum_i f_i \vec{c}_i \qquad (9) $$

Conformément à l'approche de Shan et Chen [7, 8], une vitesse modifiée est calculée, des forces d'interaction fluide-fluide $\vec{F}_{FF}$ et fluide-solide $\vec{F}_{FS}$ agissant en chaque point (10).

$$ \rho \vec{u}' = \rho \vec{u} + \tau (\vec{F}_{FF} + \vec{F}_{FS}) \qquad (10) $$

Les forces sont déduites de pseudo-potentiels d'interaction entre particules (11,12).

$$ \vec{F}_{FF}(\vec{r}) = \psi(\vec{r}) \sum_i G_i^{FF} \psi(\vec{r} + \delta t \vec{c}_i) \qquad (11) $$

$$ \vec{F}_{FS}(\vec{r}) = \psi(\vec{r}) \sum_i G_i^{FS} S(\vec{r} + \delta t \vec{c}_i) \qquad (12) $$

La fonction $\psi = 1 - e^{-\rho}$ peut être interprétée comme une densité équivalente. Le paramètre $S(\vec{r})$ adopte la valeur 1 si le point $\vec{r}$ est dans la phase solide, et la valeur 0 dans le cas contraire. Les valeurs des paramètres $G_i^{FF}$ et $G_i^{FS}$ dépendent de la direction (13,14).

$$ G_i^{FF} = \begin{cases} 2G_i^{FF} & \|\vec{c}_i\| = 1 \\ G_i^{FF} & \|\vec{c}_i\| = \sqrt{2} \end{cases} \qquad (13) $$

$$ G_i^{FS} = \begin{cases} 2G_i^{FS} & \|\vec{c}_i\| = 1 \\ G_i^{FS} & \|\vec{c}_i\| = \sqrt{2} \end{cases} \qquad (14) $$

Au-delà d'une valeur critique de $G^{FF}$, une transition liquide-vapeur se produit spontanément, donnant naissance à un système diphasique. D'autre part, $G^{FS}$ fixe l'angle de contact au niveau de la ligne triple liquide-vapeur-solide. La vitesse modifiée $\vec{u}'$, calculée via l'équation (10), remplace la vitesse $\vec{u}$ dans les équations (5-7). En pratique, ce type d'approche se révèle être particulièrement bien adapté aux problématiques des milieux poreux [2, 3]. Le modèle de Shan et Chen a été implémenté dans un code de calcul parallèle développé au laboratoire et installé sur un cluster de 64 processeurs sous réseau Infiniband.

## 4 Résultats

L'anatomie 3D est décrite sur un réseau $123 \times 146 \times 360$. Afin que les forces capillaires soient les seules forces motrices en présence, la taille du réseau a été doublée dans la direction $z$, l'échantillon virtuel proprement dit n'occupant que la moitié de la fenêtre $z \ge 361$. Une bande de liquide occupant l'autre moitié du domaine de simulation est simplement mise en contact avec la phase solide, et envahit spontanément les vaisseaux (Fig. 5). La simulation est donc réalisée sur un domaine $123 \times 146 \times 720$. Dans le cas présent, notre code de calcul réalise 10000 itérations par heure de calcul.

Si, du fait d'une force motrice plus importante, les vaisseaux étroits sont envahis plus rapidement au début de la simulation, la viscosité contrebalance cet effet à plus long terme à mesure que la longueur de la colonne de liquide augmente. La variabilité naturelle du réseau vasculaire du bois induit donc une variabilité de la vitesse globale d'imbibition d'un vaisseau à l'autre.
<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Simulation Lattice Boltzmann de l’imbition capillaire d’un échantillon virtuel de bois de hêtre. A gauche : après 10 000 itérations, à droite : après 74 000 itérations.</figcaption></figure>

## Conclusion

Ce travail propose une approche complète permettant de passer d'un échantillon de bois à la simulation diphase de l'invasion de cet échantillon par un fluide mouillant. La méthode utilisée pour décrire la morphologie 3D est destructive et fastidieuse, mais elle permet d'avoir la très bonne résolution de la microsocopie optique avec un coup d'équipement relativement réduit. Une résolution équivalente en micro-tomographie nécessite en effet les meilleurs équipements disponibles dans les synchrotrons. L'approche Lattice Boltzmann diphasique retenue ici [7] présente quelques défauts, mais elle est très bien adaptée au calcul parallèle, indispensable pour une représentation suffisamment précise d'une morphologie 3D complexe. De nombreuses applications de ce travail sont envisageables dans tout type de milieu poreux : imbibition, séchage, simulation d'essais de porosimétrie mercure...

## Références

[1] Whitaker S. Simultaneous heat, mass, and momentum transfer in porous media : A theory of drying. Advances in Heat Transfer, 13, 119-203, 1977.

[2] Sukop M. C. and Or D. Invasion percolation of single component, multiphase fluids with lattice boltzmann models. Physica B, 338, 298-303, 2003.

[3] Hatiboglu C. U. and Babadagli T. Pore-scale studies of spontaneous imbibition into oil-saturated porous media. Phys. Rev. E, 77, 066311, 2008.

[4] Perré P. Meshpore : a software able to apply image-based meshing techniques to anisotropic and heterogeneous porous media. Drying technology Journal, 23, 1993-2006, 2005.

[5] Kwon C. S. W., G. and Lee K. J. Automatic generation of tetrahedral meshes from medical images. Computers and Structures, 81, 765-775, 2003.

[6] Succi S. The lattice boltzmann equation for fluid dynamics and beyond. Clarendon Press, Oxford, 2001.

[7] Shan X. and Chen H. Lattice boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E, 47, 1815-1820, 1993.

[8] Martys N. S. and Chen H. Simulation of multicomponent fluids in complex three dimensional geometries by the lattice boltzmann method. Phys. Rev. E, 53, 743-750, 1996.