Modélisation numérique du rejet des saumures en mer

Modélisation numérique du rejet des saumures en mer

L'objectif de cette partie est d'étudier le rejet dans l'Océan Atlantique de ces eaux de saumures. Il s'agit ici de s'assurer que la salinité de l'eau autour de la buse de rejet n'augmente pas de manière trop importante et que la dilution se fait de manière efficace.

Nous réaliserons tout d'abord une étude théorique de ce rejet afin d'en déterminer les caractéristiques physiques. Puis nous réaliserons une étude numérique avec le logiciel de modélisation STAR-CCM+ en s'appuyant sur notre étude théorique.

Cadre de l'étude et position du problème

Cadre de l'étude et position du problème

Le dimensionnement de l'usine montre qu'elle devrait rejeter $40824m³/jours$ d'une salinité de $71,7g/l$. L'objectif de cette partie est d'étudier le rejet en mer de ces eaux de saumures. La configuration du rejet est représentée sur la figure 1:

Figure 1 : Schéma du mode de rejet

De plus nous ferons les hypothèses suivantes pour notre étude:

Tableau 1 : Hypohèses

Débit du rejet $40824m³/jour$
Vitesse du rejet $5m.s^{-1}$
Profondeur du rejet $15m$
Température du rejet $15°C$
Salinité du rejet $71,7 g.l^{-1}$
Température de l'eau au fond $15°C$
Température de l'eau en surface $20°C$
Salinité de l'eau au fond $35g.l^{-1}$
Salinité de l'eau en surface $30g.l^{-1}$

Ces hypothèses nous permettent de calculer les masses volumiques des différents fluides en fonction de leur salinité, leur température et des conditions de pression dans lesquels ils se trouvent. Pour cela on utilise la méthode ies80 de l'UNESCO qui se présente sous la forme d'un script Matlab (voir en annexe) : $\rho=ies80(s,t,p)$.

On obtient ainsi les masses volumiques suivantes:

$\rho_{fond}=1026 kg.m^{3}$, $\rho_{surface}=1021kg.m^{3}$ et $\rho_{rejet}=1055kg.m^{3}$

Etude théorique

Etude théorique

A présent que nous avons défini le cadre de notre d'étude, nous pouvons réaliser une étude théorique de la physique du rejet. Dans un premier temps on considère que l'ensemble du rejet s'effectue au travers d'un seul injecteur. L'un des objectifs de cette partie sera de déterminer la hauteur maximale du rejet afin d'adapter notre géométrie pour la modélisation numérique.

Calcul du diamètre de l'injecteur

$Q$ le débit d'éjection des saumure est défini de la manière suivante:

$Q=\dfrac{\pi D^2 V}{4}$

Où $D$ est le diamètre de l'injecteur et $V$ la vitesse d'éjection.

Le débit étant imposé, $D=40824~m^{3}/j$, on choisit une vitesse d'éjection $V=5~m.s^{-1}$ permettant une bonne dilution. On en déduit ainsi le diamètre de l'injecteur : $D=35~cm$

 Calcul des flux spécifiques

  • Flux spécifique de quantité de mouvement :

$M=\int v(v.n)\, \mathrm dS$

  • Flux spécifique de flottabilité :

$F=\int v.g'\, \mathrm dS$

Où g' est la gravité réduite définie par : $ g'=\dfrac{\rho_{0}-\rho_{a}}{\rho_{\infty}}g $

Avec:

$\rho_{0}$ la masse volumique de l'eau de mer.

$\rho_{r}$ la masse volumique  du rejet.

$\rho_{\infty}$ la masse volumique de référence.

$g$ l'accélération de la pesanteur.

Notre problème étant un problème stationnaire et, de plus, la vitesse $V$ étant normale à la surface de rejet on peut écrire que:

$M=Q.V$

$F=V.g'$

Ainsi, $M=2,36~m^{4}.s^{-2}$ et $F=-0,13~m^{4}.s^{-3}$

Influence de la stratification

On désire à présent évaluer l'influence de la stratification de l'océan sur notre rejet. Pour cela on calcul la fréquence de Brünt-Vünsala en faisant l'hypothèse d'une stratification linéaire :

$N^2=-\dfrac{g}{\rho_{0}}\dfrac{\partial \rho}{\partial z}$

Ainsi, $N=0.057~s^{-1}$

On calcule alors le nombre sans dimension $\dfrac{MN}{B}$ :

$\dfrac{MN}{B}=-1.05$

Calcul du nombre de Richardson initial

Le nombre de Richardson est un nombre adimensionnel qui compare l'énergie cinétique d'un fluide avec son énergie potentielle gravitationnelle. Les longueurs $l_{Q}$ et $l_{M}$ permettent de caractériser les régions où notre rejet à un caractère de jet idéal ou de panache idéal :

$l_{Q}=\dfrac{Q}{M^{\frac{1}{2}}}$

$l_{M}=\dfrac{M^{\frac{3}{4}}}{F^{\frac{1}{2}}}$

Ainsi, $l_{Q}=0.31~m$ et $l_{M}=5.3~m$

On en déduit que :

  • Pour $0 \leq h \leq l_{Q}$ l'écoulement est contrôlé par la géométrie de la buse d'éjection. On remarque que $l_{Q} \approx D$.
  • Pour $l_{Q} \leq h \leq l_{M}$ l'écoulement est proche du jet idéal.
  • Pour $l_{M} \leq h \leq 15~m$ l'écoulement est proche du panache idéal.

$Ri_{0}=\dfrac{l_{Q}}{l_{M}}=\dfrac{Q.F^{\frac{1}{2}}}{M^{\frac{5}{4}}}$

  D'où, numériquement : $Ri_{0}= 0,058$

L'écoulement est donc initialement contrôlé par les effets d'inertie $Ri_{0} \leq 1$.

Modélisations numériques avec Star-CCM+

Après avoir compris les mécanismes physiques qui régissent le rejet nous allons réaliser une étude numérique à l'aide du logiciel STAR-CCM+.

Le but est ici de s'assurer que dans les conditions de rejets que nous avons établis la concentration en saumure dans le domaine ne devienne pas trop importante et ne vienne pas polluer l'environnement.

Préambule

Présentation de STAR-CCM+

STAR-CCM+ est un logiciel de modélisation numérique pour la mécanique des fluides développé par CD-adapco depuis 2004. CCM signifie "Computational Continuum Mechanics". Contrairement à d'autres logiciel de modélisation pour la mécanique des fluides, STAR-CCM+ permet de réaliser l'ensemble d'une étude numérique depuis la création de la géométrie et du maillage jusqu'au post-traitement. De plus, il permet de résoudre simultanément les équations de transport de matière et de chaleur ce qui augmente la précision tout en réduisant les temps de calcul.

Géométrie et maillage

On choisit de réaliser une géométrie simple pour étudier le rejet de saumure : un cube de $15~m$ de côté au centre duquel se trouve un cylindre "fictif" de $3~m$ de diamètre et de $12~m$ de hauteur. Ce cylindre nous permettra de raffiner le maillage dans la partie de la géométrie qui nous intéresse le plus, là où se trouveront les vitesses et les concentrations les plus importantes. Il ne s'agit pas d'un cylindre physique. Une surface circulaire de $35~cm$ de diamètre est placée à la base du cylindre, elle modélisera la buse d'éjection des saumures.

On maille à présent notre géométrie (voir Figure 2 et Figure 3). On choisit le modèle "Polyhedral Mesher" pour sa précision et sa facilité de construction. On ajoute "Surface Remesher" afin d'améliorer la qualité du maillage à proximité des parois et d'optimiser le maillage volumique polyhédral. Le cube et le cylindre sont maillés séparément afin d'avoir une résolution supérieur dans le cylindre. La taille de base d'une maille étant de $30~cm$ dans le cube et de $6~cm$ dans le cylindre.

Figure 2 : Vu en coupe du maillage

Figure 3 : Vue du maillage en trois dimensions

On obtient ainsi  $144091$ cellules dans le cylindre et dans $681559$ le cube, soit $825650$ cellules au total.

Modèles physiques

On choisit de réaliser une modélisation en trois dimensions en régime turbulent permanent. On choisit d'utiliser le modèle "Multi-Component Liquid" avec l'option "Non-reacting" afin de pouvoir caractériser les propriétés de nos deux fluides : l'eau de mer et les saumures. Le modèle de résolution de la turbulence choisi est le modèle "K-Epsilon" (de type Reynolds-Average Navier-Stokes). Ce modèle est choisi pour sa robustesse, son économie et sa relative précision. On ajoute en outre l'option "Realizable K-Epsilon Two-Layer", qui combine le modèle "Realizable K-Epsilon" avec l'approche "two layers" qui permet d'appliquer le modèle "K-Epsilon" dans une sous-couche visqueuse. Ainsi, proche d'une paroi, la valeur du coefficient de diffusion turbulente $\epsilon$ et la viscosité turbulente $\mu_{f}$ sont des fonctions de la distance à la paroi.

On impose en outre l'accélération de la pesanteur $g=-9,81~m.s^{-2}$ selon la verticale.

Simulation de référence

Étude du champ de vitesse

On réalise en premier lieu une simulation simple. L'injecteur est modélisé par la condition à la limite "Velocity Inlet" avec une valeur de $5~m.s^{-1}$ et une masse volumique de $\rho_{rejet}=1055~kg.m^{3}$ , le fond de l'océan est modélisé par une frontière de type "Wall" et on impose aux autres faces du cubes des frontières de type "Pressure Outlet" avec la même pression de référence $P=1013~hPa$. L'océan est initialisé avec une vitesse nulle et une masse volumique de $\rho_{océan}=1020~kg.m^{3}$. On ne tient pas compte ici des effets de la stratification.

Cette simulation servira de référence au reste de l'étude.

Le résultat obtenu pour le champ de vitesse est le suivant :

Figure 4 : Champ de vitesse

On va comparer les résultats obtenus par simulation numérique (voir Figure 4) avec ceux donnés par la théorie. Comme nous l'avons dans la partie théorie des jets, notre écoulement va d'abord avoir un comportement proche du jet idéal, puis à partir d'une certaine hauteur, un comportement proche du panache idéal.

On étudie tout d'abord l'allure générale du rejet. Pour cela on s'intéresse à la grandeur $b_{w}$ qui caractérise la largeur du rejet. Théoriquement, on a :

$b_{w}=a_{0}z$ dans le cas du jet idéal avec $a_{0}=0.1$

$b_{w}=b_{0}z$ dans le cas du panache idéal avec $b_{0}=0.107$

Figure 5 : Comparaison de la valeur de $b_{w}(z)$

Sur la figure ci-dessus on constate que la théorie est la simulation sont relativement proches (voir Figure 5). L'écart observé provient notamment du fait que le résultat théorique est obtenu en faisant l'hypothèse d'une source ponctuelle alors que notre rejet s'effectue à partir d'une buse de $35~cm $ de diamètre.

On s'intéresse à présent aux profils de vitesses transverses.

Dans le cas du jet idéal, la théorie nous donne la relation suivante :

$w(r,z)=w_{m}(z).exp\{\dfrac{-r^{2}}{b_{w}^{2}z^{2}}\}$

Avec, $a_{1}=7$ une constante, $w_{m}(z)$ la vitesse sur l'axe du rejet donné par :

$w_{m}(z)=a_{1}\dfrac{\sqrt{M}}{z}$

Dans le cas du panache idéal, la théorie nous donne la relation suivante :

$w(r,z)=w_{m}(z).exp\{\dfrac{-r^{2}}{b_{w}^{2}z^{2}}\}$

Avec, $b_{1}=4,7$ une constante, $w_{m}(z)$ la vitesse sur l'axe du rejet donné par :

$w_{m}(z)=b_{1}(\dfrac{F}{z})^{\frac{1}{3}}$

On compare les profils de vitesses transverses théorique à $2~m$ de hauteur et à $14,5~m$ de hauteur avec les résultats de la simulation :

Figure 6 : Profils de vitesses transverses à $2~m$

À $2~m$ de hauteur les deux profils sont similaires (voir Figure 6). On note cependant que la vitesse théorique est légérement suprérieure à la vitesse obtenue par simulation. Elle est même supérieure à la vitesse initiale d'éjection. Ceci est dû au domaine de validité de la formule théorique utilisée, qui n'est valable que pour des valeurs assez grandes de z. On estime toutefois ici que l'on est à la limite du domaine de validité de la formule.

Figure 7 : Profil de vitesse à $14,5~m$

À $14,5~m$ le rejet a un comportement proche de celui du panache idéal (voir Figure 7). La théorie est cette fois-ci encore en accord avec les résultats de la simulation. On note cependant que la vitesse théorique est plus importante que la vitesse obtenue numériquement.

La figure suivante compare la décroissance de la vitesse le long de l'axe. On obtient des résultats assez moyens entre $5~m$ et $10~m$ d'altitude (voir Figure 8).

Figure 8 : Décroissance de la vitesse axiale

Etude de la concentration en saumure

Le profil de concentration en saumure obtenu avec simulation précédante est le suivant :

Figure 9 : Champ de concentration en saumure

On s'intéresse au profil de concentration en saumure le long du rejet afin d'étudier la dilution en fonction de la hauteur (voir Figure 9). On trace trois profils de concentration, l'un à $1~m$, un autre à $7~m$ et un dernier à $14,5~m$ au-dessus de la buse de rejet (voir Figure 10).

Figure 10 : Profils transverse de concentration en saumure

À $1~m$ au dessus de la buse la concentration en saumure est très forte.

À $7~m$ au dessus de la buse la concentration est presque divisée par deux, le profil gaussien s'applatit.

Enfin, en sortie (à $14,5~m$) le profil continu à s'applatir et la dilution est relativement bonne : on a divisé la concentration en saumure par dix.

Étude de l'influence d'un courant transverse

Étude de l'influence d'un courant transverse

On modifie la condition à la limite de l'une des faces du cube. À la place de "Pressure Outlet" on choisi "Velocity Inlet" afin de modéliser un courant marin.

On réalise des simulations avec des courants d'une vitesse différentes. Ceci se traduit par une condition à la limite différente sur l'une des faces du cube (sur le côté gauche).

Avec un courant de $0,1~m.s^{-1}$ on obtient le résultat suivant :

Figure 11 : Champ de concentration en saumure avec  courant transverse de $0.1~m.s^{-1}$

Avec une vitesse de $0,1~m.s^{-1}$ la trajectoire du rejet est peu dévié par rapport à la simulation de référence, cependant la dilution augmente de manière significative (voir Figure 11).

Pour un courant de $0,5~m.s^{-1}$ on obtient le résultat suivant :

Figure 12 : Champ de concentration en saumure avec  courant transverse de $0.5~m.s^{-1}$

Avec un courant ayant une vitesse de $0,5~m.s^{-1}$ le rejet est cette fois dévié de manière importante puisque la "sortie" à cette fois lieu sur le côté du domaine et plus sur le haut du domaine (voir Figure 12). En outre la dilution des saumure ce fait très bien puisque à $7~m$ de l'entrée du rejet on a dilué environ dix fois notre rejet.

Pour un courant de $1~m.s^{-1}$ on obtient le résultat suivant :

Figure 13 : Champ de concentration en saumure avec  courant transverse de $1~m.s^{-1}$

Avec un courant ayant une vitesse de $1~m.s^{-1}$ le rejet est cette fois-ci presque horizontal et la dilution est réalisée de manière efficace (voir Figure 13).

Pour un courant de $5~m.s^{-1}$ on obtient le résultat suivant :

Figure 14 : Champ de concentration en saumure avec  courant transverse de $5~m.s^{-1}$

Un courant de $5~m.s^{-1}$ la dilution reste bonne (voir Figure 14). Toutefois, on voit dans cette simulation que le courant plaque le rejet au fond de l'océan. Ceci peut s'avérer être un problème car il peut alors y avoir une possibilité de stockage des saumures au fond de l'océan.

La figure suivante montre les profils de concentration en sortie du domaine (sur le côté droit) pour des courants de vitesses $0,5~m.s^{-1}$, $1~m.s^{-1}$ et $5~m.s^{-1}$ :

Figure 15 : Profils de concentration en saumure en sortie de domaine

Ainsi, on peut voir que plus le courant est fort et plus la dilution est importante, comme on pouvait s'y attendre (voir Figure 15). Un courant de $0,5~ms^{-1}$ (en rouge sur la figure) permet néanmoins de diluer presque dix fois la concentration du rejet au bout de seulement $7~m$ de transport horizontal.

Conclusion de l'étude numérique

Les résultats obtenus ici montrent que le rejet de saumure dans l'Océan Atlantique est envisageable. En outre la présence de courant marin augmente de manière significative la dilution des rejets de saumures.

L'étude numérique réalisée ici a été faite avec une géométrie simple : une seule buse de rejet a été étudiée. Pour aller plus loin, il conviendrait de réaliser une étude avec plusieurs buses de plus petits diamètres qui augmenteraient la dilution. Enfin, cette étude pourrait être complétée par une étude réalisée sur un domaine plus vaste et sur un temps plus long avec un logiciel adapté tel que TELEMAC 2D.