\documentclass[11pt]{article}
\usepackage{graphicx, algorithm,pdfsync, bm,color}
\usepackage[utf8]{inputenc}
\usepackage[french]{babel}
\usepackage{hyperref}


\def\LINE#1{\hbox to \hsize{#1}}
\def\p{\partial}
\def\d{\hbox{d}}

\title{Projet Master 2, IMSP FreeFem++,  2018: Pollution du lac Nokoué}

\author{ Fr{\'e}d{\'e}ric Hecht}

\begin{document}
\maketitle

\tableofcontents

\section{G{\'e}n{\'e}ralit{\'e}s}
\paragraph{G{\'e}om{\'e}tries}
On s'int{\'e}resse {\`a} la mod{\'e}lisation et la simulation de l'effet d'une
rivi{\`e}re polluante sur le lac Nokoué. Les rivi{\`e}res ne sont pas dans le
domaine de calcul, seules leurs embouchures dans le lac le sont;
la g{\'e}om{\'e}trie du probl{\`e}me est alors constitu{\'e}e d'un lac connect{\'e} {\`a} 3 rivi{\`e}res $\{\Gamma_i\}_{1,2,3}$ ; on suppose par exemple que l'eau entre dans le lac par $\Gamma_1$ et $\Gamma_2$ et sort par $\Gamma_3$.
Comme nous faisons une approximation bidimensionnelle
du probl{\`e}me le lac est mod{\'e}lis{\'e} par un domaine plan $\Omega$ de
fronti{\`e}re $\Gamma=\Gamma_0\cup\Gamma_1\cup\Gamma_2\cup\Gamma_3$ o{\`u}
$\Gamma_0$ d{\'e}signe les rives du lac (voir figure \ref{figone}).
\begin{figure}
\label{figone} 
\LINE{\hss
\includegraphics[width=5cm]{Graphic1}\hss
%\hskip-8cm
\includegraphics[width=5cm]{Lac-Noukoue}\hss
%\vskip-1cm
}
\caption{G{\'e}om{\'e}tries pour les calculs.}
\end{figure}

\paragraph{\bf L'eau}
L'{\'e}coulement dans le lac est suppos{\'e} incompressible, non-visqueux et sans tourbillon (irrotationnel), ce qui permet l'existence d'une fonction potentiel d'{\'e}coulement $\phi$ telle que
\begin{equation}\label{0}
 u=\nabla\phi,~~~-\Delta\phi=0,~~ {\partial\phi\over\partial n}|_{\Gamma_i}=a_i,~i>0,~
~{\partial\phi\over\partial n}|_{\Gamma_0}=0
\end{equation}
 o{\`u} $u$
est la vitesse et $p$ la pression dans l'eau. Les constantes $a_i$ sont des donn{\'e}es d{\'e}terminant la
quantit{\'e} d'eau qui entre/sort du lac par ces fronti{\`e}res.

\paragraph{\bf Les pollutions solides}
Les pollutions solides (papiers, sacs plastiques, feuilles mortes...) sont mod{\'e}lis{\'e}es
par des particules ponctuelles de type ``scalaires passifs",
c'est {\`a} dire un ensemble de M points de $R^2$ convect{\'e}s par l'eau.  La position $X_j(t)$ au temps $t$  d'une particule $j$ qui {\'e}tait en $X_j^0$ {\`a} l'instant initial est donc
donn{\'e}e par
\begin{equation}\label{1}
    \dot{X_j}(t) = u(X_j(t))\quad X_j(0) = X_j^0,~~~ j=1..M
\end{equation}
On peut approcher ces {\'e}quations par
\begin{equation}\label{1001}
    X_j^{n+1} = X_j^n + u(X_j^n)\delta t,\quad X_j^0 \hbox{ donn{\'e}}.
\end{equation}

 Si
toutes les particules n'ont pas le m{\^e}me effet polluant (taille du
sac plastique par exemple) il convient d'associer un poids $w_j$ {\`a}
la j-{\`e}me particule; celui-ci reste constant dans le temps sauf si la pollution est bio-d{\'e}gradable.

\paragraph{\bf Densit{\'e} de pollution}

Soit  $\{X_j(t),w_j\}_1^M$ un ensemble de $M$ particules
ponctuelles de poids $\{w_j(t)\}_0^M$ convect{\'e}es par l'eau mais aussi
agit{\'e}s par un mouvement brownien (prise en compte de la turbulence dans l'eau du lac),
alors
\begin{equation}\label{eds}
     X_j^{n+1} = X_j^n + u(X_j^n)\delta t + \sqrt{2\sigma\delta t}{\cal N}^n_j,
     \quad X_j^0 \hbox{ donn{\'e}}.
\end{equation}
 o{\`u} ${\cal N}^n_j$ est une r{\'e}alisation de ${\cal N}$ qui est une variable al{\'e}atoire gaussienne centr{\'e}e en 0, de variance unit{\'e}
 et {\`a} valeurs dans $R^2$.
On suppose donn{\'e}e la loi de variation des poids $w$:
\[
    \dot{w}_j = a(X_j(t)) w_j(t) + b(X_j(t))
\]
On pose
\begin{equation}\label{2}
    \rho(x,t)=\lim_{M\to\infty,~\epsilon\to 0}
    \sum_{\{j~:~X_j(t)\in B(x)\}} w_j(t)
\end{equation}
o{\`u} $B(x)$ est un disque de centre $x$ et de rayon $\epsilon$
petit. Dans la limite $M=O(\epsilon^{-2})$, $\epsilon\to 0$ et avec
quelques hypoth{\`e}ses suppl{\'e}mentaires on peut
montrer que
\[
     \partial_t\rho+u\cdot\nabla\rho  - \sigma\Delta\rho =a\rho+b,
     \qquad \rho(x,0)=\lim_{\epsilon\to 0}\sum_{j:X^0_j\in B(x)} w_j(0)
\]
Pour que le probl{\`e}me soit bien pos{\'e} il faut se donner $\rho$ ou $\p_n\rho$ sur les bords.
Si on suppose que les polluants rentrent par $\Gamma_1\cup\Gamma_2$, sortent par $\Gamma_3$ et glissent sur les
bords du lac alors  on prendra
\[
\rho=\rho_\Gamma\hbox{ sur } \Gamma_i~,i=0,1,2 ~~~ \p_n\rho=0\hbox{ sur }\Gamma_3
\]
avec $\rho_\Gamma$ donn{\'e} mais nul sur $\Gamma_0$.

\paragraph{\bf Polluants liquides}
La pollution liquide consiste en g{\'e}n{\'e}ral  soit de produits chimiques dissous dans de l'eau (ex: sels), soit de polluants liquides non miscibles {\`a} l'eau (ex: mazout).  Chacun est mod{\'e}lis{\'e} par une concentration $\{\rho_i\}_1^I$ dans l'eau du lac et en l'absence de r{\'e}action chimique on peut calculer les $\rho_i$ par
\begin{equation}
     \partial_t\rho_i+u\cdot\nabla\rho_i  - \sigma_i\Delta\rho_i =0,
     \qquad \rho(x,0)=\rho^0_i|_\Gamma,\quad\rho_i|_\Gamma\hbox{ ou }\p_n\rho_i\mbox{ donn{\'e}s}\label{eq 6}
\end{equation}

o{\`u} $\sigma_i$ est le coefficient de diffusion de l'esp{\`e}ce $i$.

\paragraph{\bf La chimie}
Prenons le cas de la d{\'e}composition d'un d{\'e}tergent chlor{\'e}, le
chlorure d'hydrog{\`e}ne qui se combine avec un hydrochlorure de
sodium pour donner du sel et de l'eau.
\begin{equation}\label{chem}
                \texttt{HCl}+\texttt{NaOH} \to \texttt{NaCl} + \texttt{H}_2\texttt{O}
\end{equation}
Nous avons {\`a} g{\'e}rer 3 esp{\`e}ces dissoutes dans l'eau, chacune ayant
sa concentration, $\rho_H,\rho_N,\rho_S$ pour
\texttt{HCl,NaOH,NaCl}. On suppose que le taux de r{\'e}action est
proportionnel aux densit{\'e}s des composants ; par ailleurs la
concentration en sel ne nous int{\'e}resse pas. On aura donc pour les
deux autres esp{\`e}ces
\begin{eqnarray}\label{model}&&
    \partial_t\rho_H+u\cdot\nabla\rho_H - \sigma_H\Delta\rho_H
     =-\alpha\rho_H\rho_N
     \cr&&
    \partial_t\rho_N+u\cdot\nabla\rho_N -
    \sigma_N\Delta\rho_N
     =-\alpha\rho_N\rho_H
\end{eqnarray}
En effet la r{\'e}action (\ref{chem}) nous dit qu'il faut une quantit{\'e}
{\'e}gale des deux esp{\`e}ces pour les faire dispara{\^\i}tre. S'il n'y avait
pas de courant ni de diffusion on aurait simplement
\[
    \dot\rho_H(t)=-\alpha\rho_H(t)\rho_N(t),\qquad
    \dot\rho_N(t)=-\alpha\rho_N(t)\rho_H(t)
\]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Projet 1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%On s'int{\'e}resse {\`a} l'{\'e}volution d'un polluant de densit{\'e} $\rho(x,t)$ dans une rivi{\`e}re assimil{\'e}e {\`a} un segment de droite $\Omega=(0,L)$.  La vitesse de l'{\'e}coulement est constante en $x$ et donn{\'e}e $u(t)$ en fonction du temps,  la diffusion du polluant $\sigma$ est une constante donn{\'e}e; $\rho(\cdot,0)$ est donn{\'e}  sur $(0,L)$ ainsi que $\rho(0,t)$ et  $\rho(L,t)$ pour tout temps $t\in(0,T)$.
%
%\begin{description}
%\item[Q1] Calculer par diff{\'e}rences finies en temps et en espace et avec un sch{\'e}ma d'Euler implicite la solution de
%\[
%\p_t\rho-\sigma\p_{xx}\rho=0\hbox{ dans }(0,L)\times(0,T),~~\rho(0,t)=\rho(L,t)=0,~\forall t\in(0,T)
%\]
%lorsque $\rho(x,0)=\sin(\pi\frac{x}L), L=1,\sigma=0.1, T=1$.
%On calculera l'erreur $e$ avec la solution analytique et $\rho_e$ et on fera avec \texttt{gnuplot} un graphique de $\log(e)$ versus $\log N$
%o{\`u} $N$ est le nombre de points de discr{\'e}tisation de $(0,L)$.  La programmation doit {\^e}tre faite en C++.
%
%\item[Q2]
%Rajouter dans l'EDP le terme $u\p_x\rho$, discr{\'e}tis{\'e} par un sch{\'e}ma aux diff{\'e}rences finies d{\'e}centr{\'e} selon le signe de $u$.  Faire les simulations avec les m{\^e}mes param{\`e}tres sauf $\rho(x,0)=\exp(100(x-0.5)^2)$, et $u=\sin(10 t)$.
%
%\item[Q3] Un polluant solide, simul{\'e} par une seule particule en position $X^0=L/2$ {\`a} l'instant $t=0$,  attaqu{\'e} par le polluant liquide de concentration $\rho$ a pour {\'e}quations discr{\`e}tes:
%\begin{eqnarray}&&
%     X^{n+1} = X^n + u(t^n)\delta t + \sqrt{2\sigma\delta t}{\cal N}^n,
%\cr&&
% w^{n+1}=w^n- w^n\rho(X^n,n\delta t)\delta t,~~w^0=1.
%\end{eqnarray}
%o{\`u} ${\cal N}^n$ est une r{\'e}alisation de ${\cal N}$ g{\'e}n{\'e}r{\'e}e par la formule de Cox-Muller:
%\[
%{\cal N}^n=\sqrt{-2\log Z^n_1}\cos(2\pi Z^n_2)
%\]
%avec $Z_1,Z_2$ 2 nombres al{\'e}atoires g{\'e}n{\'e}r{\'e}s par la \texttt{(1.+rand())/(1.+RAND\_MAX)}.
%Tracer $t\to w(t)$.
% \end{description}
%
% \pagebreak

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Le projet}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

On s'int{\'e}resse maintenant {\`a} l'{\'e}volution d'un polluant de densit{\'e} $\rho_1(x,t)$ dans le lac Noukoué repr{\'e}sent{\'e} par le domaine $\Omega$. 

\subsection{Calcul de l'{\'e}coulement fluide}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nous rappelons que
l'{\'e}coulement dans le lac est suppos{\'e} incompressible, non-visqueux et sans tourbillon (irrotationnel). Il existe donc une fonction potentiel d'{\'e}coulement $\phi$ telle que
\begin{equation}\label{eq-u}
 {\bf u}=\nabla\phi
\end{equation}
o{\`u} ${\bf u} =(u_x, u_y)$ est la vitesse. La conservation de la masse s'{\'e}crivant pour le fluide incompressible
\begin{equation}
\label{eq-masse}
  \nabla. {\bf u} =0,
\end{equation}
nous obtenons de (\ref{eq-u}) l'EDP suivante pour le potentiel $\phi$:
\begin{equation}\label{eq-phi}
 -\Delta\phi=0\quad  \mbox{dans } \Omega, 
\end{equation}
avec des conditions aux limites qui seront pr{\'e}cis{\'e}es plus bas.

L'{\'e}coulement fluide peut-{\^e}tre d{\'e}crit {\'e}galement en introduisant la fonction de courant $\psi$ par
\begin{equation}
  \label{eq-upsi}
  u_x = \frac{\partial \psi}{\partial y}, \quad u_y = - \frac{\partial \psi}{\partial x}.
\end{equation}
La conservation de la masse (\ref{eq-masse}) est ainsi satisfaite. L'{\'e}coulement {\'e}tant consid{\'e}r{\'e} sans rotationnel (sans tourbillon), on peut {\'e}crire:
\begin{equation}
  0 = \nabla \times {\bf u} = \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y},
\end{equation}
et finalement, en utilisant (\ref{eq-upsi}), l'{\'e}quation pour $\psi$:
\begin{equation}
  \label{eq-psi}
  -\Delta\psi=0\quad  \mbox{dans } \Omega.
\end{equation}



Consid{\'e}rons pour commencer le domaine de calcul $\Omega$ repr{\'e}sent{\'e} sur la figure \ref{fig-rect}. La fronti{\`e}re du domaine est
\begin{equation}
  \Gamma = \Gamma_0 \cup \Gamma_1 \cup \Gamma_2 \cup \Gamma_3 ,
\end{equation}
o{\`u} $\Gamma_1$ et $\Gamma_2$ mod{\'e}lisent la section transversale de deux rivi{\`e}res se jetant dans le lac (fluide entrant) et $\Gamma_3$ est la fronti{\`e}re par laquelle le fluide sort du domaine (fluide sortant). Les autres fronti{\`e}res repr{\'e}sentent les berges. 

\begin{figure}
\centerline{\includegraphics[width=0.8\textwidth]{sketch_dom_rect}}
\caption{Domaine de calcul $\Omega$.}
\label{fig-rect}
\end{figure}

 
 



 \begin{description}

\item[{\bf Q1}] Construire avec FreeFem++ le maillage repr{\'e}sent{\'e} sur la figure \ref{fig-rect}. On gardera comme param{\`e}tres les longueurs des segments
    d{\'e}finissant les fronti{\`e}res $\Gamma_1, \Gamma_2, \Gamma_3$ qu'on va noter par  ${\ell}_1, {\ell}_2$ (largeurs des rivi{\`e}res) et ${\ell}_3$ (largeur du canal de sortie). Prendre comme valeurs par d{\'e}faut:
    \[ {\ell}_1 = {\ell}_2 = 5, \quad {\ell}_3=10 \]

\item[{\bf Q2}]  Les conditions aux limites du probl{\`e}me seront impos{\'e}es par rapport au champ de vitesses : la composante normale est nulle sur les berges 
    \begin{equation}
      {\bf u}. {\bf n} = 0 \mbox{  sur  } \Gamma_{10} \cup \Gamma_{11} \cup \Gamma_{12},
    \end{equation}
    et le fluide entrant est mod{\'e}lis{\'e} par une vitesse normale constante sur chaque section
    \begin{equation}
      {\bf u}. {\bf n} = -U_1 \mbox{  sur  } \Gamma_1, \quad {\bf u}. {\bf n} = -U_2 \mbox{  sur  } \Gamma_2,
    \end{equation}
o{\`u} ${\bf n}$ est la normale vers l'ext{\'e}rieur du domaine. 
\begin{enumerate}
  \item \`A partir de la conservation de la masse (\ref{eq-masse}) d{\'e}duire la relation de compatibilit{\'e} qui permet de calculer la vitesse normale {\`a} travers la section de sortie $\Gamma_3$.
  \item  D{\'e}duire les conditions aux limites pour l'{\'e}quation (\ref{eq-phi}) du potentiel $\phi$.      
      \item En utilisant ces conditions aux limites, {\'e}crire la formulation variationnelle de (\ref{eq-phi}).
\end{enumerate}

\item[{\bf Q3}] \'Ecrire un programme FreeFem++ pour r{\'e}soudre l'{\'e}quation (\ref{eq-phi}) pour $\phi$ en utilisant la m{\'e}thode des {\'e}l{\'e}ments finis $P^1$. Utiliser comme valeurs:
    \[U_1 = 2, \quad U_2 =1.\]
    Visualiser les lignes de niveau de $\phi$.\\
    Calculer les deux composantes de la vitesse et visualiser le champ de vecteurs (utiliser la commande \texttt{plot(phi, [ux,uy], wait=1)}).\\
    Commenter les r{\'e}sultats obtenus pour 
     \[U_1 = 1, \quad U_2 =10.\]
     
    
    

\item[{\bf Q4}] \'Ecrire les conditions aux limites et la formulation variationnelle pour la fonction de courant $\psi$. (Indication: utiliser la d{\'e}finition (\ref{eq-upsi}) et les conditions sur les vitesses -- en d{\'e}duire que $\psi=const$ sur les fronti{\`e}res mod{\'e}lisant les berges.) 

\item[{\bf Q5}] \'Ecrire un programme FreeFem++ pour r{\'e}soudre l'{\'e}quation (\ref{eq-psi}) pour $\psi$. Pour simplifier le probl{\`e}me, on consid{\`e}re:
    \begin{equation}
      \left.{\partial\psi\over\partial n}\right|_{\Gamma_1 \cup \Gamma_2 \cup \Gamma_3}=0.
    \end{equation}
    Visualiser les lignes de niveau de $\psi$. Commenter.
    

\item[{\bf Q6}] R{\'e}soudre les m{\^e}me probl{\`e}mes (calcul de $\phi$ et de $\psi$) sur le maillage du  lac Noukoué (voir figure \ref{fig msh noukoue}). Ce maillage  est stock{\'e} dans le fichier \texttt{lac-noukoue.msh} qui sera lu par le script FreeFem++.
 
 `{\color{blue}       
La fronti{\`e}re $\Gamma_1$ repr{\'e}sente l'embouchure du fleuve  Ouémé , $\Gamma_2$ l'embouchure de la lagune de Porto-Novo, $\Gamma_3$ le canal vers la mer (on négligera les effets dus aux marées), el les autres
frontières sont  defini comme suit: 
\begin{description}
\item[$\Gamma_4$ ] est comprise entre la lagune et le canal vers la mer
\item[$\Gamma_5$ ] est comprise entre le canal vers la mer et et l'embouchure du fleuve dans le lac. 
\item[$\Gamma_6$ ] est comprise entre  l'embouchure du fleuve dans le lac. et la lagune de Proto-Novo.
\end{description}
}

Le d{\'e}bit  du fleuve Ouémé est de  $170 m^3/s$ pendant le saison des pluies (moment de l'étude), le fleuve  Ouémé  se  sépare en deux branches que l'on supposera de même débit qui
se jetent respectivement  dans le lac  et  dans la lagune de Porto-Novo. Nous  n{\'e}gligerons tous les autres l'apport en eau pour fermer le modèle et  nous  supposerons  que la profondeur du lac constante et égale à $1m$.

\begin{figure}[htb]
\begin{center}
\centerline{\includegraphics[width=0.9\textwidth]{Th-noukoue}}
\caption{Maillage du lac Noukoué. \label{fig msh noukoue}}
\end{center}

\end{figure}
 \end{description}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{\'Evolution du polluant}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nous allons supposer que le polluant ${\rho_1}$ est de type liquide et miscible et que son évolution est d'écrite par l'equation (\ref{eq 6}).
Le but de cette partie est d’étudier l’influence d’une pollution provenant de la lagune de Porto-Novo et de tracer cette courbe de pollution dans l’année suivant le début de celle ci.
On négligera les apports en eau des affluent, on suppose que le lac à une profondeur constante  de 1 m.
Donc, les conditions aux limites du polluant $\rho_1$ sont sur 
sur $\Gamma_1$ et sur $\Gamma_2$ où les fonctions $g_1$, $g_2$ vont  modéliser respectivement le flux de polluant provenant  du fleuve  Ouémé et la lagune de Porto-Novo entrant dans le lac par jour ..
{\color{blue} C'est à dire que les conditions aux limites du polluant seront respectivement sur  $\Gamma_4$ ,$\Gamma_5$ ,$\Gamma_6$ :  $ \partial {\rho_1}/ \partial \bf{n} =0$, sur $\Gamma_1$ : ${\rho_1}=g_1$, sur $\Gamma_2$: $\rho_1=g_2$  et sur $\Gamma_3$ :  $ \partial {\rho_1}/ \partial \bf{n} =0$,(sortie  libre).

De plus le polluant rentre avec le debit de l'eau ( nous allons supposer que la pollution est  équi répatie dans le fleuve et la lagune) donc $g_1$ égale $g_2$ et sont des constantes en espace. La quantité de polluant entrant par unite de temps est donné par 
$ Q = \int_{\Gamma_1} g_1 \,\partial \phi/\partial \bm{n} + \int_{\Gamma_1} g_2 \, \partial \phi/\partial \bm{n}$ ce qui permet de calculer $g=g_1=g_2$ connaissant $Q$.
}

 \begin{description}
\item[{\bf Q7}] Ecrire le schéma aux différences finies implicite en temps pour la résolution de l’équation (\ref{eq 6}), puis écrire la formulation variationnelle en espace. Ecrire le schéma final correspondant la discrétisation en espace-temps (différences finies en temps et éléments finis en espace).
\item[{\bf Q8}] Calculer les valeurs des constantes pour imposer les conditions aux limites pour le potentiel de vitesses $\phi$ et pour le polluant $\rho_1$. On supposera que l’épisode de pollution dure 3 jours et a un débit de $Q=500 tonnes/jour$.

\item[{\bf Q9}]Pour un coefficient de diffusion bidimensionnel $\sigma_1 = 0.1km^2/jours$, calculer avec Freefem++, la courbe du polluant sortant à la mer  et entrant de la lagune de Porto-Novo. c’est-`a-dire tracé les deux courbes
{\color{blue} $\int_{\Gamma_2} \rho_1 \,\partial \phi/\partial \bm{n} $ et $\int_{\Gamma_3} \rho_1 \,\partial \phi/\partial \bm{n} $.}

\item[{\bf Q10}] (Question subsidiaire) Comment faire pour prendre en compte la variation du debit du fleuve au cours des 4 saisons ainsi que l'évaporation.(voir le document 
\href{http://horizon.documentation.ird.fr/exl-doc/pleins_textes/pleins_textes_4/hydrologie/16568.pdf}{Hydrologie de L'Ouémé} ).
 \end{description}
\end{document}
