26 février 1997


    BES TELEMAC 2D   
par
Julie Mautin et
Steeve Champagneux


Sommaire:

0.Préambule
I. Etude d'une rupture de barrage sur fond sec
  I.1.Position du problème et Objectifs
  I.2. Moyens mis en oeuvre et expériences numériques
  I.3. Résultats et Conclusion
II. Etude d'une rupture de barrage sur fond mouillé
  II.1.Position du problème et Objectifs
  II.2. Moyens mis en oeuvre et expériences numériques
  II.3. Résultats et Conclusion
III. Etude d'un cas réel : Modélisation numérique du comportement hydrodynamique de la baie du Mont Saint Michel

Bibliographie:

Toutes les documentations relatives à l'utilisation du système de modélisation TELEMAC 2D.



retour au sommaire des BES


0.Préambule

Nous avions pour tâche, dans ce BES, d'étudier dans des situations où l'on peut déterminer des solutions analytiques simples, l'influence de certains paramètres numériques tels que la discrétisation temporelle et spatiale ou le type de schéma numérique utilisé.

Rappel : le système Telemac2d résout les équations de Barré de Saint Venant, (Navier-Stokes intégrée sur la hauteur d'eau) 2D, incompressibles, instationnaires. Il s'agit d'un système d'équations hyperboliques. Nous ne traiterons ici que des fluides parfaits.

I. Etude d'une rupture de barrage sur fond sec.

I.1.Position du problème et Objectifs.

Ce problème est analogue à celui du tube à choc étudié dans le BES SCL.
La hauteur d'eau forme à t<0 un créneau comme ce cela est décrit dans l'intitulé, puis à t=0 le barrage est retiré, on étudie alors la propagation du fluide dans le domaine sec.

Il s'agit d'un problème 1D que nous traitons ici en 2D afin d'observer une éventuelle anisotropie de la solution introduite par le schéma numérique.

Ce problème admet une solution analytique qui est donnée par la méthode des caractéristiques. Cette solution étant fournie nous ne la démontrerons pas ici.

C0**2 = g*H0

H = 4*(X/T+C0)/3

U = 2*(C0-X/2T)**2/9g

Codage de la solution analytique : hauteur d'eau et vitesse

Conditions initiales imposées : elle sont prises égales à la solution analytique en t=0.9s

Les objectifs sont ici de qualifier le comportement de la solution numérique vis-à-vis des paramètres numériques que sont le pas de temps, le pas d'espace, et le schéma numérique employé.

Remarque : nous ne présenterons pas de diagramme type H(x), U(x), H(t), U(t) car l'étude des erreurs nous a semblée plus importante. Vous pouvez toujours vous référer à un autre BE.

I.2. Moyens mis en oeuvre et expériences numériques.

Calcul des erreurs Ex et Ey

Un outil d'analyse de l'erreur commise par le code de calcul par rapport à la solution analytique est la détermination des erreurs EZx et EZy définie comme les bilans des erreurs sur la quantité Z intégrées selon les direction y et x. Ainsi, Ex permet de quantifier l'erreur selon l'axe x est Ey selon x (notations pas géniales , mais bon.....!!)

EZy =1/NBX* SQRT(SOMME SUR NBX(Zanalytique-Znumérique)**2)

REMARQUE : Il s'agit en fait de la distance entre la quantité analytique et la quantité numérique pour la norme 2 classique dont on détaille ici le calcul en FORTRAN. Il faut noter que l'ajout de ces nouvelles variables calculées entraîne une profonde modification des programmes principaux à savoir l'ajout et la modification des routines NOMVAR et PRERES ce qui a demandé un certain temps étant donné le manque de documentation claire lorsque l'on dépasse le nombre de variables utilisateur autorisées.

Moyenne temporelles

Il est clair que les ces erreurs sont instationnaires, et une fois que l'on a dégagé les phénomènes lié à leur propagation, qui sans être triviaux noyent plus les informations intéressantes qu'elles n'apportent d'éléments d'analyse, nous avons décidé de prendre la moyenne temporelle entre les temps tinitial et tfinal où tfinal est atteint lorsque la singularité qui apparaît au niveau du ressaut a remonté l'écoulement jusqu'à la limite gauche du domaine (t=1.8 secondes).

expériences numériques réalisées :

* 3 maillages différents
      50*20, 100*20, 200*20 en x,y

* 3 pas de temps différents fixés arbitrairement (Il eut certes été plus judicieux de mener une étude préliminaire afin de déterminer l'intervalle de pas de temps sur lequel le comportement de l'erreur varie de facon significative).

Les divers pas de temps envisagés sont dt=0.005, dt=0.01, dt=0.02.

* 2 schéma différents
Dans les deux cas envisagés, le schéma utilisé pour la vitesse est de type méthode des caractéristiques. On ne modifie que le schéma portant sur la hauteur d'eau.

premier cas : schéma semi-implicite centré+décentrement SUPG

second cas :  Méthode F.C.T.S.

paramètres numériques

Il sont explicités dans le fichier CAS.
À titre indicatif, le solveur utilisé est de type gradient conjugué sur équation normale avec un préconditionnement diagonal.

I.3. Résultats et Conclusion.

Nous avons donc, dans cette étude de rupture sur fond sec, étudié l'importance de différents paramètres :

Cette étude se fait en observant la moyenne temporelle de l'erreur.

attention : notation : tempehx est la moyenne temporelle de l'erreur sur la hauteur d'eau intégrée sur l'axe y.

D'autre part, la notation "caractéristiques" est à remplacer par Schéma semi-implicite centré+SUPG"

les résultats graphiques concernant tempehx et tempeux ne sont pas à prendre en compte directement car nous avons commis une erreur de normalisation. Ainsi il faut multiplier par 2.5 pour le maillage 50, par 5 pour le maillage 100 et par 10 pour le maillage 200.

D'autre part nous avons noté l'existence d'un défaut du maillage qui introduit des erreurs sur EZy. En effet, la condition limite de sortie est dissymétrique (le point 51 est solide alors que le point 71 est liquide au moins pour le maillage 50) ce qui introduit de facon naturelle une certaine anisotropie qui influe sur la solution au centre du domaine bien que les erreurs calculées ne prennent pas en compte les points frontières.

a/ Etude de l'influence du pas de temps

Nous nous sommes ici intéressés à l'influence du pas de temps pour les différents schémas et les différents maillages.

_______ dt=0.005s       _______ dt=0.01s      _______ dt=0.02s

_______ dt=0.005s       _______ dt=0.01s      _______ dt=0.02s

_______ dt=0.005s       _______ dt=0.01s      _______ dt=0.02s

_______ dt=0.005s       _______ dt=0.01s      _______ dt=0.02s

_______ dt=0.005s       _______ dt=0.01s      _______ dt=0.02s

_______ dt=0.005s       _______ dt=0.01s      _______ dt=0.02s

Dans tous les cas étudiés, nous remarquons que les erreurs sont assez faibles. Les résultats obtenus numériquement sont donc assez proches de la solution exacte. En fait ces erreurs sont surtout importantes près des limites du domaines. Les conditions aux limites choisies nt donc une forte influence sur la solution. Pour étudier les résultats sans tenir compte des effets des conditions aux limites nous devons nous éloigner des limites du domaine, il faut essayer de dégager des paliers, ce qui n'est pas toujours évident surtout pour l'erreur eux.

Nous observons que pour le maillage 200, l'erreur diminue bien lorsqu'on diminue le pas de temps et ceci pour les 2 schémas étudiés. Par contre un comportement différent est observé pour les maillages moins raffinés. En effet, pour les erreurs ehy et euy (c'est-à-dire les erreurs sur la hauteur d'eau et la vitesse intégrées sur Y) l'erreur minimale correspond au plus grand pas de temps dans le cas du maillage 50. Pour le maillage 100, il semble que notre pas de temps de 0.005s correspondent à un pas de temps critique quelque soit le schéma utilisé. Ainsi nous observons une erreur 10 fois plus importante que pour les pas de temps moins faibles.

b/ Influence du pas d'espace :

Nous avons représenté ici, les erreurs pour chacun des schémas utilisés en fonction du pas de discrétisation en X. Nous allons comparer l'influence de ce pas d'espace.

_______Maillage 50       _______Maillage 100      _______Maillage 200

_______ Maillage 50       _______ Maillage 100      _______Maillage 200

Pour les deux schémas utilisés nous observons bien une diminution de l'erreur lorsqu'on raffine le maillage. Il faut prendre en considération le fait que nous avons choisi d'effectuer cette comparaison pour un pas de temps ne posant pas de problème. Si nous avions choisi de nous placer à un pas de temps de 0.005s nous aurions retrouver les comportements mis en évidence dans le paragraphe précédent.

c/ Influence du shéma :

 _______ Schéma FCT      ________Schéma caractéristiques

Nous remarquons que quelque soit le schéma utilisé, les résultats sont à peu près équivalent. On remarque cependant que le schéma aux caractéristiques tend à minimiser les erreurs intégrées sur y, le comportement contraire étant observé pour les erreurs intégrées sur x.

estimation de l'ordre de précision des schéma sur H

Nous avons tenté d'estimer un palier de l'erreur au niveau du milieu du domaine afin de déterminer l'odre de précision du schéma en temps et en espace en cherchant la pente des courbes LOG(E)= alpha*LOG(dx,dt) à l'aide de 3 points d'où une approximation importante.

les résultats obtenus sont pour le premier schéma :

log E = alpha*log dx

dt          alpha

0.005     1.03
0.02       0.92
0.01       0.90

log E = alpha*log dt

dt          alpha

50          0.54 trés mauvaise estimation
100        1.54
200        1.45

pour le second schéma :

log E = alpha*log dx

dt          alpha

0.005     0.66
0.02       0.90
0.01       0.93

log E = alpha*log dt

dt          alpha

50          0.097 trés mauvaise estimation
100        0.84
200        1.13

Ainsi si l'on en croit nos résultats, nous pouvons conclure que le premier schéma est précis en temps à l'ordre 1 et en espace à l'odre 1,5 !!!! et que le second schéma est précis en temps à l'ordre 1 et en espace à l'odre 1.

II. Etude d'une rupture de barrage sur fond mouillé.

II.1.Position du problème et Objectifs.

Le problème étudié ici est analogue au premier cas à la différence prés que le coté droit du barrage est mouillé sur une hauteur de 1m dans notre cas. A t<0, on à donc un créneau, puisà t>0 le barrage est supprimé et on étudie la propagation de l'eau dans le domaine. Le problème est toujours hyperbolique avec en plus la présence d'un choc. Faute de temps nous n'étudierons que l'influence des schéma sur un maillage (le 200) et pour un pas de temps (dt=0.01s). Le calcul est effectué à partir de t=0

Détermination de la solution analytique

L'utilisation de la méthode des caractéristiques, de l'invariant de Riemann au niveau du ressaut et des relations de Rankine-Hugoniot au travers du choc conduisent à la résolution du système d'inconnues s, ug, yg qui s'écrit avec les notations du texte et aprés quelques transformations

c1 = SQRT(g*y1)

2*c1 = ug + 2*SQRT(g*yg)

s*(y2 - yg) = -ug*yg

s*ug*yg = -g/2*(y2**2-yg**2)+ug**2*yg

que l'on peut écrire aprés transformations

c1 = SQRT(g*y1)

s*(y2 - yg) = -ug*yg

g/2*(yg-y2)**2*(yg+y2) = y2*yg*ug**2

g*yg = c1**2+ug**2/2-c1*ug

système que l'on résout de facon numérique en déterminant l'intersection des deux dernières équations.

on obtient les valeurs numériques suivantes pour y1 = 3 m et y2 = 1 m

A l'aide de ces paramètres, on code alors la solution analytique sur fond mouillé comme suit sachant que la forme des équations est analogue à des temps prés á la solution analytique sur fond sec dans la zone comprise entre l'amont au ressaut et l'amont immédiat du choc.

Pour plus de détails, le fichier principal FORTRAN est disponible : PRINCI.F

Codage de la solution analytique : hauteur d'eau et vitesse

Conditions initiales imposées : elle sont prises égales à la solution analytique en t=0s soit un créneau avec à gauche H = 3 m et à droite H = 1 m.

II.2. Moyens mis en oeuvre et expériences numériques.

Les moyens mis en oeuvres sont les mêmes que dans le premier cas, à savoir le calcul des erreurs et la moyenne temporelle.

Les schéma étudiés sont :

* Méthode des caractéristiques pour U dans tous les cas.

* Schéma 1 : semi-implicite centré + SUPG égal à 1

* Schéma 2 : semi-implicite centré + SUPG égal au nombre de courant

* Schéma 3 : FCTS

* Schéma 4 : PSI

II.3. Résultats et Conclusion.

Nous n'avons pas pu effectuer les mêmes analyses que sur le cas fond sec surtout faute de temps.

Une évolution temporelle pour un schéma donné est représentée, ainsi que l'évolution de l'erreur sur H et U pour les quatre schéma numériques.

Le choc est correctement propagé vers l'aval, des oscillations probablement dues à l'aspect dispersif du schéma apparaissent au sommet du choc.

Les erreurs présentées doivent être corrigées d'un facteur multiplicatif 10. Les erreurs obtenues sur H sont de l'odre du pourcent alors que les erreurs sur U sont de l'odre de la dizaine de pourcent.

en noir, le schéma 1
en rouge, le schéma 2
en bleu, le schéma 3
en orange, le schéma 4

On constate que le schéma psi (4) minimise clairement l'erreur commise sur H et U.

Remarque : Ce BES nous a permis de prendre contact avec le système Telemac2d et de d'en maîtriser les principes fondamentaux. Ainsi, nous avons pu mettre en application les connaissances acquises ici sur un cas réel et complexe : la modélisation numérique du comportement hydrodynamique de la baie du Mont Saint Michel.



retour au sommaire des BES