ANALYSE DES EQUATIONS DE NAVIER-STOKES EN BASSIN PEU PROFOND ET DE L'EQUATION DE TRANSPORT THESE présentée à la faculté des sciences, pour obtenir le grade de docteur es sciences, par Pascal AZERAD UNIVERSITE DE NEUCHATEL Institut de Mathématiques Emile Argand 11 2007 NEUCHATEL (Suisse) ANALYSE DES EQUATIONS DE NAVIER-STOKES EN BASSIN PEU PROFOND ET DE L'EQUATION DE TRANSPORT THESE présentée à la faculté des sciences, pour obtenir le grade de docteur es sciences, par Pascal AZERAD UNIVERSITE DE NEUCHATEL Institut de Mathématiques Emile Argand 11 2007 NEUCHATEL (Suisse) Remerciements. Je tiens à remercier avant tout Olivier Besson, pour son magistère ma- thématique, son soutien indéfectible et sa disponibilité bienveillante tout au long de mes recherches. Je remercie Pierre Lesaint pour m'avoir initié au calcul scientifique et m'avoir donné le goût de l'analyse numérique. Je suis également reconnaissant à messieurs les professeurs A. Robert, J. Simon et A. Valli d'avoir bien voulu faire partie du jury, de leur lec- ture attentive du manuscrit et de leurs précieuses suggestions. J'ai tiré grand profit de discussions avec R. Stenberg. Ma collabora- tion avec Pierre Perrochet et Jérôme Pousin a été (et sera j'espère) très fructueuse. L'ambiance conviviale qui règne à l'Institut de Mathématiques de l'Uni- versité de Neuchâtel est pour beaucoup dans la réalisation de ce travail, et je suis heureux de témoigner des bonnes conditions dont j'ai bénéficié, en particulier en tant qu'assistant du cours de théorie des probabilités des professeurs Alain Valette et Monique Graf. Ma reconnaissance va également au Fonds National de la Recherche Scientifique de la Confédération Helvétique (requêtes 20-31162.91 et 20-37651.93) pour le financement de mes travaux. Ma gratitude va enfin à ma famille qui a su m'entourer de son affection et me stimuler parfois, à mes collègues qui m'ont soutenu par leur ami- tié, leur humour dévastateur (et prodigué des conseils en WTffî) durant ces dernières années à Neuchâtel: Cédric Béguin, Pierre-Alain Cherix et Claude Fuhrer. Je dédie cette thèse au lac de Neuchâtel et à ma sirène. 3 IMPRIMATUR POUR LA THÈSE Analyse des équations de Navier-Stokes en bassin peu profond et de l'équation de transport de M. Pascal Azérad UNIVERSITÉ DE NEUCHÂTEL FACULTÉ DES SCIENCES La Faculté des sciences de l'Université de Neuchâtel sur le rapport des membres du jury, Messieurs O. Besson, A. Robert, P. Lesaint (Besançon), J. Simon (Clermont-Ferrand) et A. Valli (Trento) autorise l'impression de la présente thèse. Neuchâtel, le 8 janvier 1996 Le doyen: R. Dändliker V, Table des matières 1 Motivations. 9 1.1 Mathématiques lacustres................... 9 1.2 Conventions typographiques................ 9 1.3 Notations mathématiques.................. 10 1.4 Les fluides géophysiques................... 11 1.5 Plan de la thèse....................... 12 1 L'approximation hydrostatique. 13 2 La modélisation asymptotique. 15 2.1 La faible profondeur..................... 15 2.2 Équations gouvernant l'écoulement............. 16 2.3 Changement d'échelle et renormalisation.......... 17 2.4 L'hypothèse hydrostatique.................. 18 3 Le problème de Stokes hydrostatique. 21 3.1 Préliminaire d'analyse fonctionnelle............. 21 3.2 Analyse du problème continu................ 24 3.2.1 Equations du modèle................. 24 3.2.2 Formulation variationnelle du problème...... 25 3.2.3 Existence, unicité, stabilité de la solution..... 28 3.3 Résolution du problème approché............. 30 3.3.1 Critères pour une discrétisation bien posée..... 30 3.3.2 Elément fini hydrostatique.............. 33 3.3.3 Vérification des inégalités de stabilité........ 34 3.3.4 Estimation d'erreur................. 37 5 6 TABLE DES MATIÈRES 4 Écoulement en eau peu profonde. 39 4.1 Modèle Renormalisé..................... 39 4.2 Modèle "Hydrostatique"................... 40 4.3 Espaces fonctionnels et notations utilisés.......... 40 4.4 Forme faible renormalisée.................. 42 4.5 Forme faible hydrostatique................. 43 4.6 Majorations à priori..................... 44 4.6.1 Forme faible spatiale................. 44 4.6.2 Égalité d'énergie................... 44 4.7 Existence de la solution ( cas linéarisé)........... 46 4.7.1 Passage à la limite.................. 47 4.8 Obtention de la pression................... 48 4.9 Sens de la condition initiale et unicité........... 50 4.9.1 Unicité de la solution................ 53 II Le transport. 57 5 Moindres carrés dans l'espace-temps. 59 5.1 Une idée simple sur un exemple simple........... 59 5.2 Principe des méthodes aux moindres carrés........ 60 5.3 Le cadre spatio-temporel................... 62 5.4 La méthode STILS...................... 63 5.4.1 Forme variationnelle aux moindres carrés..... 63 5.4.2 Advection = diffusion anisotrope.......... 65 5.4.3 Retour à l'exemple.................. 67 5.5 Comportement numérique.................. 68 5.5.1 La discrétisation de Galerkin............ 68 5.5.2 L'analyse de Fourier-von Neumann et la matrice d'amplification.................... 69 5.6 Exemples de simulations................... 71 6 Inégalité de Poincaré courbe 75 6.1 Cadre fonctionnel....................... 75 6.1.1 Espaces de Sobolev obliques............. 75 6.1.2 Écoulement remplissant............... 76 TABLE DES MATIÈRES 7 6.1.3 Inégalité de Poincaré courbe............. 78 6.2 Équivalence entre les formulations............. 82 A 87 A.l Élément fini hydrostatique.................. 87 A.2 Le code STILS........................ 91 7 Bibliographie 8 TABLE DES MATIÈRES Chapitre 1 Motivations. 1.1 Mathématiques lacustres. Commençons par le début. Il y a bien longtemps, les hommes choisirent d'établir leurs cités au bord de l'eau, des fleuves, des mers ou des lacs. Pour des raisons pratiques bien sûr, mais ausi peut-être parce que le chant des nymphes, des sirènes et des ondines les attirait à la claire fontaine ... Lequel d'entre nous enfant n'est pas resté interdit, absorbé dans la contemplation de l'onde courant sous le pont? Et plus tard, son coeur n'a-t-il pas vibré à l'écho du vers de Baudelaire: "Homme libre, toujours tu chériras la mer"? Comme le lac de Neuchâtel est à l'origine de cette thèse, notre humble travail s'inscrit quelque part dans une provisoire sous-section d'un chapitre d'une partie d'un tome de cette histoire, volume rangé dans une des cellules hexagonales de Celle précisément où l'immense bibliothèque de Babel [12]. tu te trouves, cher lecteur. 1.2 Conventions typographiques Pour avoir nous-mêmes maintes fois sué sang et eau à la lecture d'arides et absconses mathématiques, nous avons tenté d'aérer la mise en page, par une ample marge, permettant ainsi à notre cher lecteur de croquer quelque chimère, ou de vérifier quelque assertion au gré de son humeur. Il est dans la nature des mathématiciens d'exposer de façon deductive, logiquement ordonnée, ce qu'ils ont trouvé de façon inductive, à tâtons et dans le plus grand désordre. Ils agissent ainsi par souci de concision 9 10 CHAPITREl. MOTIVATIONS. Nous avons délibérément transgressé cette règle tacite, et de temps en temps en marge du texte nous livrons quelque indication heuristique. et d'élégance, peut-être aussi par goût du secret, lointaine survivance de l'époque pythagoricienne?... A part cette marge élargie, la numé- rotation suit la convention habituelle, hormis quelques variantes que le lecteur attentif détectera. 1.3 Notations mathématiques. Nous explicitons au fur et à mesure des besoins les notations que nous utilisons. Mentionnons seulement que C désigne une constante flottante au gré des courants ( et des équations). 1.4. LES FLUIDES GÉOPHYSIQUES. 11 1.4 Les fluides géophysiques. Les résultats que nous allons présenter sont centrés autour du problème de l'hydrodynamique du lac de Neuchâtel. Ce problème, bien que très particulier, revêt pourtant un intérêt "universel". Presque tous les fluides géophysiques présentent en effet la caractéristique suivante: Le rapport étendue horizontale/ profondeur est très petit. Région Diamètre horizontal d Hauteur h Rapport h/d Europe 5000 km 10 km 0.2% Atlantique Nord 5000 km 5 km 0.1% Lac de Neuchâtel 40 km 100 m 0.25% flaque d'eau 4 m 1 cm 0.25% Le lac de Neuchâtel est donc un excellent modèle réduit de la mer par exemple. Pour lui aussi, l'effet de Coriolis est très sensible (les marées en revanche sont imperceptibles), et nous disposons d'un appareillage permettant des mesures plus denses que pour l'Atlantique. En plus de ce caractère général des fluides géophysiques, la bathymétrie du lac révèle une colline sous-marine affleurant à 8m sous la surface, ce qui ajoute une difficulté à la simulation numérique et confère à notre lac une nature exemplaire pour l'hydrodynamique des fluides géophysiques [42]. Compte tenu des grandes variations de profondeur, une discré- tisation par éléments finis est préférable à une méthode de différences finies pour simuler plus précisément les phénomènes. Les premiers travaux effectuant l'analyse mathématique de l'hypothèse hydrostatique comme modèle asymptotique pour les fluides peu pro- fonds sont [11, 10]. Ils traitent respectivement du problème de Stokes et des équations de Navier-Stokes stationnaires. Un autre aspect important de l'hydrodynamique du lac concerne la stra- tification, c'est à dire le caractère inhomogène de la densité du fluide. A cause des phénomènes thermiques, une interface appelée thermocline sépare la couche de surface formée de fluide chaud plus léger de la partie inférieure formée de fluide froid, à 4 degré. Un modèle adéquat du lac devrait prendre en compte la densité variable [I]. Cela nous a amené à nous intéresser à la résolution de l'équation de conservation de la masse, c'est à dire de l'équation de transport linéaire. 12 CHAPITRE 1. MOTIVATIONS. 1.5 Plan de la thèse. Dans la première partie, nous nous intéressons aux équations de Stokes puis de Navier-Stokes en bassin peu profond, sous l'approximation hy- drostatique qui fournit un modèle asymptotique particulièrement per- tinent. Nous montrons que le problème de Stokes hydrostatique est bien posé, construisons un élément fini stable pour sa résolution numé- rique [5]. Nous prouvons ensuite que les équations de Navier-Stokes sous l'approximation hydrostatique admettent une solution faible au sens de Leray-Hopf [33, 26]. Malheureusement nous n'avons pas pu résoudre le cas non-linéaire, à cause du manque de contrôle sur la régu- larité de la vitesse verticale, nous traitons uniquement le cas linéarisé. Dans la seconde partie, nous traitons le problème de l'équation de trans- port. Nous avons développé, en collaboration avec P. Perrochet [44], une méthode simple et stable de résolution en éléments finis, basée sur une formulation parabolique de l'advection. Ceci est surprenant à tous égards, mais nous montrerons que grâce à la méthode des moindres carrés écrite dans l'espace-temps, l'équation de transport prend la forme d'une équation de diffusion [6]. Ceci nécessite cependant l'usage d'es- paces de Sobolev anisotropes et d'une inégalité de Poincaré "courbe" que nous prouvons [7]. Partie I L'approximation hydrostatique. 13 Chapitre 2 Le modèle asymptotique des fluides géophysiques. 2.1 La faible profondeur. La circulation des masses d'air en météorologie ou la circulation des eaux en océanographie et en limnologie sont décrites de manière géné- rale par les équations de Navier-Stokes. Pour prédire de manière satis- faisante le mouvement de ces fluides, les météorologues et les océano- graphes utilisent depuis longtemps des modèles asymptotiques pour les équations de Navier-Stokes (voir par exemple [42]). Ces modèles sont tous basés sur l'observation fondamentale suivante. Les dimensions horizontales du bassin à étudier sont beau- coup plus importantes que ses dimensions verticales. Dans tous les cas on remarque que le quotient t = profondeur / diamètre horizontal est de l'ordre de quelque pour mille (voir table dans le chapitre d'in- troduction). Tous les modèles asymptotiques actuels pour étudier le mouvement d'un fluide dans un tel contexte tiennent compte du fait que c est très petit. Le plus simple est l'approximation hydrostatique des équations de Navier-Stokes [11, 10]. Cette méthode utilise de ma- nière essentielle un modèle de viscosité turbulente pour le liquide, basé sur l'observation que la différence entre les échelles horizontale et ver- 15 16 CHAPITRE 2. LA MODÉLISATION ASYMPTOTIQUE. d Figure 2.1: Le beau lac de Bâle. ticale induit une turbulence différente dans les deux directions (voir [42]). 2.2 Équations gouvernant l'écoulement. Soit fi£ C R3 le domaine borné à frontière suffisamment régulière (di- sons C1 par morceaux) défini par fi£ = {x = (ii,I2, z); (Zi1I2) € T5, -£ • A(Xi1X2) < z < 0} où F, est un domaine borné de R2 représentant la surface du lac et h : Va —> R est une fonction C1 par morceaux. On note alors Ti = dùt \ T, le fond du bassin, où diïe désigne la frontière du domaine fit. Le domaine fit modélisant un lac ou une mer est occupé par un liquide dont le mouvement est engendré par une traction horizontale (induite par le vent) sur la surface T4 et ce mouvement est soumis à la force de Coriolis induite par la rotation de la Terre. On suppose la densité constante égale à l'unité. Avec les hypothèses et les notations précédentes, l'écoulement est régi par les équations de Navier-Stokes incompressible, avec viscosité anisotrope. v, +v Vv = A„v-2w Av- VP+ g dansO, (2.1) divv = 0 dansOt (2.2) v = 0 sur T6 (2.3) v,-£-¦= T1, I/*-r- = T2, U) = O suri, (2.4) OZ OZ v(.,t = 0)=0 dansO£. (2.5) 2.3. CHANGEMENT D'ÉCHELLE ET RENORMALISATION. 17 V = (U11U21Iu) Dans Eq. (2.1) A„ est défini par A„ = v\-gpi + "2§f-ï + ^ Jp-- où v = (ui,u2,u;) représente la vitesse du liquide, u = (O1O1U3) la vitesse angulaire du bassin, P est la pression totale, g = (0,0,— g) est l'accélération de la pesanteur, u = (v\,vi',i/t) le vecteur de viscosité cinématique turbulente. Dans l'équation (2.4) r„ i = 1,2 représentent les tractions horizontales exercées par le vent sur la surface T, du liquide. Remarque. On a de plus imposé w = 0 sur T5 pour éviter la résolution d'un problème à frontière libre, qui ne serait pas justifié pour le lac de Neuchâtel. En effet les seiches les plus hautes sont inférieures à 5 cm, ainsi qu'en t'emoignent les mesures effectuées par Léon du Pasquier [47]. D Pour rendre le domaine fît indépendant de e, effectuons les 2.3 Changement d'échelle et renormalisa- tion. fi( = {x = (I11I21J); (I11I2) € ra, -e • A(I11I2) < z < 0} On pose alors: z = £ • X3 W = £ ¦ U3 U = (uuUi,U3) Cl= {x = (I11I21I3); (I11I2) 6 r„, -Zi(I11I2) < I3 < 0} On suppose de plus que: vz = t2 ¦ u3 T, =t-9i, i = 1,2 Nous effectuons cette renormalisation pour des motifs physiques et aussi des impératifs mathématiques: • Pour le vent, il est naturel de le supposer proportionnel à la pro- fondeur, sinon lorsqu't —? 0, on simulerait un ouragan! 18 CHAPITRE 2. LA MODÉLISATION ASYMPTOTIQUE. m Quant à la viscosité, elle contrôle la taille des tourbillons, et au vu de la faible profondeur, il est normal de supposer vz très petit. Une dépendance quadratique en t est nécessaire pour aboutir à l'approximation hydrostatique, comme le montrent Besson-Laydi [10]. Posons enfin p = P + g ¦ z = P + g ¦ ex3 la pression hydrodynamique. 2.4 L'hypothèse hydrostatique. Avec ces changement d'échelle et renormalisation, les équations (2.1- 2.5) deviennent: -^- + u • Vu1 - A„ui - Ju2 + -TT- = 0 dans fi (2.6) Ot OXi -p- + u • Vu2 - A„u2 + /U1 + -^- = 0 dans fi (2.7) Ot 0X2 £M^r + u -Vu3- A„u3} + P- = 0 dans fi (2.8) Ot OX3 divu = 0 dans O (2.9) u = 0 sur T6 (2.10) «*£^ = *i. "3^ = 02' "3 = 0 surT, (2.11) u(.,i = 0) = 0 dans fi (2.12) L'approximation hydrostatique consiste alors à négliger les termes en e2 pour aboutir au modèle suivant: -^- + u • Vu1 - A„ui - /u2 + —^- = 0 dans fi (2.13) Ot OX1 -^ + U-Vu2- A„u2 + /u, + —^- = 0 dans Ji (2.14) Ot OX2 -^- = 0 dansfî (2.15) 0x3 divu = 0 dans fi (2.16) U1 = U2 = u3n3 = 0 sur Tt, (2.17) "3^- = 01, V3^1 = 02, u3 = 0 surT. (2.18) Ox3 0x3 u,(.,t = 0) = u2(.,t = 0) = 0 dansfî (2.19) 2.4. L'HYPOTHÈSE HYDROSTATIQUE. 19 L'approximation hydrostatique des équations de Stokes est étudiée au chapitre 3, et celle des équations de Navier-Stokes instationnàires linéa- risées, avec effet de Coriolis, est analysée au chapitre 4. 20 CHAPITRE 2. LA MODÉLISATION ASYMPTOTIQUE. Chapitre 3 Le problème de Stokes hydrostatique. Dans ce chapitre nous étudions le problème de Stokes sous l'approxima- tion hydrostatique, nous montrons qu'il est bien posé et nous construi- sons un élément fini stable, c'est à dire n'exhibant pas de modes para- sites de pression, ni de vitesses verticales. 3.1 Préliminaire d'analyse fonctionnelle. Dans cette section, nous donnons des résultats qui, bien qu'essentielle- ment classiques, mériteraient d'être mieux connus. Il s'agit en fait de la généralisation du théorème de Lax-Milgram au cas non-coercitif, ou encore de la généralisation en dimension infinie de la notion de forme bilinéaire non dégénérée, qui est la bonne propriété caractérisant un système linéaire bien posé, le cas défini positif n'étant qu'un cas particu- lier. Alors que le théorème de Lax-Milgram ne donne qu'une condition suffisante, le théorème suivant donne une condition nécessaire et suf- fisante pour qu'un problème variationnel soit bien posé. Au préalable définissons la notion généralisant la coercitivité ou "V-ellipticité": Définition 1 Soit U un espace de Banach, V un espace de Banach réflexif. Une forme bilinéaire a sur U x V est dite U-non dégénérée si (i) l'application v i-¥ a(.,v) est injective 21 22CHAPITRE 3. LE PROBLÈME DE STOKES HYDROSTATIQUE. (ii) il existe une constante C > 0 telle que a{u,v) Remarque. Dans le cas symétrique, (i) est superflue. Nous dirons alors que a est persuasive, plus français que faiblement coercive, (voir [4I]) D V désignant le dual topologique de V, considérons le problème va- riationnel : (P) Etant donnée / € V, trouver u e U tel que a(u,v) =< f,v > Vw € V Théorème 1 Soit a bilinéaire continue sur U x V, le problème (P) est bien posé si et seulement si a est U-non dégénérée. Preuve. Il suffit de modifier légèrement la preuve de Necas [39], voir également Oden-Reddy [41] ou Roberts-Thomas [46]. Nous donnons la preuve cependant, car les références citées ne traitent pas le problème dans toute sa généralité. Soit l'opérateur A: U -+ V u t-> a(u,.) a étant bilinéaire continue, A est un opérateur linéaire continu. (P) est bien posé si et seulement si A est un isomorphisme, c'est à dire injectif, surjectif et ouvert. D'après le théorème de Banach, l'ouverture est conséquence de la surjectivité. D'autre part A est surjectif si et seulement si son image imA est dense et fermée. Or kn~A = Ker(A*)L, où A* désigne l'adjoint de A. Donc nous avons établi le Lemme 1 (P) est bien posé si et seulement si A injectif, imA fermée et A* injectif. 3.J. PRÉLIMINAIRE D'ANALYSE FONCTIONNELLE. 23 Nous pouvons maintenant prouver le théorème. =>: V étant réflexif, A*: V -> U' v i-¥ a(.,v) L'injectivité de A* est donc exactement la condition (i). Par le lemme 1 nous avons d' autre part A injectif et imA fermée, donc imA est un espace de Banach et l'opérateur (encore noté A abusivement) A : U -> imA est un isomorphisme, grâce encore au théorème de Banach. Son inverse est donc continu, il existe ainsi une constante C > 0 telle que WM > CIMI or par définition < Au,v> a{u,v) WAu\\ = sup-----rr-7.-----= sup -t-tt- v/o IMI ,*o IMI ce qui fournit la condition (ii). ¦*=: Si a est {/-non dégénérée, en particulier a est non dégénérée. En effet la condition (ii) entraîne que Vv € V, a(u,v) = 0=*-u = 0 donc A et A* sont injectifs. En vertu du lemme 1 il reste à prouver que imA est fermée. Nous avons vu précédemment que la condition (ii) traduit \\Au\\ > C\\u\\. (3.1) Soit maintenant une suite (Aun)„ telle que lim Aun = v* n-*oo La suite (Aun)n est de Cauchy donc par (3.1) la suite (unj, également donc, comme U est complet, il existe u 6 U tel que lim Un = u n-»oo Enfin par continuité de A, v* = Au € imA. ¦ 24CHAPITRE 3. LE PROBLÈME DE STOKES HYDROSTATIQUE. Remarque. I.Babuska [8] et F.Brezzi [13] ont prouvé des théorèmes analogues pour des espaces de Hilbert, mais sous les deux hypothèses : sup^l > CNI et sup ^l > CHI. */o IMI u*o IMI Soit Uh (resp.VJ,) sous-espace fermé de U (resp.V) et considérons le problème variationnel discret associé à (P) : (Ph) Etant donné / G V/, trouver Uh € Uh tel que a(uh,Vh) =< /, u/, > w„ e Vh Théorème 2 Soit a une forme bilinéaire continue sur UxV1 U-non dégénérée et telle que sa restriction à Uh x VJ1 soit Uh-non dégénérée. Soit u (resp.Uh) solution de (P) (resp.ÇP^)). On a l'estimation d'er- reur : \\u-Uh\\<(l + M/Ch) inf ||u — wjkll où les constantes C (resp.Ch) proviennent des inégalités de stabilité de la définition 3.1 dans les espaces U, V (resp. Uh, Vh) et M désigne la norme de a. Ce théorème est dû à I.Babuâka [8]. Remarque. Si Ch est indépendante de h, le théorème précédent fournit une estimation optimale de l'erreur ||u — uj,||. D 3.2 Analyse du problème continu. 3.2.1 Equations du modèle. Soit f! C R3 le domaine à frontière suffisamment régulière (disons C1 par morceaux) défini par n = {i = (Xi) e R3, (ii,X2) e r„ -/i(i,,i2) < i3 < o} où T5 est un domaine borné de R2 représentant la surface du bassin et h : T5 —f R+ est une fonction C1 par morceaux. Le domaine fl 3.2. ANALYSEDU PROBLÈME CONTINU. 25 (représentant un lac ou un océan) est occupé par un liquide incompres- sible en mouvement lent et stationnaire. Sous la condition naturelle en géophysique que la profondeur est faible devant l'étendue, on fait l'hypothèse hydrostatique [42] ^ = 0 où p désigne la pression dynamique, i.e. débarrassée de l'influence de la pesanteur. Notant dorénavant ô, = gf-, les équations du problème, où nous prenons des conditions de Dirichlet homogènes pour simplifier la présentation, sont alors [11] -Au, + Oj p = /i dans fi (3.2) -Au2+ d2p =/2 dans fi (3.3) d3 p = 0 dans fi (3.4) divu = 0 dans fi (3.5) U1 = U2 =0 sur dû (3.6) U3-Ii3 =0 surdfi (3.7) où n désigne Ia normale sortante sur dû. Remarque. Par rapport au problème de Stokes usuel, la disparition du terme -Au3 de (3.4) confère à (3.2)-(3.7) un caractère dégénéré. La vitesse verticale U3 est donc à priori moins régulière que les vitesses horizontales U\ et U2, c'est pourquoi la condition limite (3.7) diffère de (3.6). Nous verrons dans la section suivante à quels espaces appartien- nent les vitesses horizontales et verticales. O 3.2.2 Formulation variationnelle du problème. Nous désignons par (suivant [16], [57]) l?(Û) l'espace de Hilbert des fonctions de carré intégrable sur fi muni du produit scalaire usuel noté (.,.) et de sa norme associée notée II ||o qui pourront également dénoter les produit scalaire et norme produit sur L2(û)" avec n convenable. 26CHAPITRE 3. LE PROBLÈME DE STOKES HYDROSTATIQUE. H'(Q) l'espace de Sobolev d'ordre s muni de sa norme usuelle notée II II j H0-(Q) = {

(Vu2, Vv2) - (p, d2v2) = < /2, U2 > -(PtO3V3) = -(divu,g) = (3.8) ou encore en rassemblant (P) Trouver u € H tel que a(u, v) =< /, u > Vu £ W où / £ //' donné. avec a(u,u) = (Vu15Vu1) + (Vu21Vu2) - (p,divu) - (divu,q), 4 1=1 Remarque. Le problème (3.2)-(3.7) est en fait un cas particulier du problème (3.8), où /3 et /4 sont nuls. D 3.2. ANALYSE DU PROBLEME CONTINU. 27 La formulation faible (P) traduit un principe variationnel mixte qui s'énonce Théorème 3 L'élément u — (U11U21U31P) Ç H est une solution de (P) si et seulement si u est point selle de la fonctionnelle £(ti„u2,U3,p) = 1/2(1IVm1HS+ ||Vu2||2)-(p,divti) - < /l,«l > - < /2,"2 > - < /3,"3 > - < ft,P> en ce sens que pour tout (vi,v2,v3,q) € H C(uuu2,u3,q) < £{uuu2,u3,p) < C{vuv2,u3,p) C(uuu2,v3,p) < C(uuu2,u3,p) < C(vuv2,u3,p). Preuve. Soit u = (ui,u2,u3,p) € H, et un accroissement Su = ((5Ui1(JUj1(JuS1(Jg). On calcule, en utilisant la bilinéarité symétrie de a. SC = C(u + Su) - £(u) = a(u, Su)- < f,Su> +-a(Su, Su) (3.9) =>: Supposons u solution de (P). Alors, (3.9) devient SC = -a(Su,Su) Il suffit alors de prendre des accroissements particuliers. Pour (Ju = (0,0,0,9— p), cela donne £(ui,U2,U3,¢) = ) + t2a(v,v)/2. La condition de point selle entraîne la stationnante de cette fonction SC = 0 de i, donc la nullité de a(u,v)— < /, u >. ¦ Remarque. La pression p (resp. la vitesse verticale U3) s'interprète comme le multiplicateur de Lagrange associé à la contrainte —divu = /4 (resp. d3p = J3). a 28CHAPITRE 3. LE PROBLÈME DE STOKES HYDROSTATIQUE. 3.2.3 Existence, unicité, stabilité de la solution. Dans ce paragraphe, nous démontrons le résultat suivant (voir aussi [H]) Théorème 4 Le problème de Stokes hydrostatique (3.8) est bien posé Preuve. D'après le théorème 1 et vu la symétrie de a, il s'agit de prouver que o(.,.) est une forme bilinéaire continue persuasive. La bilinéarité est triviale, la continuité découle aisément de ||di toi + O2W2Ho ^ HVHIo, Vtü = («"t.«*) e HKa) * H^tI). (3.10) Il reste à vérifier la condition de stabilité ^ > CII-II (3.11) sup v/0 ||t> où C désigne une constante strictement positive. Fixons u = (ui,u2,U3,p) € H et construisons w e H tel que a(u,w) > C||ii||||u)||. La condition inf-sup continue (voir [23]) en dimension 3 »P ^r1SCiMi. (312) *etf,j(n)3-{o} Il V\U permet de choisir v 6 Hg (fi)3 tel que (p,div«)> BiHiolMI, (3-i3) avec B constante convenable. Par homogénéité on peut prendre v tel que NI. = HpIIo- (3.14) Puis on choisit q tel que ^^>IMo, (3.15) ll9l|o pour cela prenons q = d3u3 (3.16) 3.2. ANALYSE DU PROBLEME CONTINU. 29 qui est bien dans Ll(U) vu la condition limite (3.7). Cherchons w de la forme Oj(U11U21U3, -p) - ß{vi,v2,v3,0) -7(0,0,0,g). On a a(u,w) = a(||Vu,||S +HVtI2IlS) -/3((Vu1, V«,) + (Vu21Vv2) - (p.divv)) +7((Oi Ui + O2 U2, q) + (¾ U3,9)) que l'on peut minorer à l'aide des inégalités (3.10), (3.13) de Cauchy- Schwarz, et des égalités (3.14), (3.16) a(u,w) > 0(1IVu1IIS+IIVu2IIS) -0IIpIIo(I|Vu,||S + 1IVu2IiS)1/2 + /3BIIpIIS +71Ia3U3IIS-IIIO3U3IIO(IIVu1IIS + IiVu2IiS)1/2. Pour conclure on utilise le lemme élémentaire Lemme 2 // existe a, ß, 7 € R et C > 0 tels que V(x, y, î)eRJ, Qi2 + ßBy2 + 722 - ßxy - 712 > C(x2 + y2 + z2). Preuve. Comme en dimension finie toutes les normes sont équivalentes, il suffit de s'assurer que la forme quadratique ax2 + ßBy2 + 722 — ßxy — 712 est définie positive. Pour cela, ajustons a,ß et 7 pour que les mineurs principaux de la matrice associée soient strictement positifs . On critère bien connu de constate que q = 1,/3 = B, y = 1 jouent. ¦ définie positività [55] En remplaçant alors x (resp. 3/,2) par (1IVu1HS + ll^^llo)1/2 (resp. HpIIo, 11¾ U3IIo) on obtient U(U1U1)^C(IIVU1IIS + IIVU2IIS + IIo3U3IIS + IIpII2). (3.17) Or d'après la forme de w Hl < en (3.18) d'où le résultat. ¦ 30CHAPITRE 3. LE PROBLÈME DE STOKES HYDROSTATIQUE. 3.3 Résolution du problème approché. 3.3.1 Critères pour une discrétisation bien posée. Soit Hh = Uh x Ph où Uh = Vh x Wh et 14 (resp. Wh, Ph) désigne le sous-espace de dimension finie des vitesses "horizontales" v\, u2 (resp. vitesses "verticales" U3 et pressions). En reprenant la preuve précé- dente, on constate qu'elle repose uniquement sur les inégalités (3.15) et (3.13). De plus, la démonstration est encore valable si (3.15) est remplacée par IMIo 1 Ce qui fournit le Théorème 5 Le problème de Stokes hydrostatique discret correspon- dant est bien posé si et seulement si sont vérifiées - l'inégalité de Brezzi (voir [23], [14]) (ph,divvh) SUP ~~iïïTiï— - CWI°> KhEV)1-(O) ||«ft||l - l'inégalité "hydrostatique" (d3u3h,qh) ...J, ., SUP ll„ M > 11¾ MIo- ?<.<=n-{o) libilo Jusqu'à très récemment, les techniques de vérification de l'inégalité de Brezzi étaient disparates et délicates (voir [23], [14]). En revanche, la technique de macro-élément de Stenberg [53], [54] est basée sur un principe simple: un macro-élément M est la réunion d'un ou plusieurs éléments adjacents par une face, une arête ou un sommet. Deux macro-éléments M et M' sont dits équivalents s'ils sont homologues par une transformation continue F appliquant les éléments internes de l'un sur ceux de l'autre de façon isoparamétrique. Soient Ch un maillage formés d'éléments A', Tj1 la collection des faces des éléments de Ch contenues dans l'intérieur du domaine H et h = maxAgc,, {diarn(K)}. Notons 1Il suffit alors de prendre 7 = A i dans le lemme. 3.3. RÉSOLUTION DU PROBLÈME APPROCHÉ. 31 Pm = {p\m |p€ Ph) Vo,m = {v\m \v e Uh, v = 0 suröM} Nm = {p e Pm \ (d\vv,p)M = O, Vu € V0,M}. Rolf Stenberg prouve alors le Théorème 6 Supposons que le maillage satisfasse l'hypothèse de ré- gularité inverse diam(K) > Ch, VA' £ Cj1. Supposons qu'il existe une collection finie de classes d'équivalence £,, i = 1,..., N de macro- éléments, un entier L et un ensemble Mh de macro-éléments tels que (Ml) pour tout macro-élément M, l'espace Nm est de dimension un et se réduit aux fonctions constantes (MS) chaque macro-élément M € Mh appartient à l'une des classes Si, I-=I,..., JV (MS) chaque élément K G Ch est contenu dans au moins un et au plus L macro-éléments de Mh (M4) chaque T Ç Tj, est contenu dans l'intérieur d'au moins un et d'au plus L macro-éléments de Mh alors l'inégalité de Brezzi est satisfaite. Remarque. La seule hypothèse vraiment contraignante est (Ml), pour (M2)-(M4) il suffit de construire un maillage qui les satisfait. L'hypo- thèse (Ml) s'interprète physiquement, elle garantit l'absence de modes parasites locaux depression (spurious modes ou checkerboard pressures, voir [14]). La puissance virtuelle des forces de pression (divu,p) = — < Vp, v > ne devant être identiquement nulle que si la pression est constante, ü 32CHAPITRE 3. LE PROBLÈME DE STOKES HYDROSTATIQUE. noeuds V1, v2 Q noeuds V3 |___J noeuds p Eléments R2/R1/R0.5 Figure 3.1: Eléments hydrostatiques de référence 3.3. RÉSOLUTION DU PROBLÈME APPROCHÉ. 33 3.3.2 Elément fini hydrostatique. Considérons un maillage Ch dont les éléments sont des prismes (resp. des hexaèdres) verticaux, base et toit de ceux-ci n'étant pas nécessaire- ment parallèles. Rt(X1,X1,-. .,In) désigne P* ou Q* (voir Ciarlet [16]) selon que les éléments sont prismatiques ou hexaédriques. Le produit tensoriel des espaces vectoriels E et F est noté E ® F. Soit Fk l'appli- cation transformant l'élément de référence K en l'élément K de façon isoparamétrique. Pour x Ç K notons f(x) = /(Fa(x)). Soit Cfc/2 le maillage obtenu en subdivisant chaque élément de C/, en huit, en coupant chaque arête en son milieu. Les espaces des vitesses sont vk = {<è€ c(A)2n//,Un)2I va-e c/2, ^ €(*,(£;,*;) ®p,(ïs))2}, wh = {x e c(n)n/Z0(O3,n)| va- e ck, xjî e a,(£;,ï;)® P1OS)J. Pour la pression, on ne demande pas la pleine continuité, on exige seule- ment la continuité aux interfaces verticales : pour toute F face verti- cale commune à deux éléments A'j et A'2, Va G F, \\mI-ta%xeKl p(x) — l\mx^a,xeKj p(x) = p(a). L'espace des pressions est Ph = { V-€ Ij(O)I VA'6 C*, ^-e Ri(xi,x-2) ® P0(Si) et 0 continue aux interfaces verticales }. Les degrés de liberté sont décrits sur 1' élément de référence (voir fig. 3.1). Remarque. Au lieu de découper les éléments en huit, pour les vitesses horizontales on peut prendre vh = { e C(A)2 n //0'(n)2| va' e C, ^- e (a2(îï,X2^p2(S3))2). D Nous donnons en annexe un programme FORTRAN de l'élément fini de pression pour l'hexaèdre. UCHAPlTRE 3. LE PROBLÈME DE STOKES HYDROSTATIQUE. propriété caractéristique de ce que nous appelons un élément fini hydrostatique. 3.3.3 Vérification des inégalités de stabilité. Dans la suite, nous supposons que le maillage Ch est régulier et qu'il satisfait l'hypothèse inverse "verticale" £""""- > n > O1 où 2/îmm (resp. 2/imai) désigne la plus petite (resp. grande) arête verticale des élé- ments de Ch- Nous donnons les preuves pour l'élément prismatique, elles s'étendent aisément au cas hexaédrique. Proposition 1 L'élément fini hydrostatique vérifie l'inégalité hydro- statique. Preuve. L'idée de base est que si U3 € Wh alors 63 U3 est "presque" dans Ph, En effet, soient ZJ11 B[, B" les points de K homologues aux Ai, A\,A", i = 1,2,3 de K (voir fig. 3.1); comme les éléments sont tous verticaux, BjBj' = 2Zi^e3, i = 1,2,3, par dérivation composée on a alors du3 du3 Ox3 Ox3 où /ia(xi,X2) = (1 — x"i — x2)/io + xihi -f x"2/i2. Ainsi min(Zio,ht Zìa- < ITIaX(ZiO1Zi11Zi2). Soit u3 € Wh, sur l'élément K prenons -7-,. dûs 9'*(x) = ôxy on voit aisément que ¢6 Ph, en particulier q est continue aux interfaces verticales. Minorons (3.19) M< (3.20) (d3 U3, q) = Y, / . ^3 "3 Q àx d3 u3 \\l .||&«s|| (3.21) (3.22) par (3.19) et (3.20) (d3u3, q) > h„ IkIIo < h„ De (3.21) et (3.22) on déduit (d3u3,q) hmi„ ||9||0 hmax L'hypothèse de régularité "verticale" du maillage garantit enfin que jf"^ > C constante strictement positive indépendante de la taille du maillage. ¦ 3.3. RÉSOLUTION DU PROBLÈME APPROCHÉ. 35 Remarque. On peut voir que l'élément de Hood-Taylor R2/Ri, ainsi que l'élément de Bercovier-Pironneau 8Ri/Ri ne vérifient pas l'inégalité hydrostatique et peuvent donner lieu à des vitesses verticales parasites. D Pour chaque sommet So intérieur à H, formons le macro-élément M composé de tous les éléments le contenant (voir fig. 3.2). Nous suppo- sons aussi que le maillage Ch considéré vérifie les conditions (M2)-(M4) du théorème 6, qui sont aisées à respecter dans la pratique, ainsi que 1' hypothèse de régularité. Il reste à vérifier (Ml) pour prouver la Proposition 2 L'élément fini hydrostatique vérifie l'inégalité de Brezzi. Preuve. Faisons la preuve pour le macro esquissé fig.3.2, les autres cas sont analogues. Notons p(Ai), p(Bi), i = 0,...,5 les degrés de liberté de p sur les douze arêtes verticales du macro-élément. Il s'agit de prouver Vu € Vo1M, (div,p)M = 0 =*• p = constante (3.23) (d'ivv, p)m = — 5Z / ^7P ' v^x + / pu • nds KcMJK JBK à l'aide de la formule de Green. Comme p est continue aux interfaces verticales, les intégrales de bord correspondantes se détruisent mutuel- lement. De plus U = O sur les parois externes, donc ne subsistent que les intégrales sur le plancher du premier étage (contenant So) (divu,p)M = - Yl / Vpudx + 5Z / [p+-p_]u' L-mJK „ Ja nds où p_(resp. p+) concerne le rez de chaussé (resp. le premier étage), n désigne la normale dirigée vers le bas et a les faces composant le plancher du premier étage. Il reste à choisir des vitesses virtuelles révélatrices. Soit u nul sur tous les sommets sauf en So où u = (0,0,1) On a Vp-v = 0 car (¾ p = 0 donc (divu, p)m = H / [P+ - P-K • n3ds j=0 J"l mCEAPlTRE 3. LE PROBLÈME DE STOKES HYDROSTATIQUE. B1 \/ s J \ ^ A5, y y y y y S5 4 1S N I5" < y V, / S0 y* B5 y y y y y -—rr/ y I Bq J3 \ y Bc Figure 3.2: Exemple de macro-élément 3.3. RÉSOLUTION DU PROBLÈME APPROCHÉ. 37 où (Tj désigne la face (So,Sj,Sj+i). En utilisant une formule de qua- drature, on obtient 5 (divu,p)M = 53Wj(P(Aj) -P(Bj)I Lo1 > 0. (3.24) J=O Puis on choisit v nul sur tous les sommets sauf en Ij où v = AoAj. On a (divu,p)jw = / Vp-udx-f / Vp ¦ udx = C(pj — po). Jk3-I Jk3 où Kj désigne l'élément posé sur la face cr,-, et C > 0. D'après (3.23) ceci doit être nul donc P(Aj) = PMo) i = 0,...,5. En travaillant sur le premier étage, on obtient de même P[Bi) = P(B0), i = 0,...,5. Enfin en portant dans (3.24) on obtient p = constante Des propositions 1, 2 et du théorème 5 on déduit le Théorème 7 La discrétisation par éléments finis hydrostatiques de (3.8) est bien posée. 3.3.4 Estimation d'erreur. Introduisons les espaces de Sobolev anisotropes Hk(d3,il) = {w€ ff*-'(n)| d3we H*'1 (il)} D'après le théorème 2, nous pouvons assurer le Théorème 8 On a l'estimation d'erreur ||«J -Ulftlll + ||U2 — U2fc||l + ||£feU3-â3U3k||o+ ||p-p*||o < C inf { IK - £3, rendant l'hypothèse (4.16) superflue mais il n'en est rien. On peut en effet construire \j> dans #0(¾), sans (xux2,x3) = (x\,x2) ¦ x3(/i(xi,Z2) - X3), où <£(xi,X2) est dans L2(ù) mais pas dans L3(fi). O On se donne également ¢,,¢,€!"(0,7-^-"2^.)). (4.17) Remarque. Les fonctions O1, O2 correspondent aux tractions horizon- tales du vent sur la surface du lac, il est donc naturel de les supposer essentiellement bornées. O Même en cas d'oura- gan. Remarque. On a évidemment V <-> W avec injection dense. En effet Cf(Q) x Cf(fi) x ©(0) C V. D Remarque. Nous avons donné à a la régularité minimale, afin d'englober autant que faire se peut le cas non-linéaire. Q Nous aurons besoin du Lemme 3 / (a • Vu1)Vj = - / Uj(aV - Vj) Vu, v6ff|x H\ x H0(O3) Ja Jn Preuve. Il suffit d'intégrer par parties en utilisant que a • Vuj = div(uja) grâce à la divergence nulle de a et les conditions limites au bord du domaine. ¦ 42 CHAPITRE 4. ÉCOULEMENT EN EAU PEU PROFONDE. Iu1 O ON Rappelant (chap. 2) que A1^ = div(i/Vy>), où v = J 0 U1 0 \ 0 0 u3 ) définissons naturellement V'„p = t/1/2VvJ. Une simple intégration par parties permet d'établir le Lemme 4 Soient u;,u, 6 //£, i = 1,2 vérifiant 1/383 u; = ¢,-, i = 1,2, alors (-A„u,,i>i) = (V„u;, V„u,)- < 0„u, >w-i/2(pi)iHi/2(ri) . 4.4 Forme faible renormalisée. Supposons que (4.1-4.7) admette une solution suffisamment régulière, en intégrant par parties, en utilisant les conditions limites, initiales et les lemmes précédents, on obtient la formulation faible suivante: Etant donnés O1, A2 6 L°°(0,T; /T1'2^,)), trouver u € £2(0,T; V) rT dv, rT rT rT / (u,, -757) dt+/ («,, a- Vu1) dt+/- / (u2, U1) dt-/ (V„u,,V„t;,)dt Jo ai ^o jo jo + J0 (u2,-^)dt+y (u2,a-Vuj)dt-/-y (ui,u2)dt-jf (V„u2,V„t>2)dt +£2{j[ (u3."^)dt + / (U3, a ¦ Vu3) dt - I (V1^u3, V„u3)dt} = -/ <ôuv{> dt- f < ô2,u2 >dt (4.18) Jo jo Vv = (U11U21U3) 6 C(O1T; V); v(T) = 0. Il est bien connu [32]-[35]-[58] que (4.18) admet une solution u = u' pour chaque e > 0. De plus u dt-/ <6»2,i>2>dt (4.19) 7o Jo Vv = («,,Vj1W3) e C(O1T; W); v(T) = 0. Vérifions que tous les termes figurant dans la formulation ci-dessus ont bien un sens. Les seuls posant problème sont rT / (u,,a- Vt>i), dt i = l,2. Jo Or par le théorème d'injection de Sobolev Ui e L2(0, T; //'(H)) -> L2(0,T; L6(fi)), pour la même raison aj e L2(0,T;H\n)) ^ L2(0,T;L3(U)), j = 1,2 l'hypothèse (4.16) assurant la même qualité à a3, d'autre part t>i, V2 € C(O1T;//1 x //') donc Vu1, Vv2 € L°°(0,T; L2(ü)), l'inégalité de Helder donne alors: /T(u,-, a • Vu1) dt < fT 1Iu1II6IIaII3IIVu1H2 dt Jo Jo < ||«i|Uî(o,r;t'(n))HaIlL2Io1^LS(O))IlVuiHioojo.TiL^in)). Le paragraphe suivant est consacré à l'étude de la convergence lorsque e —> 0 d'une solution de (4.18) vers une solution du problème (4.19) . Nous utilisons la méthode de compacité [35], basée sur des estimations à priori d'énergie. Nous aurons besoin pour cela d'une autre formulation faible, déduite de (4.18), que nous présentons ci-après. 44 CHAPITRE 4. ÉCOULEMENT EN EAU PEU PROFONDE. 4.6 Majorations à priori. 4.6.1 Forme faible spatiale. En prenant des fonctions tests de la forme {t)v{x),teV{o,T),vev dans (4.18), on obtient la formulation suivante. Etant donnés ¢,, O2 e L°°(0,T; H-1^(T,)), trouver u G L2(0,T; V) tel que ^1'"!) (u„a- V«,) - / • («2,«,) + (v„u„ V„t>,) dt d(u2,v2) dt (u2, a • Vu2) + / • (U1, V2) + (V„u2, V1^2) +£2. {^3p) _ Ka . Vws) + (Vl/u„ V,«,)} =< O11U1 > + <ô2,u2>, dans #(0, T) Vv = (vuv2,v3) e V. (4.20) u = u£ désigne dans ce qui suit une solution du problème renormalisé. 4.6.2 Égalité d'énergie. Prenant v := u dans (4.20), par symétrie les termes de Coriolis dispa- raissent, du fait de la divergence nulle les termes non-linéaires également donc il reste Hu1II2 + ||u2||2 + £2||u3||2) 2 dt + IIV.UjH2 + ||v„u2||a + E2IIV11U3II2 = / U1O, ds / u202ds (4.21) Jr, Jr, Remarque. Ceci s'interprète comme un bilan d'énergie ': Energie ciné- tique + Energie dissipée par viscosité = Energie fournie par le vent. D Nous allons prouver la Proposition 3 U11U21CU3 sont bornées dans Z,°°(0,T; L2(Çl)), et dans L2(0,T;H\iï)). Ml s'agit en fait de densité de puissance 4.6. MAJORATIONS A PRIORI. 45 Preuve. Grâce à (4.17), on majore le second membre de (4.21) par INI-l/2||«l||l/»+l|02||-l/»|M,/2, D'autre part l'inégalité de Poincaré verticale, applicable puisque u, = 0 sur Ti,, et la continuité de la trace donnent IKIIi/ï0 (4.23) avec ij = v3/C pour majorer le second membre de (4.22) et obtenir \jt (IhII2 + IhII2 + <2lhl|2) + ",(lift«,us + liftu2||S + £2 lift U3IiS) +^(Hftu.iiS+iiftu.iiS+^iift^i^+^diftu.iis+iiftu.iiS+^iift^ilS) (u*,u*2,w*) dans L^(O,T\L2(Q)) faible - *. (4.26) (u\,u\,eu3) -^ (u,, U2, w) dans L2(0,T; H^ x Hl x H^) faible. (4.27) u3 -- U3 dans L2(0, T; H0(d3)) faible. (4.28) Puisque L2(0,T; H^iI)) C L1 (0, T; L2(U)), de (4.26)-(4.28) on tire K, u;, w*) = (u,, u2, u;) G L°°(0, T; L2(ü)) D L2(0, T; V). D'autre part, comme divu' = 0, Vt, on a divu = 0. Donc u € L2(0,r;W). Remarque. D'après le corollaire 1, u| est bornée dans L2(0, T; L2(Q)), donc 6U3-f OdansZ,2(0,T;L2(fi)). Avec (4.26-4.27), il vient w = 0. D Considérons alors la forme faible (4.18), et fixons v = (»„»2, ws) £ C(O1T5V); v(T) = 0. Examinons la convergence de chacun des termes composant (4.18) lorsque <:->0. Par (4.26): rT ftit tT gv. j( <"">»-/>•>' .--,,2 De même: / (u2,vi)dt -> / (u2, D1) dt, Jo Jo / (uj,u2)dt -> / (U1, U2) dt. Jo Jo Par (4.27): / (V„u<> V1^1) dt -»• / (V„u,-, V17Ui) dt » = 1,2 Jo Jo 48 CHAPITRE 4. ÉCOULEMENT EN EAU PEU PROFONDE. E2 / (Vl/uc3,Vl,v3)dt = c- [ (V11CU3,V„u3)dt -> O. JO JO Pour les termes faisant intervenir l'advection, nous supposerons de plus Vi G C(O x [O1T]) i = 1,2,3, de sorte que a • Vu; G L1 (O1T; L2(U)). la régularité des Par (4.26) à nouveau: fonctions test est à .j .T discrétion. / (u|,a-Vt>;)dt-> / (u;,a • Vt);)dt, i = 1,2 JO JO De même: e2- / (u3,a-Vu3)dt = e- / (tu3, a • Vu3) dt -*¦ 0. On conclut enfin par la densité de C(H x [0,T])3 n C(O1TjV) dans C(O1TjW) que (ui,u2)u3) G L2(0,T;W) est une solution faible du problème hydrostatique (4.19). En effet C(O, T; Cf (fi) x Cf(Jl) x P(Sl)) C C(O X[O1T])3OC(O1T; V). 4.8 Obtention de la pression. Nous utilisons le lemme suivant sur les distributions à valeurs vecto- rielles dû à J. Simon [50] qui généralise le fameux lemme de G. de Rham [19]-[58]: Lemme 5 Soit g G X>\]0,T[; (tf-'(fl))3) tel que <9,v >(tf->(fl))3,(tf,>(n))3= 0 Vu 6 (H0-(Sl))3, div« = 0. Alors il existe q G V(]0, T[; L2(fl)) tel que 9 = V9. Si de plus pour s € R11 < r < oo, on a g G W(O1T; (tf-'(fl))3) alors on peut choisir q G W"-r(0,T; L2(ù)). De ce lemme nous allons déduire le Théorème 10 II existe p G W-^""(0,T; L* (Sl)) telle que -^- + a • Vu1 - A„u, - /u2 + 7— = 0 dans fl dt dxi -^- + a • Vu2 - A„u2 -r Ju1 + -- - 0 dans fl Ot OX2 -—- = 0 dans fl OJJ3 4.8. OBTENTION DE LA PRESSION. 49 Preuve. Soit v = (i>i,u2,i>3) 6 (//0(H))3, divv = 0. En prenant des fonctions tests de la forme *(t)v(*),*€ U(O1T) dans (4.19) on obtient la forme faible spatiale suivante d{uuvi) dt (U1,a- V Wi) - f ¦ (u2,vi) + (V„uuVuvi) + ';"2) - (u2, a • Vv2) + / • (u,, D2) + (V„u2, V„t>2) CEI = 0 dans V(O, T). (4.29) Posons Sf = {-Qf + a-Vu, - A„ui -/u2,-7r-î + a • Vu2- A„u2 + /u,,0). On ape iy(]0,T[J(^-1Jn))3), et < 5,v >(W-.(f,))3,(^(n))3= 0 Vu G (W1J(Q))3, divu = 0. Par le lemme précédent, il existe donc p € 2?'(]0, T[; L2(fl)) tel que g = Vp. Il reste à montrer que g e VF-1^(O1T; (//"'(fi))3). Par définition: W-l"°{0,T; H-1W)) = L0 + ^i. I 4 g L~{0,T;H-1W)] . Or u, € Z,~(0, T; L2(fì)) C ^(0,TiW(Ii)) donc u.-et^ieiv-^o.rjff-^n)). D'autre part ||A„u,-||„_1(n) = ||V„u,-||t,(n). Or u, <= L2(0,T;//'(O)) donc ||V„u,-||t2(n) e L2(0,T) C 1'(0,T). Donc A„u, e P(0,T; H^(U)). Lemme 6 Soit X un Banach, L1(0,T;X)CW-1-'X,(0,T;X) En effet, soit rp e Ll(0,T;X), posons ip(t) = /0'V>(s)ds. Alors tp e L°°(0, T; X) donc V = ^f € W"''00^ T; X). D Du lemme on déduit que A„u,- G W-1^(O, T; H~l(ü)). Considé- rons maintenant les termes a • Vu1 = div (u,a). En vertu du théorème 50 CHAPITRE 4. ÉCOULEMENT EN EAU PEU PROFONDE. d'injection de Sobolev u, € L2(0, T; H1W)) <-* L2(0, T; L6(U)). Pour la même raison aj e L2(0,T; L3W)) j = 1,2. Par l'hypothèse (4.16), a3 6 L2(0,T; L3W)). Donc par l'inégalité de Helder: u^eL^T;!2^)), j = l,2. Ainsi div(u.a) 6 Ll(0,T; H-1M) C W-l'°°{Q,T; H^[Si)), par le lemme précédent. Tous les termes constituant g sont donc dans W-1^(O1T; H-1W)). m Remarque. On ne peut se dispenser de l'hypothèse (4.16). Sans celle-ci on a seulement a3 e Z,2(0,T; L2(0)) et u< e L2(0,T; L6(fi)), donc par Helder um e L1 (0, T; L312W)) i I'(0, T; L2JO)). Q 4.9 Sens de la condition initiale et unicité. Définissons WH = ((U11U2) 6 H\ x W011 Bu3 £ /Z0(O3); u = (u„u2,u3) € w) Cet espace est en fait le projeté horizontal de W. Munissons le de la semi-norme produit sur Hl x Hl. Proposition 4 W^ est un espace de Hilbert. Preuve. Considérons (uj, UJ)716n une suite de Cauchy dans Wn- Grâce à la divergence nulle on a: O3 (u3 — U3") = —di (u" — u™) — O2 (uj — U2"). Donc (u3)neN est de Cauchy dans H0(O3). Ainsi il existe U3 6 Ho(O3) tel que uj —> U3. D'autre part divu" = 0, Vn. Donc en faisant n —> oo, on obtient pour u = (ui,u2,u3): divu = 0, Vn. ¦ Nous étudions maintenant la continuité de Uh = (ui,u2). Pour cela nous aurons besoin de la 4.9. SENS DE LA CONDITION INITIALE ET UNICITÉ. 51 Proposition 5 ^ €i'(O1T; H^) Preuve. On a Uh = (ut,u2) 6 L2{0,T;Wh), donc Pour t € (0,T), notons $(<) la forme linéaire vH = (vuv2) € WH •-> (u,,a- Vu,) + /• (u2,u,) - (V^u1, V„u,) +(u2, a- Vu2) - /• (u,, V2) - ( V„u,, V„i>2)+ < O1, V1 > + < 02, V7 > I <||«i||-l/»INI 1/2- < q|0.ll-./2INI«,;. Les autres termes ne posant pas de difficulté, on obtient IWOHn < cdiaïuiiuHii^ + ||uw|kH + !!¢,11-,/, + 11¾!!.,/,). Comme a e L2(0, T; L3(fl)) et uw e L2(0, T; W„) \\»Mv»\\wa € Ll(0,T). De plus \\uH\\wH € L2(0,T) C L1 (0,T) et ||ô.||_,/2 € L°°(0,T) C V(0,T). Ainsi ||$(<)l|ivi £ Ll(0,T), autrement dit $ £ Ll(0,T;W'H). Soit W// = (U)11W2) € Wh, l'application vh = (vi,v2) 6^h (Wh1Vh) = (tu,,mk» + (w2,v2)v e R 52 CHAPITRE 4. ÉCOULEMENT EN EAU PEU PROFONDE. est linéaire continue, donc par la théorie des distributions à valeurs vectorielles {*£tVH) = d(uuVi) + M € mn Ww e Wh En prenant des fonctions tests de la forme ip(t)v(x), V> 6 P(O1T), ve W dans (4.19) on obtient la forme faible spatiale suivante —--------^----------= (ui,a • Vu1) + / • (u2lui) - (VkMi1V1^1) +(u2l a -Vu2)-/-(U11U2)-(VkU21V1^2) dans î>'(0, T). (4.30) Avec la définition de $ (^p,vH) =<*(0,VH >wi,Hr. dans 15'(O1T)1 Vv^eW», donc (%Sv„) 6 Z-1 (O1T),et ^€£'(0,r;WH). ¦ Comme uw € L2(0,T;H'//), la proposition précédente donne, après modification éventuelle sur une partie négligeable de (O1T)1 uw € C(O1 T; W11). Proposition 6 Pour tout (ui,u2) € Wh la fonction *-> (U1(O1U1) + (U2(O1U2) est absolument continue sur [O1T]. Preuve. En effet par la proposition précédente <««„*.)+ (*„«*) gL, (or) 4.9. SENS DE LA CONDITION INITIALE ET UNICITÉ. 53 On peut dès lors préciser le sens de la condition initiale. La forme faible spatiale (4.30) s'écrit: ^pI - (Ul,a ¦ Vu1) - / ¦ (U21U1) + (V11U11V1,!,,)- < B11V1 > + j, - (u2, a • Vu2) + / • (Ui1U2) + (V„u2, V„Vî)- < 02lu2 > al = 0 dans V{0, T) Vv=(U11U21U3)CW Prenant tp € Cl[0,T], vérifiant 0(0) = 1 et xj>(T) = 0. Par absolue continuité de (U1(^)1Ui) + (u2(t),u2) [T U(UuV1) d(Ui,V2)\ Jo \ dt dt ) v = - / dt((U1,U1) + (U21U2)) t/>'(<) +MO),«) + (U2(O)1U) Jo = / dt {(u,,a-Vt>,)-/- (u2,vi) + (V^u1, V^u1)- < 6uvt >} V JO [T + / dt {(u2la- Vu2) + /(U11U2) + (V^u21V17U2)- <02lu2 >} V Jo En comparant avec la forme faible spatio-temporelle (4.19), il vient (U1(O)1U1) + (U2(O)1U2) = 0, Vv = (vuv2) € WH Donc U1(O) = u2(0) = 0 dans W'„. Remarque. Pour la vitesse verticale u3l on ne peut rien dire en général de la continuité, ni a fortiori de la condition initiale. O 4.9.1 Unicité de la solution. On suppose de plus pour simplifier ae r°()o,:r[xn). Cette régularité supplémentaire permet d'établir comme précédemment que """ eL\Q,T; W'„). (4.31) dt On a alors le 54 CHAPITRE 4. ÉCOULEMENT EN EAU PEU PROFONDE. Théorème 11 Le problème hydrostatique linéarisé (4-19) admet une solution unique u 6 £2(0,T;W) Preuve. Avec (4.31), u„ e L2(0, T; WH) et ^ g L2(0, T; ^). Par le lemme 1.2 p 261 de Temam [58], %Ì = 2 < u', u > et u„ 6 C(O, T; L2(H)). ai Donc la fonction < ^ IMOIlS+ IMOIlS est continue sur [0,T], après modification éventuelle sur une partie né- gligeable. Supposons l'existence de deux solutions u et v, et montrons qu'elle coincident. Posons w = u — v. En prenant w = (wi, 102,1113) dans la forme faible spatiale (4.30), en tenant compte des simplifications pour raison de symétrie et de divergence nulle, on obtient ^+^+(V.»,, V^)+(V^, V^2) = 0 .UnSD-(O1T). Donc %l + «o0) U(I1O) =p(x), *€]0,1[ 59 60CHAPITRE 5. MOINDRES CARRÉS DANS L'ESPACE-TEMPS. u(0,t) =h{t), <6]0,l[ multiplions par Elevons au carré l'opérateur différentiel vdx-rdt pour obtenir l'équation l'adjoint -{udx + dt) de diffusion. U2U1x + 2vuxt + utt = 0, (i, t) 6]0, l[x]0, l[ (5.2) avec la condition de Dirichlet U(S1O) =g(x), xe]0,l[ u(0,0 =h(t), t6]0,l[ Le problème est que pour une équation parabolique, il faut une condi- tion limite sur tout le bord du domaine, alors que pour une équation hyperbolique, il faut seulement une condition limite sur le bord à flux entrant. Le fait étonnant est qu'avec des conditions limites adéquates, les deux problèmes sont équivalents. De plus, comme nous allons voir, la mé- thode de Galerkin standard appliquée à l'équation de diffusion fournit un bon algorithme pour l'équation d'advection. 5.2 Principe des méthodes aux moindres carrés pour les EDP. Des applications variées de la méthode des moindres carrés aux EDP sont présentées par exemple dans [24]. La méthode des moindres carrés permet en effet de "convexifier" des problèmes non coercitifs. L'idée de base est la même qu'en dimension finie: on transforme un système non dégénéré A ¦ x = b forme normale. en un système défini positif ATA ¦ x = ATb. Récemment la méthode des moindres carrés a permis d'améliorer la fameuse méthode SUPG (Streamline Upwind Petrov Galerkin) [15] en la méthode GLS Galerkin Least Squares [28]. Dans un même esprit les 5.2. PRINCIPE DES MÉTHODES AUX MOINDRES CARRÉS. 61 célèbres méthodes Streamline Diffusion et Discontinuous Galerkin [29, 30] ont muté pour donner la méthode CSD Characteristic Streamline Diffusion [25]. Le principe commun à ces méthodes est de stabiliser le problème en ajoutant une petite quantité de résidu quadratique des termes advectifs. Considérons par exemple l'équation-type Au = a • Vu = / dans fi (5.3) u =0 sur T- (5.4) où T- = {x 6 r| a • n < 0} désigne la frontière à flux entrant. Une formulation faible de ce problème aux limites est: B(u,4>) = L(),VeV (5.5) où V = {^eHl\=0onT-}, LW = / }, Jn B(u,)= /(a-Vu)Oi. Jn Soit Vh un sous espace de dimension finie de V (par exemple obtenu par une discrétisation en éléments finis, la méthode de Galerkin revient à: Trouver Uj1 € V/, tel que B(uh,) = L(), V G V/,. Malheureusement, il est bien connu le schéma résultant est instable et donne lieu à des oscillations parasites (voir par exemple [29]). Les méthodes (GLS, CSD) présentant une meilleure stabilité sont défi- nies par les formes bilinéaires et linéaires Bs(uh,t) = Ls(), Bs(uk, ) = B(uh, ) + r(Auh, A), les instabilités clas- siques des schémas centrés pour les pro- blèmes advectifs. Ls() = L(t) + T(f,A). (,2CHAPlTRE 5. MOINDRES CARRÉS DANS L'ESPACE-TEMPS. Le désavantage est qu'il faut maintenant choisir le paramètre r. De plus ce n'est pas la méthode des moindres carrés stricto sensu, définie elle par les formes LLS(t) = (LM)- Nous allons voir que cette méthode est plus simple et aussi efficace que les méthodes GLS, CSD. Nguyen and Reynen [40] sont les pre- miers à notre connaissance à avoir appliqué simplement la méthode des moindres carrés dans l'espace-temps, mais leur article ne semble pas avoir trouvé l'écho qu'il méritait. 5.3 Le cadre spatio-temporel. L'analyse précédente concerne un problème stationnaire. Nous traitons maintenant le cas instationnaire. Considérons l'équation-type: Au = u, + a-Vu =/ dans Ox]O, T[ (5.6) u =g surrx]0,T[ (5.7) u(.,t=0) = h sur fi (5.8) Comme il n'y a pas de raison de traiter la variable temps différemment des variables spatiales, introduisons les objets spatio-temporels suivants (voir fig. 5.1): Le domaine Q = fix]0,T[ Le bord dQ = diïx}0,T[ + Q x({T}-{0}) ( (n,0) sur dfix]0,T[, La normale sortante n = < (0,1) sur fi x {T}, { (0,-1) sur fix {0} La vitesse advective S = (a, 1) La frontière-influx G~ = {{x, t) 6 dQ \ à • n < 0} = T" x]0, T[+fi x {0} Le gradient V = (V, dt). Maintenant le problème est "stationnaire": Au = à-Vu = f dans Q (5.9) u =p sur G~ (5.10) 5.4. LA METHODE STILS. 63 Distance [m] Figure 5.1: domaine, frontière influx, normale sortante spatio- temporels. Remarque. Le cadre spatio-temporel unifie condition sur le bord entrant et condition initiale. O L'idée est maintenant de conjuguer méthode aux moindres carrés et formulation spatio-temporelle pour donner la méthode STILS (Space Time Integrated Least Squares method) [44, 43, 6]. 5.4 la méthode STILS (moindres carrés intégrés dans l'espace-temps). 5.4.1 Forme variationnelle aux moindres carrés. A partir de maintenant, nous nous plaçons dans l'espace-temps. Reé- crivons l'équation (5.9) sous la forme équivalente: ir(u) = l/2 / [Au-J)2 Jq (5.11) où les u admissible doivent remplir la condition de Dirichlet sur le bord éclairé (5.10). Le minimum zéro correspond à la solution de (5.9-5.10). La condition de stationnarité se traduit par la forme faible: Euler- Lagrange 64CHAPITRE 5. MOINDRES CARRÉS DANS L'ESPACE-TEMPS. Sn = / (â • Vu - /)) = /" /(a • V<£) (5.21) Ja Ja De la théorie des équations linéaires hyperboliques [18], nous savons qu'il y a une et une seule solution au problème advectif (5.19). Si cette solution est dans L2, elle sera aussi solution du problème diffusif (5.21), comme on le voit aisément en effectuant le produit scalaire L2 de (5.19) par a • V^ pour un^e H0(a, fi, T-) arbitraire. Nous donnons ci-dessous une condition suffisante pour que le problème de diffusion (5.21) soit bien posé, à savoir l'inégalité de Poincaré "cour- be" que nous étudions au chapitre suivant: 3C tel que V<£ € %(a, O, T" ), \\\\0 < C||a • V<£||0 (5.22) Théorème 12 Si f £ L2(Q) et si l'inégalité de Poincaré "courbe" (5.22) est satisfaite, le problème (5.21) est bien posé. Preuve. Sur V = //o(a, fi, T-), définissons la forme bilinéaire B(^1V)= /(a -V*)(a -VV-), Jn et la forme linéaire LW = / /(a • VtA). Jq Elles sont naturellement continues sur V, et la forme bilinéaire est coer- citive grâce à l'inégalité (5.22): B(«,«)= ||a-V«||S > CIIHII2. Par le lemme de Lax-Milgram, le problème (5.21) admet une solution unique et |||«||| ) 5.4. LA MÉTHODE STILS. 67 par le lemme de Céa et la théorie standard de l'interpolation [16], où p est le degré des polynômes utilisés, en particulier pour des éléments spatio-temporels Q1 cela donne l'estimation d'erreur |||u—u/,||| = 0(h), du même ordre que celle obtenue dans [28]. Le chapitre suivant est consacré à l'inégalité de Poincaré "courbe". Nous prouvons en particulier le Théorème 13 Si la vitesse advective a est bornée et de classe C1 au voisinage de Q, si a • k > o > 0 pour un vecteur unitaire fixe k et si Çï est borné dans la direction k alors l'inégalité de Poincaré "courbe" est vraie. L'important est que Remarque. Si le problème est instationnaire et si la vitesse advective a es "9nes de cou- „. ,, . ,. „ _ , , , .„ ... _ rant du flot associé est C, 1 inégalité est vraie. 11 surfit de prendre k = (0,1) ! D . , . . Q ai x, 11 tic sotenz pas trop "sauvages", que le flot soit "Ia- 5.4.3 Retour à l'exemple. minaire" en quelque Considérons à nouveau l'équation d'advection scalaire transiente (5.1) avec une vitesse constante v et l'équation de diffusion associée (5.2). a = Kl),[a®a] = (^ *) La frontière à flux sortant est {x = IJx]O1T = l[+]0,i = l[x{T = 1), La normale sortante est, suivant la portion du bord considérée n = (l,0)ou (0,1). La conormale est (i^2,i/) colinéaire à (i/, 1) ainsi la condition limite na- turelle (5.18) devient VUx + ut=0 Dans ce cas, les hypothèses du théorème 13 sont vérifiées, validant ainsi la mise en oeuvre de la méthode STILS. sorte. 68CHAPiTRE 5. MOINDRES CARRÉS DANS L'ESPACE-TEMPS. Distance [m] Figure 5.2: Tranche spatio-temporelle. 5.5 Comportement numérique. Nous restreignons notre analyse à l'exemple ID (i.e. 2D en espace- temps). Pour la discrétisation, nous prenons des éléments Qj. Pour pouvoir changer le pas de temps sans devoir remailler le domaine, nous posons ( = flt et résolvons le problème dans le plan (x, Q avec Ai = AÇ. Pour mettre les équations sous forme adimensionnelle , nous utilisons le nombre de Courant Cr = i/.At/Ax = vjß. 5.5.1 La discrétisation de Galerkin. En dimension 2 or 3 d'espace, le système spatio-temporel est onéreux à résoudre, car le nombre d'inconnues est élevé. En particulier en dimen- sion spatiale 3, il faut mailler un domaine spatio-temporel de dimension 4. Nous évitons ces difficultés en discrétisant par tranches successives (voir fig. 5.2) flx](,i + At[, pour aboutir à un système en (u',ut+At). La matrice de rigidité par élément est: (Cr2/3 +Cr/2+1/3 -Cr2/3+l/G -Cr7/6 - Cr/2 - 1/6 Cr2/6 - 1/3 Cr2 /3 -Cr/2+ 1/3 Cr2/6-I/3 -Cr2/6+Cr/2 - 1/6 • . Cr2/3 + Cr/2+ 1/3 -Cr2/3+l/6 Cr2/3- Cr/2+ 1/3 5.5. COMPORTEMENT NUMÉRIQUE. 69 La jc équation est typiquement: (1 - 2Cr2)u$ +(1 - 3Cr + Cr2JuJ+1 Remarque. C'est un schéma implicite, mais il nécessite seulement l'in- version d'une matrice tridiagonale. Avec un développement de Taylor, on obtient aisément que c'est un schéma du second-ordre. O l'erreur de dis- _ ... , , , ., crétisation est du Cette equation a déjà ete obtenue dans [40]. troisième degré en Ai, At 5.5.2 L'analyse de Fourier-von Neumann et la ma- trice d'amplification. La réponse de l'équation discrète à une harmonique simple e'*1 est donnée par l'amplification complexe: _ 2 - Cr2 + ( 1 + Cr2) cos6 - i ¦ 3Cr sin 0 7_ 2 + 2Cr2+ (1-2Cr2) cosò où 8 appelé angle de phase vaut fc ¦ Ai. Le schéma est inconditionnellement stable: \f\ < 1, VB, Cr. Nous avons calculé l'amortissement et le déphasage après une oscillation complète dans la table (5.1), où le retard de phase (en unité de longueur d'onde) est A=I- arccos(»(7/|7|))/0Cr, l'amortissement après une oscillation est g=|7r/9cr. Ceci montre que STILS se comporte convenablement au moins dans la région des discrétisations spatiales réalistes 0 < n/4, et se comporte aussi bien que les autres méthodes pour Cr < 1 (bien qu'étant un peu moins précise) [45]. Elle se révèle meilleure pour Cr > 1 (simplement parce que les autres méthodes explosent alors!). Mais l'analyse de Fou- rier n'est qu'un aspect des choses. L'étude de la matrice d'amplification donne un éclairage complémen- taire. L'algorithme d'avance en temps conduit au système linéaire Su'+* = Ru', IOCHAPITRE 5. MOINDRES CARRÉS DANS L'ESPACE-TEMPS. Cr e AG 0.2 tt/4 0.0061 0.9717 tt/2 0.0544 0.7462 3tt/4 0.2923 0.2707 0.5 tt/4 0.0251 0.9237 tt/2 0.0977 0.5220 3tt/4 0.2479 0.0910 0.9 tt/4 0.0690 0.8447 tt/2 0.1825 0.4030 3tt/4 0.2371 0.1074 1.2 tt/4 0.1096 0.7803 tt/2 0.2485 0.3775 3tt/4 0.2926 0.1703 1.6 tt/4 0.1663 0.7009 tt/2 0.3288 0.3795 3tt/4 0.3880 0.2708 2.0 tt/4 0.2211 0.6371 tt/2 0.3976 0.4000 3tt/4 0.4721 0.3618 Tableau 5.1: Amortissement et déphasage. 5.6. EXEMPLES DE SIMULATIONS. 71 Figure 5.3: Norme de la matrice d'amplification. posant G = S'1 R nous obtenons: u,+A< = Gu*. G est la matrice d'amplification du schéma. Son spectre contrôle l'am- plification des erreurs de troncature au cours des itérations. P. Per- rochet [44] a calculé la norme spectrale ||G|| en fonction de Cr pour les schémas standards CNG (Crank-Nicolson), SUPG, CNTG (Crank- Nicolson-Taylor-Galerkin) [45], avec des maillages de plus de 20 élé- ments, que nous avons reporté dans la figure (5.3). Notre schéma montre de bonnes propriétés. 5.6 Exemples de simulations numériques. La simplicité conceptuelle de STILS se reflète aussi dans sa program- mation, le code étant similaire à celui d'un laplacien, avec pour seule différence la procédure d'assemblage de la matrice de rigidité, dont le terme général est (a-V^>,,a-V^) au lieu de (Vj, S7j). Nous donnons en appendice, un exemple de code. 12CHAPlTRE 5. MOINDRES CARRÉS DANS L'ESPACE-TEMPS. Nous référons à [44] pour des exemples issus de problèmes hydrogéo- logiques, où STILS se révèle particulièrement, performant, notamment dans les cas délicats d'écoulements en milieu karstique. Nous reprodui- sons ci-après les simulations conduites par Pierre Perrochet, hydrogéo- logue à l'EPF Lausanne. Injection dans un milieu hétérogène fissuré. Considérons un domaine vertical composé de deux couches à haute perméabilité (K = 10-3m/s) connectées par une fissure à travers une couche quasi-imperméable (K = 10_8m/s), voir fig. 5.4. Ce domaine en (x,z) est discrétisé en 20x20 éléments biquadratiques (Ax = Az — Im) et la fissure au moyen d'éléments quadratiques ID à haute trans- missivité (A'e = 10~2m2/s où e désigne l'ouverture de la fissure. La figure 5.4 montre le champ de vitesse a, laplacien stationnaire, engen- dré par les conditions de flux nul sur les faces latérales et de charge hydraulique imposée aux coins inférieurs ($ = 100m) et supérieurs (ty = 50m). On fixe une concentration u = 1 (resp. u = 0) de pol- luant au coin inférieur gauche (resp. droit). La figure 5.5 illustre le domaine 3D spatio-temporel. La figure 5.6 représente les solutions après 10000s et 20000s obtenues après seulement 5 et 10 pas de temps (AC = 2m, Ai = 20003,/3 = 10~3m/.s). Les tentatives pour simuler ce problème avec le schéma d'ordre plus élevé CNTG (Crank-Nicolson-Taylor-Galerkin) ont échoué, malgré des pas de temps jusqu'à 500 fois plus faibles. En comparaison les oscil- lations numériques parasites visibles sur la fig. 5.6 sont stables et ne dépassent pas 7-8%, ce qui est tout à fait acceptable vu la difficulté du problème. Ces résultats réalistes permettent une interprétation rai- sonnable. Comme des quantités égales de fluide entrent par les deux coins inférieurs, on constate bien qu'après environ 20000s le mélange est atteint: passé cette durée la fissure agit comme une source ponctuelle avec u = 0.5, et le front dans la couche supérieure est beaucoup plus régulier à cause de l'effet de mélange. 5.6. EXEMPLES DE SIMULATIONS. 73 Figure 5.4: (a) Aquifères séparés par une roche fissurée, (b) Equipo- tentielles ty et lignes de courant de l'écoulement. r+ r+ k .^..Jlïtcy i i-1 *~r-: a(20,0,0=0 ^r-: u(x, z, O)=O Figure 5.5: Représentation 3D du milieu et des caractéristiques. 74CHAPJTRE 5. MOINDRES CARRÉS DANS L'ESPACE-TEMPS. \ / X -¦-¦--- -x K=I U=O It=I «=0 Figure 5.6: Solutions après (a) 5 pas de temps = 10 000s, (b) 10 pas de temps = 20 000s. Chapitre 6 Inégalité de Poincaré courbe pour le traitement variationnel de l'équation de transport. Il est possible de donner à l'équation de transport un traitement va- riationnel analogue à celui des équations elliptiques. Pour obtenir une forme bilinéaire coercitive, nous introduisons un espace de Sobolev ap- proprié, des fonctions tests spécialement adaptées au transport et nous travaillons dans l'espace-temps. La méthode repose sur une inégalité de Poincaré "courbe", dont nous étudions la validité. 6.1 Cadre fonctionnel. 6.1.1 Espaces de Sobolev obliques. Soit il un domaine de M" de frontière C1 par morceaux. Soit V un voisinage ouvert de fi. Soit a : V —? E" un champ de vecteurs C1. On définit Jï(a,fi) comme le complété de T)(Cl) muni de la norme UNII2= HS+ 1Ia-VuIIS, où ||.||o désigne la norme I?. Le vecteur n désignant la normale sortante est défini sauf peut-être sur une partie négligeable Af, on note T- = {i 6 dil \ M\ a • n < 0}, 75 76 CHAPITRE 6. INÉGALITÉ DE POINCARÉ COURBE (resp. T+) la frontière à flux entrant (resp. sortant), supposées de mesure strictement positive, et T0 la frontière à flux nul. On définit V(U, r-) = { G V(U) I i> = 0 sur T-) . On désigne par H0(a, 0, T~) (resp. %(a, H1T+)) la fermeture de V(U, T~) (resp. V(Ù,T+)) dans tf(a,fi). Remarque. On peut définir ces espaces différemment, comme au cha- pitre précédent, en utilisant un opérateur de trace comme dans [9, 22, 18] et les définitions en sont équivalentes pour un domaine de frontière C1 par morceaux. O 6.1.2 Écoulement remplissant. Soit l'écoulement ou flot associé au champ a, décrit par les trajectoires ou courbes intégrales maximales £ : (.s,x) G [U1)T1] xfti-t £(s,x) G 0 solutions de g = a(0 (6.1) C(O, x) = x (6.2) Définition 1 Un écoulement est appelé iï-remplissant si il existe T > 0, une partie négligeable V pour la mesure de Lebesgue sur f), tels que pour tout x de n \ V il existe x0 € T- et 0 < t < T tels que x = Ç(i,x0). Remarque. Cette définition peut se reformuler ainsi: les trajectoires issues du bord à flux entrant remplissent Q, sauf peut-être une partie négligeable, en une durée finie, bornée par un nombre fixe T. ? Une condition suffisante de régularité est donnée par la Proposition 7 si le champ a est C1 et borné sur un voisinage V de fi, s'il existe une direction donnée par un vecteur unitaire fixe k, un nombre a > 0 tels que Vx G A, a(x) • k > a (6.3) 6.1. CADRE FONCTIONNEL. 77 et que le domaine soit borné dans cette direction alors l'écoulement est Çl-remplissant. Preuve. Par le théorème de Cauchy-Lipschitz appliqué dans V, tout point x de Û est situé sur une courbe intégrale maximale. Le domaine est borné dans la direction k. Notons diamit(Q) = sup {(x — y) ¦ k}. !,yen On a ak > a > 0, donc il existe Cx < 0 tel que Xo = f (¾, x) G dû, avec \ax\ < diami,(n)/a. Prouvons que le complémentaire du sous-ensemble {ï€n|ç( 0, a(x0) ¦ n < 0. Donc X0 € T" U T0. Enfin par le lemme de Sard, on montre comme dans [9] prop. 2.3 p.194, que {x e AUK,*) e r°} est négligeable. ¦ Remarque. Le domaine n'étant pas a priori borné, l'hypothèse que le champ a soit borné est nécessaire pour éviter que les courbes intégrales ne "s'échappent à l'infini" sans pouvoir sortir de O. O Remarque. La régularité du champ sert à assurer l'existence de courbes intégrales maximales. En fait il suffirait que a soit lipschitzien, par exemple. On pourrait même se contenter d'une régularité de type So- bolev et de solutions renormalisées (voir [20]). D 78 CHAPITRE 6. INÉGALITÉ DE POINCARÉ COURBE 6.1.3 Inégalité de Poincaré courbe. Les définitions précédentes nous permettent d'énoncer le Théorème 14 Pour un écoulement iï-remplissant, gouverné par un champ à divergence nulle, il existe une constante C telle que Vu6ffo(a,n,r-), HIo < C||a • Vti||o. (6.4) L'idée naturelle de preuve consiste à redresser globalement le champ de vecteurs gouvernant le flot et d'appliquer ensuite l'inégalité classique de Poincaré dans une direction donnée. Cette idée simple peut être réalisée [56] et nous donnons cette démonstration plus loin, mais nous présentons d'abord une démonstration "hilbertienne" qui se généralise pour le cas de champs a moins réguliers. Preuve. Soit x € T-, alors (T1 = 0 et Tx représente la durée de la trajectoire issue de x. La régularité de l'écoulement entraîne Tx (Ù, T"). ||a-VuIl0= sup {Jjl~r,------}• *6L2\{0} ||0||o Soit p fourni par le lemme 7, quitte à régulariser on peut supposer que p 6 V[U, T+). Prenons alors ~ pu. On a (pu)a- Vu = l/2(pa-Vu2). Puis, en intégrant par parties sur Si, les termes de bord disparaissent à cause de la nullité de p sur T+, de u sur T-, de a • n sur T0. Comme diva = 0, le choix de p conduit à 1/2 J pA- Vu2dx= I u2dx= ||u||jj D'autre part on sait que la norme ||pu||0 est majorée par 1IpII00IIuIIo, on en déduit IpW*- Wo -I|a-Vu|k 6.1. CADRE FONCTIONNEL. 79 ©1 Figure 6.1: Le flot prolongé. Ainsi (6.4) est démontrée avec C = ||/j||oo- ¦ Donnons maintenant une démonstration plus géométrique du théorème. Pour simplifier l'exposé, on prend comme premier vecteur de base Étape 1: Prolongement du champ à ffin. On utilise pour cela une partition de l'unité subordonnée au voisinage V de H et fic définie par une fonction A : x i-> A(x), telle que A = I sur fi et A = 0 hors de V. On pose alors w(x) = A(x)a(x) + (l-A(x))e,. Le champ w est défini sur R", et il prolonge a. En particulier ses courbes intégrales prolongent celles de a hors de fi, ne modifiant en aucune façon les frontières à flux entrant, sortant, nul. Enfin l'hypothèse (6.3) est préservée: w • ei = Aa • ei + (1 — A) est compris entre 1 et a Etape 2: Construction du difféomorphisme redressant. Comme fi est borné dans la direction ei, il existe un hyperplan normal (il = c) tel que: (xi = c)nfi = 0. 80 CHAPITRE 6. INÉGALITÉ DE POINCARÉ COURBE Paramétrons les caractéristiques par leur impact sur cet hyperplan en désignant par s>->i(s,b) la courbe intégrale maximale passant par x = (c, b). Elle est définie sur K par maximalité, de plus lim f, = ±00, car ds ¦ Wi > a > 0. Avec le théorème de différentiabilité du flot [2], on en déduit que: f : R X R"-1 -> R" (s,b) n Ç(s,b) est un C'-difféomorphisme global, redressant le champ w: autrement dit, dériver par rapport à s dans les coordonnées (s, b) re- vient à dériver dans la direction w. Etape 3: Changement de variable pour se ramener au cas rectiligne. Grâce au difféomorphisme £, les lignes de champ de w sont les lignes de coordonnées b = Cte, s € R. Notant û = u o Ç, les nouvelles coor- données z = (s, b) 6 n t-n = £(s, b) € il, cela signifie que: dû dix (x) = a • Vu(i). L' inégalité à prouver s'exprime alors par: h2 D(i D(x) dx 1/2 et valant l'unité pour r < 1/4 par exemple. Manifestement a • Vu = 0, violant ainsi (6.4). Voir fig. 6.2. On a la propriété étonnante: Proposition 8 Un écoulement instationnaire gouverné par un champ C à divergence nulle dans un voisinage V de Ù est Q-remplissant. 82 CHAPITRE 6. INÉGALITÉ DE POINCARÉ COURBE Preuve. Il suffit de prendre k = (0,1) . ¦ Nous pouvons maintenant énoncer, sous forme légèrement plus géné- rale, le théorème 12 démontré au chapitre précédent: Théorème 15 St l'écoulement est fl-remplissant, alors le problème (5.21) est bien posé. 6.2 Équivalence entre les formulations ad- vective et diffusive. Soit / 6 L2(fi), pour l'équation a • Vu = / dans 0 la ., u =0 surr" (6'6) nous examinons maintenant l'équivalence entre les formulations faibles suivantes: Trouver u € H0{a,fi, T") tel que \/ 6 #o(a,fi, T-), /(a Vu)(a-V^)= //(a -W) (6.7) Jn Jn et Trouver u e H0(a, fi, T") tel que V<£ e L2(ïl), /(a-Vti)*=//* (6.8) Remarque. Dans le chapitre précédent, nous avons vu que la méthode STILS correspondant à (6.7) revient à chercher le minimum parmi les fonctions admissibles de tt(u) = 1/2 / {Au - /)2 ->¦ min. Ja En revanche, la formulation faible de Galerkin standard (6.8), n'étant pas symétrique, ne correspond pas à un principe de minimum. O Tout d'abord toute solution de (6.8) est aussi solution de (6.7). De plus, si le problème advectif (6.8) est bien posé dans L2(Q), cela signifie que NIo e L2(fi) orthogonal à {a • V, €D(A,T-)} , montrons que 0 = 0. L'orthogonalité se traduit par: / (a-V(M = O, V<£€ D(A1T-). (6-9) Jf! Prenons en particulier é> € V(Q.) et intégrons par parties / d\v((f>-a)ip = 0. Ja 84 CHAPITRE 6. INÉGALITÉ DE POINCARÉ COURBE Ceci signifie que < VV>, 4> ¦ a >= O dans P(U). Comme Vtp est une distribution d'ordre 1, que a est de classe C1, la distribution a • Vip est bien définie et vérifie < a- Vip,d>>= 0. Donc la distribution a • V^ = 06 L2(û), ainsi, par la définition donnée au chapitre précédent, rl>eH(a,(l). On peut donc [18, 22] considérer sa trace sur le bord T+, le bord étant supposé C1 par morceaux. Prenant maintenant € 25(Q1T-) non nul sur T+, intégrant à nouveau (6.9) par parties conduit à / èé a • n ds = 0. Jr* Comme ip est arbitraire sur T+, que a • n > 0 sur F+, cela entraîne V'ir+ = 0. Appliquant enfin l'inégalité de Poincaré courbe pour le champ —a issu inégalité de Poincaré de T+, il résulte courbe pour é = 0 l'écoulement dual. m Remarque. La condition (6.3) doit être vérifiée jusqu'au bord du do- maine, sinon il peut ne pas y avoir de solution dans L2(f!): prenons par exemple a = (y,0), alors l'équation ydxu = 1 dans le carré ]0,1[2 avec u = 0 sur T- = {i = 0} admet une unique solution u = x/y gui n'est pas dans L2(ü). Voir fig. 6.3. O 6.2. ÉQUIVALENCE ENTRE LES FORMULATIONS. 85 Figure 6.3: Ecoulement de Couette. 86 CHAPITRE 6. INÉGALITÉ DE POINCARÉ COURBE Annexe A Nous utilisons dans les programmes suivants les routines de la biblio- thèque FEK (Finite Element Kernel) développée à la Hebrew University of Jerusalem par l'équipe de M. Bercovier. A.l Élément fini hydrostatique. La brique QO.5. subroutine b4shap (xl,x2,x3,r,s,t,shp,dl,dg,det,iopt,eps.ierr) ccbegin c c-------------------------------------------------------------------- c call Mshap (xl.x2,x3,r,s,t.shp.dl.dg.det,iopt,eps.ierr) c c c usage c ==s== c c quadrilateral element, in three dimensions. c linear shape functions in xl and x2, constant in x3. c c purpose c ------- C c finds shape function local and global derivatives and jacobian c of transformation, at a point specified by the coordinates (r.s.t) c . in the reference element (t = 0.0). c ccdoc c input arguments c --------------- C c xl = vector of length 8 containing first coordinates c of element nodal points. c x2 = vector of length 8 containing second coordinates 87 88 ANNEXE A. of element nodal points. x3 = vector of length 8 containing third coordinates of element nodal points. r = user specified first coordinate. s = user specified second coordinate. t = user specified third coordinate (= 0.0). mndp = row dimension of dadg and dadi in calling program. iopt = option paramater (see note 1). eps = criterion for checking zero determinant. output arguments = vector of length 4 containing shape functions. = matrix of order 4x3 containing local derivatives. = matrix of order 4x3 containing global derivatives. = jacobian of transformation. = error flag see note 2. 1. iopt = 0 shape functions only. 1 local derivatives. 2 local derivatives + jacobian. 3 local derivatives + jacobian + global derivatives. -1,-2 and -3 outputs the shape functions as well. c shp C dl e dg C det C ierr C C notes C C ierr = 0 successful execution 1 zero jacobian. -1 negative jacobian. cc end implicit none real»8 xl,x2,x3,r,s,t,shp,dl,dg,det,eps integer iopt f ierr dimension xl(8),x2(8),x3(8),shp(4),dl(4,3),dg(4,3) real*8 shape,dadi,a,b,zero,un,us2 real*8 rp,rm,sp,sm,tp,tm,cl,c2,c3,c integer ielt,i, j ,k,ioptl dimension shape(8) ,dadl(8,3) dimension a(3,3),b(3,3) data zero,us2,un / 0.OdO1O.BdO,i.dO/ ierr = 0 ielt = 8 ÉLÉMENT FINI HYDROSTATIQUE. rp = us2*(un + r) rm = us2*(un - r) sp = us2*(un + s) sm = us2*(un - s) tp = us2»(un + zero) tm = us2*(un - zero) if (iopt .gt. 0) go to 10 shape functions shp(l) = rm»sm shp(2) = rp*sœ shp(3) = rp*sp shp(4) = rm»sp shaped ) = nn*snt*tm shape(2) = rp*sm»tm shape(3) = rp*sp*tm shape(4) = rm*sp»tm shape(B) = nn*sm»tp shape(6) = rp»sm*tp shape(7) = rp*sp*tp shape(8) = rm*sp*tp il (iopt .eq. 0) return local derivatives ioptl = iabs (iopt) dl(l,l) = -us2*sm dl(2,l) = us2*sm dl(3,l) = us2*sp dl(4,l) = -us2*sp dl(l,2) = -us2*rn dl(2,2) = -us2*rp dl(3,2) = us2«rp dl(4,2) = us2*rm dl(i,3) = zero dl(2,3) = zero dl(3,3) = zero dl(4,3) = zero dadl(l,l) = -us2*sm*tm dadl(2,l) = us2*sm*tm dadl(3,l) = us2*sp*tm dadl(4,l) = -us2*sp*tm dadi(6,1) = -us2*sm*tp dadi(6,1) = us2»sm»tp dadl(7,l) = us2*sp*tp dadl(8,l) = -us2*sp»tp 90 dadi(1,2) = -us2*rm«tm dadl(2,2) = -us2*rp*tm dadl(3,2) = us2*rp»to dadl(4,2) = us2*rm*tm dadl(5,2) = -us2*nn*tp dadl(6,2) = -us2*rp«tp dadl(7,2) = us2»rp«tp dadl(8,2) = us2*rm*tp e dadl(l,3) = -us2*rm*sm dadi(2,3) = -us2*rp«sm dadl(3,3) = -us2*rp*sp dadl(4,3) = -us2*nn«sp dadl(6,3) = -dadl(l,3) dadl(6,3) = -dadl(2,3) . dadl(7,3) = -dadl(3,3) dadl(8,3) = -dadl(4.3) if (ioptl .eq. 1) return e e jacobian matrix a e do i = 1,3 ci c2 c3 do j ci c2 c3 end d a(i,l a(i,2 a(i,3 end do e e invert jacobian e do i = 1,3 j = i + 1 il (j .eq. 4) j = 1 k = j + 1 il (k .eq. 4) k = 1 b(i,i)= a(j,j)*a(k,k) b(i,j)= a(k,j)*a(i,k) b(j,i)= a(j,k)»a(k,i) end do e e find determinant of jacobian matrix. e det = a(l,l)»b(l,l) + a(l,2)*b(2,l) + a(l,3)»b(3,l) = zero = zero = zero = 1.8 = Cl = c2 = c3 O ) = Cl ) = c2 ) = c3 dadl(j,i)*xl(j) dadl(j,i)*x2(j) dadl(j,i)*x3(j) - a(k,j)*a(j,k) - a(i,j)*a(k,k) - a(j,i)*a(k,k) A.2. LE CODE STILS. 91 C c check determinant of jacobian e if (abs(det) .lt. eps) go to 70 if (det .gt. zero) go to 45 ierr =-1 call erfac (4hb4sp,ierr) 45 if (ioptl .eq. 2) return c c find global derivatives c do i = 1,3 do j = 1,4 c = b(i,l)*dl(j,l) + b(i,2)*dl(j,2) 1 + b(i,3)*dl(j,3) dg(j ,i) = c/det end do end do return c c zero jacobian c 70 continue ierr = 1 call erfac (4hb4sp,ierr) end A.2 Le code STILS. PROGRAN STYLE C C=============================================================== C PROBLEM C ------- C C THE PROBLEM IS TO SOLVE THE PARTIAL DIFFERENTIAL EqUATIOH : C U_{t} + aU_{x} = 0 REPLACIHO IT BY C U_{tt} + 2aU_{tx} + a"{2}U_{xx} = 0 on (0,Dx(O1I) C ATTEHTIOH Xl = x, X2 = t C**************************************************************** C THE IHITIAL-BOUHDARY CONDITIOHS ARE: C U(x,0) = f(x); U(O,t) = g(t) C c=============================================================== C IMPLICIT REAL*8 (A-H,0-Z) DIHEHSIOH A(1600),B(20O),MDIAG(2OO),B0UND(200),Y(2OO) DIMEHSIOH X(2,200),IEH(4,200),IA(200) DIMEHSIOH STIFF(4,4).BE(4),LH(4),XX(4),YY(4) 92 ANNEXE A. e—..........-.....-..........................—¦ C OPEH FILE FOR IHITIAL IHFORHATIOH AHD RESULTS C......--.......................................... C open (unit=l,file='co.dat',status='old') open (unit=2,file='el.dat',status='old') open (unit=4,file='topo.dat',status='neu') open (unit=7,iile='resul.dat',status='neu') open (unit=8,file='bound.dat',status='nen') READ ADVECTIVE VELOCITY print read* ;?,'EHTER ADVECTIVE VELOCITY' , alia READ HODAL C0ORDIMATES AHD BOUHDARY COHDITIOH CODE nf ixb = 0 readU.*) HUHHP do 10 K=I1HUHHP readd,*) nd, x(l,k), x(2,k) if ((x(l,k) .eq. 0.D0) .0R. (x(2,k) .eq. 0.DO)) then ia(k) = 1 nf ixb = nf ixb + 1 else ia(k) = 0 endif 10 continue C-------------------.......-.......-----.......-----.......— C READ IH FIXED BOUNDARY CONDITIONS C C NEGATIVE NUHBER IH IA FOR FIXED BOUNDARY CONDITIOH C C.............--............------------------.......-....... C ifix = 0 DO 30 nd =1, NUMNP if (ia(nd) .eq. 1) then ifix = ifix + 1 if (x(l,nd) .eq. 0.D0) then bound(ifix) = 1.D0 else bound(ifix) = 0.DO endif if ( abs(bound(ifix)).gt.l.e-8 ) ia(nd) = -ifix »rite (8,2030) HD ,bound(ifn) 2030 formatUB,1Ox,flO.6) endif 30 CONTINUE C C NUHBERING OF EQUATIONS A.2. LE CODE STILS. e SEq = o DO 150 J=I1KUHKP IF (IA(J) .LT. 0) GO TO ISO IF (IA(J) .GT. O) GO TO 140 NEQ = NEq + 1 IA(J) = NEq GO TO ISO 140 IA(J) = O IBO CONTINUE C C READ IN ELEMENT CONNECTIONS C.......-.................-.............-..............-...... C iel = 4 read(2,») NELEH DO 180 N=I1NELEH READ (2,*) ND1(IEK(I1N)1I=LIEL) 180 CONTINUE C C---.....-........-............--.............-.............. C FIND THE STRUCTURE OF THE HATRIX C......-............................................-......... C.......--.................CODE SYHETRIE HATRICE.....—....... NSYH = 0 CALL ICLEAR (HDIAG1KEq) CALL FRUSKI (IA1IEN,1,IEL,HDIAG,NELEH) CALL FSTSKI (NSYH1HDIAG1NEq1NA) C C-----......---..........-......-........-................... C COMPUTE AND ASSEHBLE LOCAL HATRICES AND RIGHT-HAKD SIDES C---.........--.................-.....-........-.....--.....-- C CALL CLEAR (A1NA) CALL CLEAR (B1NEq) C DO EOO N=I1NELEH C C FIKD COORDIHATES DEFINING ELEMENT NE AND C ELEMENT COKNECTION VECTOR LH. C DO 170 I=I1IEL ND = IEN(I1N) XX(I) = X(I1ND) YY(I) = X(2,ND) LH(I) = IA(ND) 170 CONTINUE C C CONSTRUCT ELEHEKT STIFFNESS HATRIX C 94 ANNEXE A. CALL ESTIFF (alfa,STIFF,BE,LH,IEL,XX.YY) C C ELIHIHATE FIXED BOUNDARY VALUES AT ELEHEST LEVEL. C CALL FXVGHG (STIFF,BE,LH,IEL,IEL,BOUHD) C C ASSEHBLE TOTAL STIFFHESS HATRIX AHD RIGBT-BAHD SIDE C CALL ADDSPI (3,A,B,HDIAG,STIFF,BE,LH,IEL,IEL) C 500 CONTIHUE C C--------......-----.....----------.....-..........------........ C SOLVE THE LIHEAR SYSTEH OF EQUATIONS C-----------.....-.........-----.....--..........-............... C CALL SLVSPI (A,B,B,HDIAG,HEq,1,IERR) C C----------.....--..........----.....-........-.....------....... C OUTPUT SOLUTIOH C...........-.....-.....-.......------.....-----.....------..... C DO 200 I=I1HUHHP INT = IA(I) IF (IA(I) .GT. O) THEN Y(I) = B(IHT) ELSE IF (IA(I) .LT. 0) THEN Y(I) = BOUHD(-IHT) ELSE IF (IA(I) .Eq. O) THEN Y(I) = O. EHD IF 200 COHTINUE LIG = nint( sqrt(dble(numnp))) WRITE (7,2070) ( (Y(I),I=l+(J-i)*LIG ,J*LIG),J=I1LIG) 2070 F0RHATUKF7.3)) C--...........----FERHETURE DES FICHIERS------------............. close(l) close(2) close(4) close(7) close(8) END C C---..............--..........--------.....--------------------- C------.....-.........---.....--.....--......----------......-- C SUBROUTIHE ESTIFF (A,STIFF,BE,LH,IEL,R,Z) C----------.......----........-------------------------......... C COHTRUCT ELEHENT STIFFNESS AND RIGHT-HAHD SIDE C----------...........----------......-.....-------.........---- C A.2. LE CODE STILS. 95 IHPLICIT REAL*8 (A-H1O-Z) PARAMETER (EPS=I.E-è) DIHENSIOH STIFF(4.4),BE(4),LH(4),R(»),Z(*) DIHENSION SHAPE(4),DADL(4,2),DADG(4,2),GAUS(2) DATA GAUS /.57735026918963,-.67735026918963/ C CALL CLEAR (BE,IEL) CALL CLEAR (STIFF,IEL»IEL) C C INTEGRATION ON 2 X 2 POINTS C DO 400 J=I,2 DO 400 JJ=I,2 C C FIND SHAPE FUNCTIONS LOCAL AND GLOBAL DERIVATIVES AND C JACOBIAN OF TRANSFORMATION , AT POINTS GAUS(2) IN THE C REFERENCE ELEMENT. C CALL Q4SHAP (R,Z,GAUS(J),GAUS(JJ),IEL,SHAPE,DADL,DADG,DET,-3 1 ,EPS,IERR) DO 300 K = 1,IEL DO 200 L = 1,IEL STIFF(K,L) = STIFF(K,L) + DET*(DADG(K,2)*DADG(L,2) 1 + A*(DADG(K,2)*DADG(L,1) 2 + DADG(K,1)»DADG(L,2) 3 + A»DADG(K,1)*DADG(L,1))) 200 CONTINUE C C STORE FORCES IN BE C IF (LM(K) .LE. O) GO TO 300 BE(K) = BE(K) + DET*SHAPE(K)«(0.D0) 300 CONTINUE 400 CONTINUE END Bibliographie |1] S. N. Antontsev, A. V. Kazhikhov, and V. N. Monakhov, Boun- dary value problems in mechanics of nonhomogeneous fluids, North- Holland Elsevier, 1990. [2| V. Arnold, Équations différentielles ordinaires, Mir, Moscou, 1974. |3) A. AsSEMIEN, Comportement asymptotique des équations de Navier- Stokes pour des écoulements de faible épaisseur, PhD thesis, Université Claude Bernard Lyon I, 1993. (4] A. Assemien, G. Bayada, and M. Chambat, Inertial effects in the asymptotic behavior of a thin film flow, tech. rep., Equipe d'Analyse Numérique Lyon-St Etienne, 1991. [5] P. AzÉRad, Analyse et approximation du problème de stokes dans un bassin peu profond, CR. Acad. Sci. Paris, 318 (1994), pp. 53-58. |6] P. AzÉRAD, P. Perrochet, AND J. PousiN, Space-time integrated least-squares : a simple, stable and precise finite element scheme to solve advection equations as if they were elliptic, in Journées de Metz 1995, M. Chipot, ed., Pitman, to appear. |7j P. Azérad and J. PousiN, Inégalité de poincaré courbe pour le traite- ment variationnel de l'équation de transport, soumis à CR. Acad. Sci. Paris, (1995). |8) I. Babuska, Error bound for the finite element method, Num. Math., 16 (1971), pp. 322-333. |9] C. Bardos, Problèmes aux limites pour les équations aux dérivées par- tielles du premier ordre à coefficients réels; théorèmes d'approximation; applications à l'équation de transport, Ann. scient. Ec. Norm. Sup., 3 (1970), pp. 185-233. 11OJ O. Besson AND M. R. LaYDI, Some estimates for the anisotropic navier-stokes equations and for the hydrostatic approximation, M2AN - Mod. Math. Ana. Num., 26 (1992), pp. 855-865. [11] O. Besson, M. R. Laydi, and R. Touzani, Vn modèle asymptotique en océanographie, CR. Acad. Sci. Paris, 310 (1990), pp. 661-665. 97 98 BIBLIOGRAPHIE [12] J. L. Borges, Fictions, Gallimard. [13] F. Brezzi, On the existence,uniqueness and approximation of saddle- point problems arising from lagrangian multipliers, R.A.l.R.O. Ana. Num., 8 (1974), pp. 129-151. [14] F. BREZZI and M. Fortin, Mixed and hybrid finite element methods, Springer-Verlag, 1991. [15] A. N. Brooks and T. J. R. Hughes, Streamline upwind/ petrov- galerkin formulations for convective dominated flows with particular emphasis on the incompressible navier-stokes equations, Comput. Meths Appi. Mech. Engrg., 32 (1982), pp. 199-259. [16] P. G. Ciarlet, The finite element method for elliptic problems, North- Holland Elsevier, 1978. [17] V. Cominciou, Fortran 77, introduzione e applicazioni numeriche, McGraw-Hill, 1991. [18] R. Dautray and J. L. Lions, Analyse mathématique et problèmes aux limites, vol. 9, Masson, 1988. [19] G. de Rh am, Variétés differentiates, formes, courants, formes har- moniques, Hermann, 1960. [20] R. J. Diperna and P. L. Lions, Equations différentielles ordinaires et equations de transport avec des coefficients irreguliers, in Séminaire Equations aux dérivées partielles 88-89, Ecole Polytechnique, Palaiseau, 1989. [21] R. J. DiPerna and P. L. Lions, Ordinary differential equations, transport theory and sobolev spaces, Invent, math., 98 (1989), pp. 511- 547. [22] G. Geymonat and P. Leyland, Transport and propagation of a per- turbation of a flow of a compressible fluid in a bounded region, Arch. Rat. Mech. Anal., 103 (88), pp. 53-81. [23] V. Girault and P. A. Raviart, Finite element methods for Navier- Stokes equations, Springer-Verlag, 1986. [24] R. Glowinski, Numerical methods for nonlinear variational problems, Springer-Verlag, 1984. [25] P. Hansbo, The characteristic streamline diffusion method for convection-diffusion problems, Comput. Meths Appi. Mech. Engrg., 96 (1992), pp. 239-253. [26] E. Hopf, über die anfangswertaufgabe fur die hydrodynamischen grundgleichungen, Math. Nachr., 4 (1951), pp. 213-231. BIBLIOGRAPHIE 99 (27) L. HeRMANDER, Non-linear Hyperbolic Differential Equations, Depart- ment of Mathematics, University of Lund, 1987. [28] T. J. R. Hughes, L. P. Franca, and G. M. Hulbert, A new fi- nite element formulation for computational fluid dynamics: Viii. the galerkin/least-squares method for advective-diffusive equations, Corn- put. Meth. Appi. Mech. Engrg., 73 (1989), pp. 173-189. (29] C. JOHNSON, Numerical solution of Partial Differential Equations by the Finite Element Method, Cambridge University Press, 1987. [30] ------, A new approach to algorithms for convection problems which are based on exact transport + projection, Comput. Meths Appi. Mech. Engrg., 100 (1992), pp. 45-62. [31] A. Kufner, O. John, and S. Fucik, Function spaces, Noordhoff, Leyden, 1977. (32] O. A. Ladyzhenskaya, The mathematical theory of viscous incom- pressible flow, Gordon Breach, 2nd ed., 1969. [33] J. Leray, Aspects de la mécanique théorique des fluides, C. R. Acad. Sci. Paris, 11 (1994), pp. 287-290. (34] P. Lesaint, Introduction aux systèmes de friedrichs et à V équation de transport, notes de cours. (35] J. L. Lions, Quelques méthodes de résolution des problèmes aux limites non linéaires, Dunod, 1969. [36] ------, Perturbations singulières dans les problèmes aux limites et en contrôle optimal, Lecture Notes in Mathematics 323, Springer, 1973. [37] ------, On some problems connected with navier-stokes equations, in Non- linear Evolution Equations, M. Crandall, ed., Academic Press, 1978, pp. 59-84. (38) P. L. Lions, Mathematical Topics in Fluid Mechanics, CEREMADE Université Paris-Dauphine, 1995. preliminary version. (39) J. Necas, Sur une méthode pour résoudre les edp du style elliptique voi- sine de la variationnelle, Ann. Sc. Norm. Sup. Pisa, 4 (1962), pp. 305- 326. [40] H. Nguyen and J. Reynen, A space-time least squares finite element scheme for advection- diffusion equations, Comput. Meths. Appi. Mech. Engrg., 42 (1984), pp. 331-342. (41] J. T. ODEN and J. N. Reddy, Variational methods in theoretical mechanics, Springer-Verlag, 1983. 100 BIBLIOGRAPHIE [42] J. Pedloskv, Geophysical fluid dynamics, Springer-Verlag, 1987. (43] P. Perrochet and P. Azérad, Space-time integrated least-squares: Advection = anisotropic diffusion, in Journées de Besançon 1995, J. Crolet, ed., Pitman, to appear. (44] ------, Space-time integrated least-squares: Solving a pure advection equation with a pure diffusion operator, J. Comput. Phys., 117 (1995), pp. 183-193. [45] L. Quartapelle, Numerical Solution of the Incompressible Navier- Stokes equations, Birkhaeuser, 1993. [46] J. E. Roberts and J. M. Thomas, Handbook of numerical analysis, in Mixed and hybrid methods vol II, Finite Element Method (part 1), North-Holland Elsevier, 1991, pp. 527-637. [47] E. Sarasin and L. du Pasquier, Les seiches du lac de neuchâtel, Bull. soc. sc. nat. Neuchâtel, (1895), pp. 3-9. [48] P. J. Shopov, Basis for finite element schemes for an inhomogeneous fluid, Math. Balkan., 4 (1990), pp. 113-128. [49] J. Simon, Compact sets in the space lp(o,t;b), Ann. Mat. Pura Appi., 146 (1987), pp. 65-97. (SO] ------, Sur les fluides visqueux incompressibles et non homogènes, CK. Acad. Sci. Paris, 309 (1989), pp. 447-452. [51) ------, Non-homogeneous viscous incompressible fluids: existence of ve- locity, density and pressure, SIAM Journal of Mathematical Analysis, 21 (1990), pp. 1093-1117. [52] ------, Sobolev, besov and nikolskii fractional spaces: imbeddings and comparisons for vector valued spaces on an interval, Ann. Mat. Pura Appi., 157 (1990), pp. 117-148. [53] R. Stenberg, On some three-dimensional finite elements for incom- pressible media, Comput. Meth. Appi. Mech. Engrg., 63 (1987), pp. 261- 269. |54] ------, Error analysis of some finite element methods for the stokes pro- blem, Math. Comp., 54 (1990), pp. 495-508. [55] G. Strang, Linear Algebra and its application, Harcourt Brace Jova- novich, San Diego, 3rd ed., 1998. (56] U. Suter, communication personnelle. [57] R. Temam, Sur la stabilité et la convergence de la méthode des pas fractionnaires, Ann. Mat. Pura ed Applicata, (1968), pp. 191-379. BIBLIOGRAPHIE 101 [58] ------, Navier-Stokes equations, North-Holland Elsevier, 1985. (59] A. Valli, An existence theorem for non-homogeneous inviscid incom- pressible fluids, in Differential Equations, C. Dafermos and al., eds., Marcel Dekker, New York, 1989, pp. 691-698. [60] A. Valli and W. M. Zajaczkowski, About the motion of nonhomo- geneous ideal incompressible fluids, Nonlinear Anal., 12 (1988), pp. 43- 50. (61 ) W. VelTE, Direkte methoden der Variationsrechnung, Teubner, Stutt- gart, 1976. |62] R. K. Zeytounian, Modélisation asymptotique en mécanique des fluides newtoniens, Springer-Verlag, 1994. 102 BIBLIOGRAPHIE Les mathématiques sont simples, c'est nous qui sommes compliqués. Dr Teddy 103