# Etude d'un Ecoulement dans un Milieu Poreux Saturé
# Utilisé en Substrat de Culture

Saadi Bougoul, Azeddine Soudani et André Jaffrin¹

Département de physique, Faculté des Sciences, Université de Batna, Rue Chahid Mohamed El Hadi Boukhlouf, Batna,
Algérie E-mail: s_bougoul@hotmail.com

¹URIH-INRA Route des Colles; O6410 Biot, France

**Résumé** – Un modèle numérique d'écoulements de solution saline dans un milieu poreux saturé, représentant le substrat de culture parallélépipédique, a été développé pour rendre compte des résultats de visualisations par colorants et prédire l'évolution de la concentration de la solution en chaque point du substrat. Une superposition de puits et de sources a permis de traiter analytiquement le problème de Darcy correspondant à une injection ponctuelle et un point de drainage localisé et elle a donné une illustration des lignes de courant, du champ de vitesse ainsi que la forme du bulbe d'irrigation. Pour interpréter le phénomène de l'élution de la solution ancienne par l'apport de solution nouvelle, l'équation de Navier-Stokes associée à une équation de dispersion a été résolue.

**Abstract** – A numerical study of saline solution flows in a parallelepipedic saturated porous medium, which represents a growing substrate, was developed to account for actual flow visualisations with dies and predict the solution concentration at any point of the growing substrate. In a first approach, Darcy flows corresponding to an injection point to a drainage point were analytically treated by a linear superposition of periodic sources and wells. To interpret the elution phenomenon by a new saline solution to the equation of a Navier-Stokes coupled with the dispersion equation was solved.

**Mots clés** : Substrat de culture – milieu poreux – solution saline – écoulement laminaire – élution – modèles

## 1. INTRODUCTION

Un drainage apparaît à la suite d'un apport de solution nouvelle lorsque celui-ci excède la capacité disponible dans le substrat. L'analyse de ce drainage constitue un élément d'information pratiquement en temps réel. Néanmoins, l'étude de cette solution drainée ne renseigne pas directement sur la distribution de la concentration qui existait préalablement dans le substrat : Le volume drainé est généralement faible par rapport à la capacité totale du substrat et ne représente que la faible couche de liquide initialement présente en partie basse, un phénomène de dispersion affecte le liquide déplacé et en modifie la composition à chaque déplacement. En particulier, la solution neuve ne remplace que progressivement la solution ancienne, tout en balayant un volume de pores apparent plus grand que son propre encombrement.

Dans ces conditions, l'évolution du contenu du substrat lors d'apports successifs n'est pas intuitive et nécessite une modélisation précise du processus de dispersion hydrodynamique qui permettrait alors de mieux gérer la nutrition de la plante en prévoyant la répartition spatiale des éléments minéraux lors d'irrigations successives, et de l'adapter à la demande de la plante qui obéit directement aux fluctuations de son microclimat [8].
## 2. MODELISATION

### 2.1. Modèle des puits et des sources

L'étude bidimensionnelle consiste, en supposant que le milieu étant saturé, à modéliser l'écoulement dans le substrat rectangulaire par l'équation de Darcy linéaire, et à utiliser le principe de superposition d'écoulements simples pour satisfaire les conditions aux limites à savoir une vitesse normale nulle le long des limites du domaine et donc des lignes de courant qui suivent les frontières du substrat. Pour cela, on utilise le modèle de plusieurs puits et sources emprunté à la dynamique des écoulements potentiels à deux dimensions des fluides parfaits, qui satisfont à une équation de Laplace strictement équivalente à la loi de Darcy. En effet le modèle le plus simple, basé sur une source et un puits, qui correspondent respectivement à un point d'apport et de drainage localisés et fixes, ne satisfait pas aux conditions aux limites.

Pour un fluide sans viscosité obéissant à l'équation d'Euler, il existe une classe d'écoulements particulièrement simples, les écoulements sans vorticité, pour lesquels on sait que le vecteur vitesse dérive d'un potentiel $\varphi$. L'analyse de l'écoulement peut être effectuée à l'aide de la seule fonction $\varphi$. Lorsque l'écoulement est incompressible, on peut facilement montrer que cette fonction est solution de l'équation de Laplace [11].

$$ \nabla^2 \varphi = 0 \qquad (1) $$

### 2.2. Loi de Darcy

Les écoulements en milieu poreux ne peuvent pas être considérés comme des écoulements de fluides parfaits à cause de la viscosité. Néanmoins, si l'on fait abstraction de la complexité (à l'échelle de la porosité) de l'écoulement et si l'on ne retient que le mouvement global du liquide, on est amené à une description très voisine de celle des écoulements potentiels. C'est ce que permet de confirmer la loi de Darcy qui exprime sous forme analytique une observation expérimentale [3]. La loi de Darcy exprime une relation de proportionnalité entre la vitesse de filtration et le gradient hydraulique dans la direction de l'écoulement [3].

$$ u = Km \frac{\partial h}{\partial x_i} \qquad (2) $$

$Km$ est la conductivité hydraulique, elle dépend des caractéristiques du substrat. A partir de l'équation de continuité pour un fluide incompressible dans un état stationnaire, de l'équation (2) et en introduisant le potentiel $\varphi = Km.h$, on obtient:

$$ \nabla^2 \varphi = 0 \qquad (3) $$

L'équation de Laplace est linéaire. Ce qui signifie que le principe de superposition lui est applicable. Il est donc possible de construire des solutions complexes comme superposition de solutions simples. Un écoulement potentiel élémentaire, celui d'une source et d'un puits, à deux dimensions peut être étudié d'après l'équation de continuité en coordonnées polaires. Pour le cas d'une source située à un point F(a,b) associé à un puits situé à un point E(-a,-b) en coordonnées cartésiennes, on obtient [11] :

$$ \psi = \frac{Q}{2\pi} \left( \arctg \frac{y-b}{x-a} - \arctg \frac{y+b}{x+a} \right) \qquad (4) $$

$$ u = \frac{Q}{2\pi} \cdot \left( \frac{x-a}{(x-a)^2 + (x-b)^2} - \frac{(x+a)}{(x+a)^2 + (y+b)^2} \right) \qquad (5) $$

$$ v = \frac{Q}{2\pi} \cdot \left( \frac{y-b}{(x-a)^2 + (y-b)^2} - \frac{y+b}{(x+a)^2 + (y+b)^2} \right) \qquad (6) $$
## 2.3. Application du modèle des puits et des sources à l'étude du substrat

### 2.3.1. Apport à l'aplomb du point de drainage

#### 2.3.1.1. Cas d'une source et d'un puits.

##### i) Formation des lignes de courant

Pour notre cas, la source et le puits sont respectivement situés aux points $(0, \frac{h1}{2})$ et $(0, -\frac{h1}{2})$ dans le plan des coordonnées cartésiennes $(o,x,y)$ ou $h1$ est la hauteur du substrat. En superposant la source et le puits, la fonction de courant s'écrit:

$$ \psi = \frac{Q}{2\pi} \cdot \left( \arctan \frac{y-h1/2}{x} - \arctan \frac{y+h1/2}{x} \right) \quad (7) $$

On constate que le débit n'a aucun effet sur l'allure des lignes de courant. Par raison de symétrie, on les trace sur un carré de côté $\frac{h1}{2}$ représentant les dimensions de la moitié du substrat dans le plan $(o,x,y)$. Les courbes obtenues représentent bien l'allure des lignes de courant (Fig 1.). On constate que la direction des lignes de courant ne respecte pas bien les frontières du substrat de section rectangulaire et ne pourrait s'appliquer qu'à une géométrie de section arrondie bien particulière.

##### ii) Calcul de la vitesse

Dans notre cas les expressions des composantes de la vitesse sont données par :

$$ u = \frac{Q}{2\pi} \cdot \left( \frac{x}{x^2 + (y-h1/2)^2} - \frac{x}{x^2 + (y+h1/2)^2} \right) \quad (8) $$

$$ v = \frac{Q}{2\pi} \cdot \left( \frac{y-h1/2}{x^2 + (y-h1/2)^2} - \frac{y+h1/2}{x^2 + (y+h1/2)^2} \right) \quad (9) $$

En tout point des lignes de courant, la vitesse a la valeur suivante :

$$ V = \sqrt{u^2 + v^2} \quad (10) $$

La vitesse est symétrique par rapport à l'axe des abscisses et elle est beaucoup plus grande à la position de la source et du puits. Elle est d'autant plus faible que la trajectoire est plus longue.

#### 2.3.1.2. Cas d'une somme de sources et de puits

Le modèle basé sur une source et un puits ne satisfait pas les conditions aux limites de l'écoulement à simuler. On est amené à disposer un réseau de sources et de puits répartis périodiquement dans le plan $(x,y)$, obtenus comme images les unes des autres par rapport aux frontières du domaine. L'idée est que le vecteur vitesse obéira aux symétries introduites et suivra le contour du substrat.

##### i) Formation des lignes de courant

Les puits et les sources sont disposés comme il est indiqué dans les configurations représentées en annexe. Plusieurs configurations sont à examiner. Pour le traçage de ces courbes, le principe est le même qu'auparavant, la fonction de courant étant décrite par l'équation suivante :

$$ \psi = \psi_s + \psi_p \quad (11) $$

$\psi_s$ : fonction de courant pour l'ensemble des sources, $\psi_p$ : fonction de courant pour l'ensemble des puits.
$$ \psi_s = \frac{Q}{2\pi} \sum_i \sum_j \operatorname{arctg} \left( \frac{y - y_s[i]}{x - x_s[i]} \right) \quad (12) $$

$x_s[i]$ et $y_s[j]$ représentent respectivement l'abscisse et l'ordonnée de chacune des sources par rapport à l'origine.

<figure><img src="image_1.png" /><figcaption>Fig. 1: Visualisation des lignes de courant cas d'une source et d'un puits</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Fig. 2: Visualisation des lignes de courant cas d'une somme de sources et de puits</figcaption></figure>

Les valeurs maximales atteintes par i et j dépendent de la configuration choisie:

$x_s[i] = (i - m) \cdot l$ pour i allant de 1 à n et $y_s[j] = (-1)^{j+1} \cdot (2 \cdot (j-1) + l) \cdot \frac{h}{2}$ pour j allant de 1 à m.

Le nombre entier m caractérise la configuration choisie : $n = 2m - 1$

Le nombre total N de sources où de puits est donné par : $N = m(2m - 1)$

De même pour la fonction de courant des puits:

$$ \psi_p = -\frac{Q}{2\pi} \sum_i \sum_j \operatorname{arctg} \left( \frac{y - y_p[i]}{x - x_p[i]} \right) \quad (13) $$

$x_p[i]$, abscisse des puits, est ici identique à $x_s[i]$. $y_p[i]$ est l'ordonnée de chaque puits par rapport à l'origine. $y_p[i] = (-1)^j \cdot (2 \cdot (j-1) + l) \cdot \frac{h}{2}$ pour j allant de 1 à m.

Chaque configuration génère ses propres lignes de courant. La configuration est choisie après convergence, le tracé des lignes de courant correspondantes est donné en figure 2. On constate que les lignes de courant épousent bien les contours rectangulaires du substrat.

## ii) Calcul de la vitesse

Pour chaque configuration, les expressions des composantes de la vitesse selon les axes sont données par les équations suivantes :

$$ u = \frac{Q}{2\pi} \sum_j \sum_i \left( \frac{x - x_s[i]}{d_s[i, j]} - \frac{x - x_p[i]}{d_p[i, j]} \right) \quad (14) $$

$$ v = \frac{Q}{2\pi} \sum_j \sum_i \left( \frac{y - y_s[i]}{d_s[i, j]} - \frac{y - y_p[i]}{d_p[i, j]} \right) \quad (15) $$

$$ \text{où } d_s[i, j] = (x - x_s[i])^2 + (y - y_s[j])^2 \quad (16) $$
et

$$d_p[i, j] = (x - x_p[i])^2 + (y - y_p[j])^2 \quad (17)$$

Le module de la vitesse est défini par l'équation (10). La symétrie est toujours gardée par rapport à l'axe des x et la vitesse est maximale au niveau de la position de la source et du puits.

### iii) Formation du bulbe d'irrigation

On se propose dans ce qui suit d'évaluer la géométrie de ce bulbe d'irrigation, le substrat étant supposé saturé. Pour cela on singularise un temps to au cours de l'apport continu de solution sur le substrat, et on suit la solution ainsi repérée qui chemine vers le point de drainage le long des lignes de courant connues. Sur chacune d'elles, on détermine le temps mis par une particule pour se déplacer d'une distance donnée jusqu'au point de drainage. Les incréments successifs du temps sont gérés par l'équation suivante :

$$t_{i+1} = t_i + \frac{\delta}{V} \quad (18)$$

$t_{i+1}$ : temps à la position présente, $t_i$ : temps à la position précédente, $\delta$ : distance entre deux points successifs d'une ligne de courant, $V$ : vitesse moyenne entre les deux positions

$$V = \frac{V_{i+1} + V_i}{2} \quad (19)$$

Une fois obtenue cette matrice des temps de parcours de la solution injectée à to, il suffit de relier les points de même valeur temporelle. L'ensemble des courbes ainsi obtenues permet d'évaluer l'évolution du bulbe d'irrigation au cours du temps (fig 3.), à quelques déformations près qui résultent du manque de précision dans l'allure des lignes de courant et dans la matrice des vitesses. Pour le cas envisagé ici, la solution apportée est de même densité que celle déjà présente dans le substrat, de sorte qu'aucune convection naturelle ne vient se superposer au mouvement forcé du fluide injecté au sein du substrat saturé.

Il apparaît des zones de vitesses faibles, le long de chemins anormalement longs (les angles du substrat), de sorte qu'on peut y prédire un moindre renouvellement de solution au cours du temps. Ce bulbe a été évalué expérimentalement [7] par un volume de solution de colorant apporté par gouttage sur un substrat de laine de roche (fig 4.). Pour avoir une réduction des zones mortes et un renouvellement de solution meilleur qu'auparavant, un deuxième cas est envisagé.

<figure><img src="image_3.png" /><figcaption>Fig. 3: Visualisation du bulbe d'irrigation</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 4: Visualisation expérimentale du bulbe d'irrigation</figcaption></figure>

#### 2.3.2. Points de drainage et d'apport décalés

##### 2.3.2.1.Cas d'une source et d'un puits
i) Formation des lignes de courant

Pour illustrer ce cas, nous considérerons que la source et le puits sont situés aux deux extrémités de la diagonale du substrat. La fonction de courant s'écrit:

$$
y = \frac{Q}{2\pi} \left( \frac{\arctan\left(\frac{y-hl}{2}\right)}{x+\frac{l}{2}} - \frac{\arctan\left(\frac{y+hl}{2}\right)}{x-\frac{l}{2}} \right) \quad (20)
$$

Les courbes obtenues sont représentées en figure 5.

ii) Calcul de la vitesse

Dans ce cas, les composantes de la vitesse selon les axes s'écrivent :

$$
u = \frac{Q}{2\pi} \left( \frac{x+l/2}{(x+l/2)^2 + (y-h1/2)^2} - \frac{x-l/2}{(x-l/2)^2 + (y+h1/2)^2} \right) \quad (21)
$$

$$
v = \frac{Q}{2\pi} \left( \frac{y - hl/2}{(x + l/2)^2 + (y - hl/2)^2} - \frac{y + hl/2}{(x - l/2)^2 + (y + hl/2)^2} \right) \quad (22)
$$

Le module de la vitesse est donné par l'équation (10).

<figure><img src="image_5.png" /><figcaption>Fig. 5: Visualisation des lignes de courant
cas d'une source et d'un puits</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6: Visualisation des lignes de courant
cas d'une somme de sources et de
puits</figcaption></figure>

2.3.2.2. Cas d'une somme de sources et de puits

i) Formation des lignes de courant

Comme le modèle fondé sur une source et un puits ne satisfait pas aux conditions aux limites, on est amené à superposer une somme de puits et de sources disposés de façon à contraindre l'écoulement à suivre les frontières du substrat. Les sources et les puits sont disposés comme indiqué en annexe, plusieurs configurations étant à examiner suivant le degré de précision recherché. Pour chaque configuration la fonction de courant s'écrit :

$$
\psi = \psi_s + \psi_p \tag{23}
$$

$$
\psi_s = \frac{Q}{2\pi} \sum_i \sum_j \arctan\left(\frac{y - y_s[j]}{x - x_s[i]}\right) \quad (24)
$$

$y_s[j]$ est défini comme ci dessus. $x_s[i] = (-1)^i \cdot (2i - 1) \cdot \frac{l}{2}$ pour i allant de 1 à m.

$$
\psi_p = -\frac{Q}{2\pi} \sum_i \sum_j \arctan\left(\frac{y - y_p[j]}{x - x_p[i]}\right) \quad (25)
$$

$y_p[j]$ est défini comme ci dessus et $x_p[i] = (-1)^{i+1} \cdot (2i-1) \cdot \frac{l}{2}$ pour i allant de 1 à m.
m représente la configuration choisie. Les courbes sont représentées sur la figure 6.
## ii) Calcul de la vitesse

Pour chaque configuration, les composantes de la vitesse selon les axes sont données par :

$$ u = \frac{Q}{2\pi} \sum_{i} \sum_{j} \left( \frac{x - x_s[i]}{(x - x_s[i])^2 + (y - y_s[j])^2} - \frac{x - x_p[i]}{(x - x_p[i])^2 + (y - y_p[j])^2} \right) \quad (26) $$

$$ v = \frac{Q}{2\pi} \sum_{i} \sum_{j} \left( \frac{y - y_s[i]}{(x - x_s[i])^2 + (y - y_s[j])^2} - \frac{y - y_p[i]}{(x - x_p[i])^2 + (y - y_p[j])^2} \right) \quad (27) $$

Le module de la vitesse est toujours donné par l'équation (16)

## iii) Formation du bulbe d'irrigation

Pour l'évaluation du bulbe d'irrigation, on suit la même procédure. Ce bulbe d'irrigation est représenté en figure 7. On remarque que la solution irrigue de façon plus homogène le substrat, ce qui diminue les zones mortes. Ce bulbe a été évalué experimentalement [7] par usage de colorant (fig 8.). En pratique l'étude de la concentration se fait au drainage, elle représente une valeur moyenne d'une couche de fluide au fond du substrat décalé dans le temps. Il serait utile d'avoir une idée sur la variation de la concentration de la solution en chaque point du substrat, ceci est difficile à réaliser expérimentalement. A ce niveau, une approche théorique est plus facile.

<figure><img src="image_7.png" /><figcaption>Fig. 7: Visualisation du bulbe d'irrigation</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 8: Visualisation expérimentale du bulbe d'irrigation</figcaption></figure>

Le modèle précédent, dit des puits et des sources, est trop sommaire pour nous renseigner sur l'évolution de la concentration de la solution dans le substrat, puisqu'il est à l'image de l'écoulement potentiel d'un liquide non visqueux monophasique. On est donc amené à étudier la variation de la concentration en formulant le problème d'une façon plus réaliste qui tienne compte d'effets visqueux et diffusifs.

## 2.4. Etude de la dispersion dans le substrat

La dispersion du soluté au sein d'une solution saline s'écoulant à travers un milieu poreux est un processus qui fait intervenir un couplage étroit entre le mouvement de transport au sein des espaces libres, certainement complexe à échelle locale, et le phénomène de diffusion moléculaire qui joue un rôle notable sur ces petites distances. Cette dispersion est un mécanisme qui peut être simulé globalement par l'effet d'un coefficient de diffusion effectif très amplifié dans l'équation d'évolution de la concentration. Plusieurs approches permettent de décrire ce phénomène par des modèles simples, et certaines peuvent même donner lieu à des solutions analytiques. Un tube capillaire, un ensemble de tubes capillaires, un empilement de particules, etc. sont des exemples de tels modèles [9].
On se base sur les équations de Navier-Stockes, capables de décrire le mouvement d'un fluide très visqueux et de concentration variable à travers le volume parallèlépipédique du substrat, ce qui donne des lignes de courant et à un champ de vitesses. On complète la description par une équation d'évolution de la concentration, propriété scalaire soumise à un effet de transport et de diffusion, ce qui donne une carte des concentrations dans le substrat. On perd le caractère intuitif du formalisme précédent

## 2.4.1. Mise en équation et conditions aux limites

Le fluide, newtonien, étant supposé incompressible au sens strict, les seules inconnues sont le champ de vitesse, la pression et la concentration. On restreint, par symétrie, l'écoulement au plan de repère (x,o,y) :

### 2.4.1.1. Champ de vitesses

On se donne la distribution de la vitesse dans un premier temps; celle-ci est injectée dans l'équation de la concentration, ce qui nous permet de déterminer sa variation. Dans le cas présent d'un fluide homogène, la variation de la concentration n'a pas d'effet sur la vitesse. L'écoulement considéré étant stationnaire dans sa dynamique, les équations de "Navier-Stokes" du fluide incompressible s'écrivent :

#### i) Equations de mouvement

$$u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = - \frac{1}{\rho} \frac{\partial p}{\partial x} + v \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \quad (28)$$

$$u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = - \frac{1}{\rho} \frac{\partial p}{\partial y} + v \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right) + g \quad (29)$$

#### ii) Equation de continuité

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \quad (30)$$

### 2.4.1.2. Equation de la variation de la concentration

Ici, il n'est plus vrai que les conditions sont stationnaires, puisque l'on veut considérer le cas du remplacement d'une solution par une autre de composition différente. Pour le cas d'un milieu poreux où la porosité provoque un brassage local qui favorise le travail de la diffusion, il a été vérifié [9] que tout se passe comme si le coefficient de diffusion était remplacé par un coefficient de dispersion K, beaucoup plus grand. On obtient l'équation de dispersion :

$$\frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = K \left( \frac{\partial^2 C}{\partial x^2} + \frac{\partial^2 C}{\partial y^2} \right) \quad (31)$$

Dans le cas d'un milieu poreux granulaire, le fluide injecté et le fluide résident se mélangent selon un coefficient de dispersion qui est la somme de deux composantes l'une liée à la phase stationnaire (diffusion) et l'autre liée à la phase mobile (dispersion) [3]. La solution qui progresse dans un milieu fibreux est ralentie au contact des fibres (condition de non glissement aux parois), le cisaillement du fluide dans la couche limite crée une dispersion latérale à forts gradients de concentration vers les parois. Dans cette couche limite à fort gradient, la diffusion moléculaire intervient pour niveler les concentrations. Ce phénomène explique l'accroissement du volume apparent du bulbe de solution fraîche dans un substrat fibreux ainsi que l'étalement du front de concentration observé expérimentalement à l'aide des colorants. Au total, on peut supposer que l'on a affaire dans un milieu fibreux à un coefficient de dispersion analogue à celui introduit pour les milieux granulaires.
### 2.4.2. Choix des conditions aux limites et de la condition initiale

La détermination de la vitesse et de la concentration nécessite la connaissance des conditions aux limites et de la condition initiale.

On suppose qu'on a une vitesse d'injection de la solution non nulle uniquement au milieu de la surface supérieure du substrat (zone d'apport de la solution) et un drainage dans les mêmes conditions au milieu de sa face inférieure (fente de drainage). Pour éviter la discontinuité de la vitesse qui résulterait d'une zone d'apport rigoureusement ponctuelle, ce qui a un effet déstabilisateur sur les calculs, on représente le profil de la vitesse dans le jet d'apport par une gaussienne (fig 9).

<figure><img src="image_9.png" /><figcaption>Fig. 9: Profil de la vitesse injectée</figcaption></figure>

Cette modification est justifiée car l'apport de la solution prend place, par effet capillaire, sur une zone d'une certaine extension à la surface supérieure du substrat, et avec une vitesse initiale faible.

$$x=0, y=\frac{hl}{2}, U=-Vi.\exp\left(-\lambda\left(x-\frac{l}{2}\right)^2\right) \qquad (32)$$

$\lambda$ est une constante, son choix dépend du débit de la solution injectée. A partir d'un débit $Q$ en $(\text{m}^3\text{s}^{-1})$ connu, on peut calculer $\lambda$ par : $Q = \int_{-\infty}^{+\infty} Vi.\exp(-\lambda(x - \frac{1}{2})^2) dx$.

La vitesse est maximale au milieu et tend vers zéro latéralement ; ce profil a une influence sur le flux de matière apporté sur la face supérieure. La concentration à la frontière est soumise à une condition aux limites particulière. Le flux d'apport est évacué par diffusion et par transport, ce qui s'écrit :

$$x=0, y=\frac{hl}{2}, \quad UC1=UC-K\frac{\partial C}{\partial y}. \qquad (33)$$

où $C1$ est la concentration de la solution d'apport et $U$ est la vitesse d'apport.

A la partie basse du substrat, on suppose qu'on a le même profil de vitesse. Pour la concentration, la condition est la continuité du flux massique. Comme le domaine extérieur est l'air (où le film d'emballage du substrat), il n'y a pas de diffusion vers le côté externe. Après égalisation des expressions du flux de masse à cette surface côté intérieur et celui du côté extérieur, on arrive à écrire la condition sous la forme suivante:

$$x=0, y=-\frac{hl}{2}, \quad \left(\frac{\partial C}{\partial y}\right)=0 \qquad (34)$$

Pour les parois latérales du substrat, on suppose qu'elles sont imperméables aux échanges massiques et que le fluide satisfait aux conditions d'adhérence à ces parois. A l'instant initial, on suppose que la concentration $C$ a une valeur $C_0$
## 2.5. Formulation des équations en fonction de la vorticité et de la fonction de courant

Dans le cas d'un écoulement plan, l'équation de continuité est prise en compte en définissant une fonction de courant définie par l'équation (4). Dans le cas plan, le vecteur tourbillon n'a qu'une composante $\omega$ sur la normale au plan [11].

$$ \omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \qquad (35) $$

La difficulté liée à la détermination de la pression conduit à chercher à éliminer celle-ci des équations. C'est ce qui est fait en appliquant l'opérateur rotationnel à l'équation du mouvement qui contient le gradient de la pression. Des équations (35) et (34), on obtient :

$$ u \frac{\partial \omega}{\partial x} + v \frac{\partial \omega}{\partial y} = v \left( \frac{\partial^2 \omega}{\partial x^2} + \frac{\partial^2 \omega}{\partial y^2} \right) \qquad (36) $$

$$ \omega = - \left( \frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} \right) \qquad (37) $$

## 2.6. Réduction des variables :

On définit les variables réduites :

$$ \begin{aligned} x^* &= \frac{x}{l}; \quad \psi^* = \frac{\psi}{Vi.l}; \quad v^* = \frac{v}{Vi}; \quad t^* = \frac{Vi}{l}; \\ y^* &= \frac{y}{l}; \quad \omega^* = \frac{l}{Vi}\omega; \quad u^* = \frac{u}{Vi}; \quad C^* = \frac{C}{CI}. \end{aligned} \qquad (38) $$

Des équations (42), (43) et (37), on obtient les équations suivantes :

$$ u^* \frac{\partial \omega^*}{\partial x^*} + v^* \frac{\partial \omega^*}{\partial y^*} = \frac{1}{\text{Re}} \left( \frac{\partial^2 \omega^*}{\partial x^{*2}} + \frac{\partial^2 \omega^*}{\partial y^{*2}} \right) \qquad (39) $$

$$ \omega^* = -\nabla^2 \psi^* \qquad (40) $$

$$ \frac{\partial \mathcal{C}^*}{\partial t^*} + u^* \frac{\partial \mathcal{C}^*}{\partial x^*} + v^* \frac{\partial \mathcal{C}^*}{\partial y^*} = \frac{1}{\text{Re} \cdot \text{Sc}} \left( \frac{\partial^2 \mathcal{C}^*}{\partial x^{*2}} + \frac{\partial^2 \mathcal{C}^*}{\partial y^{*2}} \right) \qquad (41) $$

avec $Re = \frac{Vi.L}{\nu}$ (nombre de Reynolds) et $Sc = \frac{\nu}{K}$ (nombre de Schmidt). On obtient aussi les grandeurs suivantes :

$$ x^* = \pm \frac{1}{2}; \quad \forall y^*, u^* = v^* = 0, \quad \frac{\partial \mathcal{C}^*}{\partial x^*} = 0. \qquad (42) $$

$$ y^* = \frac{h_1}{2l}; \quad \forall x^*, u^* = 0; v^* = -\exp\left(-\lambda\left(x^* - \frac{1}{2}\right)^2\right); \quad \exp\left(-\lambda\left(x^* - \frac{1}{2}\right)^2\right) = \exp\left(-\lambda\left(x^* - \frac{1}{2}\right)^2\right) C^* \frac{K}{Vi} \frac{\partial \mathcal{C}^*}{\partial y^*} \qquad (43) $$

$$ y^* = -\frac{h_1}{2l}; \quad \forall x^*, u^* = 0; v^* = -\exp\left(-\lambda\left(x^* - \frac{1}{2}\right)^2\right); \quad \frac{\partial \mathcal{C}^*}{\partial y^*} = 0. \qquad (44) $$

# 3. ANALYSE NUMERIQUE

Les conditions aux limites étant indépendantes du temps, la solution pour $\omega$ l'est également.
### 3.1. Méthode pseudo-instationnaire:

On résout le système suivant :

$$ \frac{\partial \omega^*}{\partial t^*} + u^* \frac{\partial \omega^*}{\partial \hat{x}^*} + v^* \frac{\partial \omega^*}{\partial \hat{y}^*} = \frac{1}{\text{Re}} \left( \frac{\partial^2 \omega^*}{\partial \hat{x}^{*2}} + \frac{\partial^2 \omega^*}{\partial \hat{y}^{*2}} \right) \quad (45) $$

$$ \frac{\partial \psi^*}{\partial t^*} - \alpha (\nabla^2 \psi^* + \omega^*) = 0 \quad (46) $$

En faisant tendre t vers l'infini, on obtient la solution stationnaire. $\alpha$ est une constante positive. Pour éviter les phénomènes d'instabilité induits par les schémas explicites, on choisit un schéma implicite. Pour la résolution de ce système, on choisit la méthode implicite des directions alternées. (A.D.I) [10]. La solution stationnaire doit satisfaire les conditions aux limites précédentes. Il faut compléter le problème par le choix des conditions initiales.

Comme c'est un cas stationnaire, le choix des conditions initiales est arbitraire:

$u = v = \varepsilon$ ; $\varepsilon$ est une petite valeur pour amorcer le mouvement, $\omega = \psi = 0$

Pour déterminer la fonction de courant et la vorticité sur les limites du substrat, on suit la démarche faite par D.Euvrard [4]. Une fois satisfaite cette condition, on passe au calcul de la concentration, en utilisant à nouveau la méthode A.D.I.

### 3.2. Visualisation des courbes

Les courbes ainsi obtenues sont représentées dans les graphiques ci-dessous. Les lignes de courant sont montrées en figure 10. On remarque que ces lignes suivent avec précision les frontières du domaine rectangulaire. Le remplacement progressif de la solution ancienne est illustré en figure 11. Cette courbe met en évidence la superposition de deux phénomènes :

* dans un premier temps, on observe que la concentration est constante (et égale à $C_o$) qui correspond au phénomène de chasse d'un volume donné de liquide initial.

* dans un second temps, on observe la montée rapide de la concentration sous la forme d'une courbe d'une fonction erfc. Cette montée est due à un effet de chasse partielle du liquide combiné à un phénomène de dispersion hydrodynamique.

<figure><img src="image_10.png" /><figcaption>Fig. 10: Visualisation des lignes de courant</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 11: Variation de la concentration en fonction du temps à un noeud situé au milieu du substrat</figcaption></figure>

La solution injectée chasse la solution ancienne comme le ferait un piston ajouré, dont l'effet se trouve atténué par le processus de dispersion. L'évolution de la concentration suit ensuite une progression asymptotique lente qui résulte du mélange progressif de deux flux de liquides différents. On peut schématiser en disant qu'il existe deux phases, l'une mobile et l'autre fixe, sachant que les échanges massiques se font par diffusion de l'une vers l'autre.

L'évolution de la concentration à la base du substrat est représentée ci-dessous (fig 12.). Pour ce qui est du débit d'apport Q nécessaire pour renouveler le volume de la solution qui
existe déjà dans le substrat, on peut le calculer de la façon suivante. On détermine le temps t nécessaire pour que la concentration à la base du substrat atteigne la valeur de celle de la solution d'apport à (par exemple) 10% près; à ce moment, on a : $Q = U.C1.t (mole.m⁻²)$

<figure><img src="image_12.png" /><figcaption>Fig. 12: Variation de la concentration à la base du substrat à un temps fixé</figcaption></figure>

Lorsqu'on fait un apport de solution sur un substrat saturé, le fluide injecté se propage plus rapidement le long des lignes de courant centrales, ce qui provoque l'apparition d'un phénomène de dispersion latérale. Le mélange entre le fluide initial et le fluide injecté est le résultat d'un compromis entre le transport et la dispersion. Le bulbe d'irrigation obtenu par le calcul est représenté ci dessous (fig 13.).

<figure><img src="image_13.png" /><figcaption>Fig. 13: Visualisation du bulbe d'irrigation</figcaption></figure>

Ce bulbe est évalué visuellement par l'extension du volume coloré lors d'un apport de colorant sur le substrat de laine de roche. On trouve ainsi que le bulbe a un volume apparent double de celui de l'apport pratiqué, ce qui confirme l'importance du phénomène de dispersion par les fibres du substrat (la part de solution "adhérente" au fibres est voisine de 50%).

## 4. DISCUSSION

Le modèle des puits et des sources donne une description macroscopique de ce qui se passe dans le substrat. Il nous permet également de voir l'allure des lignes de courant et l'hétérogénéité de la vitesse ainsi que la forme du bulbe d'irrigation. Le résultat de la modélisation reproduit assez bien ce qui a été observé lors de visualisations par colorants.

Dans le cas où la source et le puits sont à l'aplomb, on a des zones mortes aux coins et aux côtés latéraux du substrat. Par contre dans le cas décalé, les zones mortes sont réduites, la solution d'apport intéresse alors un plus grand volume du substrat, ce qui permet d'obtenir un meilleur renouvellement de la solution ancienne. Quand on alimente un substrat par une solution de concentration différente de celle qu'il contient, la solution qui pénètre déplace la solution initiale loin en aval de l'interface puisque le liquide est incompressible et agit comme un piston agissait en amont. Cela provoque au bas du substrat un drainage de solution
ancienne exclusivement. Si l'on examine plus en détail ce qui se passe à l'interface, le piston
doit être compris comme étant ajouré (à échelle fine) et un mélange des deux solutions est
provoqué par la dispersion bien au delà du volume d'encombrement strict de l'apport de
solution.

Afin de minimiser les volumes morts mal renouvelés et de provoquer un bon brassage du substrat, il suffit placer les goutteurs et les fentes de drainage à distance relative éloignée (en diagonale), car les lignes de courant seront de longueur plus uniforme et la dispersion sera plus efficace pour atténuer les écarts de concentration sur les trajets différents.

L'étude de la concentration au drainage ne permet pas de connaître sa valeur à chaque
point du substrat, elle nous donne qu'une valeur moyenne dans une couche de fluide proche
du fond du substrat, qui représente la physionomie du substrat avec un certain retard ( le délai
mis par le front de solution fraîche pour progresser à travers le milieu). Le modèle présenté
doit permettre de connaître la variation instantanée de la concentration en chaque point du
substrat.

## NOMENCLATURE

<table><tr><td>α</td><td>constante</td><td>K</td><td>coefficient de dispersion</td><td>[m<sup>2</sup>.s<sup>-1</sup>]</td></tr><tr><td>ε</td><td>erreur admise</td><td>k</td><td>pas de descrétisation</td><td>[s]</td></tr><tr><td>φ</td><td>potentiel des vitesses</td><td>km</td><td>conductivité hydraulique</td><td>[m<sup>2</sup>]</td></tr><tr><td>γ</td><td>poids spécifique</td><td>l</td><td>longueur du substrat</td><td>[m]</td></tr><tr><td>λ</td><td>constante</td><td>P</td><td>pression totale</td><td>[Pa]</td></tr><tr><td>ν</td><td>Viscosité cinématique</td><td>Q</td><td>débit</td><td>[m<sup>3</sup>.s<sup>-1</sup>]</td></tr><tr><td>ρ</td><td>masse volumique</td><td>Re</td><td>nombre de Reynolds</td><td></td></tr><tr><td>δ</td><td>distance entre deux points</td><td>Sc</td><td>nombre de Schmidt</td><td></td></tr><tr><td>ω</td><td>composante du vecteur tourbillon</td><td>t</td><td>Temps</td><td>[s]</td></tr><tr><td>ψ</td><td>fonction de courant</td><td>U</td><td>vitesse d'injection</td><td>[m.s<sup>-1</sup>]</td></tr><tr><td>Ω̄</td><td>Vorticité</td><td>u, v</td><td>composantes de la vitesse</td><td>"</td></tr><tr><td>C</td><td>concentration</td><td>V</td><td>module de la vitesse</td><td>"</td></tr><tr><td>Co</td><td>concentration initiale</td><td>Vi</td><td>vitesse d'injection maximale</td><td>"</td></tr><tr><td>C1</td><td>Concentration de la solution apportée</td><td>x,y</td><td>coordonnées cartésiennes</td><td></td></tr><tr><td>Dt</td><td>pas de temps</td><td>z</td><td>Elévation</td><td>[m]</td></tr><tr><td>H</td><td>perte de charge</td><td>i, j</td><td>incréments</td><td></td></tr><tr><td>H</td><td>pas d'espace de descrétisation</td><td>s</td><td>Source</td><td></td></tr><tr><td>H1</td><td>hauteur du substrat</td><td>p</td><td>Puits</td><td></td></tr></table>

## REFERENCES

[1] J.R. Allen and Ditsworth, "Fluid Mechanics", International Student Edition, pp. 215-252, 1972.

[2] S. Bougoul; 'Etude d'écoulements salins dans des milieux fibreux saturés utilisés en substrats de culture', Thèse de Doctorat en Sciences; Université de Nice. France, 1996.

[3] De Wiest, "Flow through porous media", Academic press New York London, pp. 1-52, 1969.

[4] D. Euvrard, "Résolution numérique des équations aux dérivées partielles", Masson Paris Milan Barcelone, 3e édition, pp. 149-175, 1994.

[5] R. Gras, "Propriétés physiques du substrat", Les cultures hors sol, Ouvrage collectif dirigé par Denis Blanc, pp 78-126, 1987.

[6] Incropera, P. Frank, De Witt and P. David, "Fundamentals of heat and mass transfer", Third edition, John Willey & Sons, pp 872-886, 1990

[7] O. Kiffer, "Réponse hydraulique d'un substrat de laine de roche de type "Cultilène"", Mémoire de DEA d'Energétique, Nice. France, 1992.
[8] A. Morisot, "Déplacement avec mélange des nitrates dans un matériau de porosité bimodale", Ann. agron., Vol. 31 (1), pp.1-17, 1980.

[9] J. B. Passioura, 'Hydrodynamic Dispersion in aggregated media, 1. Theory', Soil science, vol 111, pp. 339-344, 1971.

[10] R. Peyret, Taylor and D. Thomas, 'Computational Methods for Fluid Flow', Springer-Verlag New York Heidelberg Berlin, 1983.

[11] I. L. Ryming, "Dynamique des fluides", Un cours de base du deuxième cycle universitaire, Presses polytechniques romandes, 1985.

# ANNEXE

Les sources et les puits sont disposés de la façon suivante pour chaque cas :

1- Cas de l'apport à l'aplomb :

Chaque configuration à ses propres sources et puits :

Configuration 1:

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

Configuration 2:

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

Configuration 3:

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

2- Cas décalé :

Configuration 1:

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

Configuration 2:

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

Configuration 3:

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