SIMULATION DE LA RUPTURE D'UN BARRAGE.

I- PRESENTATION DU PROBLEME.

Nous étudions la rupture d'un barrage:

Le barrage est modélisé dans un domaine de 20m en abscisse et 5m en ordonée, le fond est plat, et le barrage est situé à l'abscisse x=10m. A l'instant initial, on ouvre les vannes. Les résultats obtenus sont de la forme ci-dessous:

Notre programme calcule les différents paramètres de l'évolution de cette rupture tel que la hauteur d'eau, la vitesse... Nous avons de plus introduit la solution analytique de cette évolution afin de la comparer à celle calculée par le code.

Nous allons tout d'abord chercher le schéma le plus précis parmi un panel de quatre schémas, puis, en se basant sur ce schéma, nous allons étudier l'influence du pas de temps, du maillage et du solveur afin d'avoir un calcul optimal avec ce schéma.

Sommaire.









II- PROGRAMMATION.

A l'instant initial, il y a une hauteur d'eau de 6m a l'amont du barrage et 2m a l'aval du barrage. Cette condition initiale est implémentée dans la subroutine CONDIN du programme PRINCI.f.

Sommaire.







Afin de pouvoir vérifier la validité des calculs effectués par le code, nous avons introduit deux nouvelles variables calculant la hauteur d'eau analytique et la vitesse analytique.

La solution analytique de la hauteur et de la vitesse est de la forme suivante:

Les deux nouvelles variables sont calculées à l'aide des subroutines NOMVAR et PRERES du programme PRINCI.f.

Dans la subroutine NOMVAR, on définit les noms et unités des variables: la variable 23 (N) contient la hauteur d'eau analytique , la variable 24 (O) contient la vitesse analytique, et les variables 25 (R) et 26 (z) contiennent les erreurs sur la hauteur et la vitesse.

Afin d'utiliser ces variables dans Rubens, il faut les appeler dans le fichier cas***.

Dans la subroutine PRERES, on implemente les calculs comme ci-dessous:

On peut alors par la suite visualiser les deux solutions pour voir leur différences comme ci-dessous:

Sommaire.








III- ANALYSE DU COMPORTEMENT DES SCHEMAS.

Afin de juger la précision des différents schémas, nous allons regarder à la fois leur influence sur la conservation de la masse ainsi que leur exactitude comparée à la solution analytique.

Nous allons étudier pour la vitesse et la hauteur d'eau deux sortes de schémas et leur couplage (les résultats ci-dessous sont effectués pour un pas de temps de 0.01 et un maillage de 200 points):

BILAN DE MASSE SELON LE SCHEMA

SCHEMA 12 12

SCHEMA 15

SCHEMA 32 12

SCHEMA 35

Temps de calcul selon le schéma:

V\H

2

5

1

7'46''

8'30''

3

7'45''

8'23''

Bilan de masse selon le schéma:

V\H

2 5

1

1.64

-0.039

3

0.487

-0.0478











VISUALISATION DES RESULTATS

SCHEMA 12 12

SCHEMA 15

SCHEMA 32 12

SCHEMA 35

Les schémas étudiés présentent peu de différences concernant la visualisation des paramètres ou des erreurs sur la hauteur et la vitesse, toutefois, ils presentent des différences au niveau des temps de calcul et des bilans de masse.

Parmi les cas présentés, les schémas les plus précis nous semble être les schémas 35 et 15, le 15 nous paraissant le plus précis concernant l'abcisse de la vague aval pour le temps représenté qui nous semble être un critère essentiel lorsque l'on simule la rupture d'un barrage. Le 35 par contre nous paraissait engendrer le moins d'erreur parmi les trois schémas restants.



Afin de mieux voir la précision de ces deux schémas, nous les avons testés avec un maillage de moindre qualité qui devrait moins gommer les erreurs dues au schéma. Nous avons donc lancé le programme avec les mêmes caractéristiques mais un maillage différent.





BILAN DE MASSE SELON LE SCHEMA

SCHEMA 15

SCHEMA 35

Temps de calcul selon le schéma:

15

35

1'31''

1'44''

Bilan de masse selon le schéma:

15 35

0.0049

0.0005












VISUALISATION DES RESULTATS

SCHEMA 15

SCHEMA 35

Cette fois, les différences sont plus flagrantes et nous choississons le schéma 35 pour sa précision.

Sommaire.






IV- SENSIBILITE AUX PAS DE DISCRETISATION ET OPTIMISATION.

Nous conservons cette fois le schéma 35, et nous allons voir l'influence des autres paramêtres afin d'optimiser le calcul avec ce schéma.

Nous avons d'ores et déjà testé les maillages de 50 et 200, nous avons donc relancé un calcul avec un maillage de 100 points en largueur et présentons les résultats ci-dessous:

Temps de calcul selon le maillage:

50

100

200

1'44''

4'00''

8'23''

Bilan de masse selon le maillage:

50

100

200

0.0005

-0.0533

-0.0478





VISUALISATION DES RESULTATS

MAILLAGE 50 POINTS

MAILLAGE 100 POINTS

MAILLAGE 200 POINTS

Le maillage donnant les résultats les plus justes graphiquement parlant est le maillage de 200 points, le temps de calcul supplémentaire nous semble justifié vu les meilleurs résultats que nous obtenons. Nous allons garder ce maillage dans la suite des calculs.

On peut penser que ce schéma n'est pas trop dissipatif car l' on peut observer des wiggles au niveau des fortes variations. Le maillage 200 a masqué cet élément par sa finesse.

Sommaire.






Cette fois, en gardant un maillage de 200 points, nous allons tester quatre pas de temps différents qui sont 0.001, 0.005, 0.01 et 0.02.

Temps de calcul selon le pas de temps:

0.001

0.005

0.01

0.02

92'47''

15'49''

8'23''

***

Bilan de masse selon le pas de temps:

0.001

0.005

0.01

0.02

0.157

-0.0032

-0.0478

***

Pour des pas de temps supérieurs à 0,01 s les calculs n'ont pas abouti. Le schéma 3 est une hybridation entre une méthode des caractéristiques et de SUPG. La méthode des caractéristiques impose peut-etre une contrainte de type CFL et donc sur le pas de temps.






VISUALISATION DES RESULTATS

PAS DE TEMPS  0.001

PAS DE TEMPS  0.005

PAS DE TEMPS  0.01

PAS DE TEMPS  0.02

On remarque que plus on réduit le pas de temps, plus les erreurs sur la hauteur et la vitesse sont faibles (les courbes approche mieux les angles x=2.5 et x=7.5), toutefois les instabilités autour de x=17 (vague aval) sont plus fortes. Le meilleur compromis entre ces deux critères et donc le pas de temps 0.01. Le calcul n'a pas convergé pour des pas de temps supérieurs à 0.01 comme le montre la dernière figure, le calcul obtenu a été effectué pour le schéma 15.

Sommaire.







Il est possible de choisir le solveur utilisé pour la résolution de l'étape de propagation.

Dans tous les calculs précédents, nous avons utilisés le solveur 3 qui est la méthode du gradient conjugé sur équation normale. Nous allons par la suite tester les solveurs suivants:

Les résultats sont les suivants:

Temps de calcul selon le solveur:

1

2

3

4

***

6'54''

8'23''

8'39''

Bilan de masse selon le solveur:

1

2

3

4

***

-0.04149

-0.0478

-0.5898







VISUALISATION DES RESULTATS

SOLVEUR  1

Le calcul pour ce solveur a planté à la deuxième itération.

SOLVEUR  2

SOLVEUR  3

SOLVEUR  4

Le solveur 2 nous semble être le plus performant pour ce calcul vu sa rapidité d'execution, les trois solveurs ayant convergés donnant des résultats assez similaires.

Sommaire.










V- CONCLUSION.

Dans le cas de notre étude, nous avons essayé d'optimiser les résultats en choissisant des options pas à pas parmi un panel non exhaustif tel que le maillage, le pas de temps, le solveur et le schéma de résolution.

Dans le cas de l'étude de la rupture d'un barrage sur fond mouillé, nous avons un calcul optimal en utilisant un schéma de pondération de la méthode des caractéristiques et de la méthode SUPG pour la vitesse et le schéma UPG appliqué à l'équation de continuité pour la hauteur associé à un pas de temps de 0.01, un maillage de 200 points et en utilisant le solveur 2.

Sommaire.