Algorithme et Schéma numérique résolu par MATLAB

Problème :

Le programme Matlab que nous avons écrit va nous permettre de trouver la hauteur de la nappe libre au sein de la digue en fonction du temps. Pour cela et comme expliqué précédemment, nous utilisons l'équation de la conservation de la charge de DUPUIT-BOUSSINESQ au sein d'un milieu poreux:

$\displaystyle {\frac {\epsilon}{2\,K}\frac {1}{h} \frac{\partial h}{\partial t}=\frac {{\partial }^2 h}{{\partial x}^2}}$

Schéma Numérique :

On discrétise l'équation de Dupuit-Boussinesq en utilisant la méthode des différences finies. En notant $h_k^n$ la hauteur d'eau au pas de temps n et au pas d'espace k, on obtient la discrétisation suivante :

$\displaystyle {{\frac{\epsilon}{2\,K} \frac{1}{h_k^n} \frac{h_k^{n+1}\ -\ h_{k}^{n}}{\Delta t}} = \Big( {\frac{h_{k+1}^{n+1}\ -\ 2h_{k}^{n+1}\ +\ h_{k-1}^{n+1}}{2\Delta x^2} + \frac{h_{k+1}^{n}\ -\ 2h_{k}^{n}\ +\ h_{k-1}^{n}}{2\Delta x^2}}\Big) }$

 

 

En isolant le terme $h_k^{n+1}$ , nous obtenons le schéma numérique Crank-Nicholson centré, consistant en temps et en espace à l'ordre 2 :

 

$\displaystyle {\displaystyle{h_k^{n+1} =  \biggl( {\frac{h_{k+1}^{n+1}\ +\ h_{k-1}^{n+1}\ +\ h_{k+1}^{n}\ +\ h_{k-1}^{n}-\ 2h_k^{n}}{2\Delta x^2} + \frac {\epsilon}{2K\Delta t}} \biggr) \frac{1}{{1\over {\Delta x^2}}+{\epsilon \over {2K\Delta t}}{ \frac 1{h_k^n}}}\; +o \bigl( \Delta x^2 \bigr) +o \bigl( \Delta t^2 \bigr)}}$

 

Le nombre de courant associé à ce schéma numérique vaut : $\displaystyle {Co=\frac {K\Delta t}{\epsilon \Delta x^2}}$

Algorithme du calcul Matlab :

Le programme que nous avons écrit avec Matlab utilise le schéma numérique présenté ci-dessus. Il rempli ainsi une matrice avec les valeur de h aux différents pas de temps et d'espace. Le schéma est explicite, c'est à dire que pour connaître la valeur de h au pas d'espace  k il nous faut, entres autres, les valeurs aux pas d'espace k+1 et k-1 du même pas de temps.


Schéma de la matrice H regroupant toutes les valeurs de hauteurs d'eau

 

Pour résoudre ce problème, nous allons remplacer la valeur de $h_{k+1}^{n+1}$ par la valeur de $h_{k+1}^{n}$ (au pas de temps précédent). La valeur de $h_k^{n+1}$  sera donc légèrement erronée. Plus on avance dans les pas d'espace, plus l'erreur sur chaque h est potentiellement élevée.

En fin de domaine, on "rattache" le calcul de h en utilisant la condition à la limite droite. En faisant tourner plusieurs fois cette boucle en espace, nous convergeons vers la solution de h réelle. Le critère de convergence permet de régler l'erreur maximale tolérable. La boucle tourne donc tant que ce critère n'est pas respecté.

Cette action est répétée à chaque pas de temps du calcul (chaque ligne de la matrice).

Le script de ce programme est disponible dans la partie suivante : programmes MATLAB.

Effet de la stratification :

La digue que nous considérons est stratifiée, c'est pourquoi il est nécessaire de prendre à chaque pas de temps et chaque pas d'espace le paramètre $\displaystyle {A={\epsilon\over{2K}}}$ correspondant. Pour cela, la fonction Avar calcule, à partir de la hauteur d'eau $h_k^n$, la valeur de A utilisée pour calculer  $h_k^{n+1}$.

Nous considérons donc une variation $\delta A << 1 \text{ avec    } \delta A = A(h_k^{n+1}) - A(h_k^n)$.

Le calcul de cette quantité $\delta A$ dans notre programme nous a montré que les valeur étaient extrêmement faibles. Nous avons donc décidé de calculer A à tous les pas d'espace mais pas pour tous les pas de temps. Le paramètre $f_{max}$ de notre programme détermine cette période de calcul de A (3 par défaut).

Voici maintenant la façon dont est calculée A dans la fonction 'Avar.m' :

Le paramètre A, et donc $\epsilon$ et K, doit être la moyenne arithmétique des A situés entre le zéro hydrographique et la hauteur h considéré. En d'autre terme, si on a le milieu stratifié suivant :


Caractérisation du paramètre A par la stratification

 

On obtient la valeur de A pour une hauteur d'eau h appartenant à la strate numéro i par la formule suivante :

 

$\displaystyle {A\ (h \ \mathcal{2} \ i^{ème} \text{strate}) = \int\limits_{0_{{h}}}^h A(z)\, \mathrm dz = \displaystyle\sum_{k=1}^{i-1} A_k (z_k-z_{k-1}) \ \ \ +(h-z_i) A_i  }$

La fonction Avar.m cherche donc d'abord à déterminer à quelle strate h appartient puis elle applique la formule ci-dessus.

Calcul des vitesses et des débits :

Le calcul des vitesses et des débits découle immédiatement des équations présentées dans la section précédente. On obtient les schémas numériques respectifs suivants :

$\displaystyle {U_k^n = - K_k^n\ \frac{h_{k+1}^{n}\ -\ h_{k-1}^{n}\ +\ h_{k+1}^{n-1}\ -\ h_{k-1}^{n-1}}{4\Delta x}}$

$\displaystyle {Q_k^n = U_k^n\ \frac{h_k^{n-1}\ +\ h_k^{n}}{2}}$

 

 

Voici donc l'algorithme général du programme Matlab 'solveur.m' :


Algorithme général du programme solveur

 


Retour : Mise en équations

Suivant : Programmes MATLAB

Accueil Modélisation des écoulements dans la digue