Solveur éléments finis quatre champs pour des écoulements quasinewtonien

Télécharger le fichier original (Mémoire de fin d’études)

La glace : histoire d’un fluide rampant

L’historique de cette section s’inspire de e.g. Rémy et Testut [149] et Blatter et al. [22]. On renvoie à ces références pour de plus amples détails sur l’histoire de la description de la dynamique de la glace.
Le terme “fluide” est employé en 1773 pour la première fois pour qualifier la glace par le naturaliste genévois André-César Bordier [23]. Il écrit que “la masse entière des glaces est liée ensemble et pèse l’une sur l’autre de haut en bas à la manière des fluides”. Cependant l’idée ne sera pas retenue. Une théorie longuement en vogue, proposée par Johann Scheuchzer en 1705 et reprise par Jean de Charpentier en 1830, explique le mouvement des glaciers par un système de fonte-regel de la glace à la base du glacier. En gelant, l’eau se dilate et soulève le glacier, provoquant ainsi son déplacement. Parallèlement, une autre théorie, avancée par Gotlieb Gruner en 1760 et reprise plus tard par Horace de Saussure, attribue le mouvement des glaciers à un glissement sur leur base. En 1840, le suisse Louis Agassiz défend ardemment les théories de de Charpentier sur l’avancée par dilatation tout en rejetant avec force les théories de de Saussure sur le glissement basal. Une variation sur ce thème de fonte-regel, conciliant le paradoxe apparent du caractère à la fois solide et ductile de la glace, proposée par John Tyndall et Michael Faraday en 1873 et soutenue par des faits expérimentaux, deviendra la seule théorie consensuelle durant plusieurs décen-nies. Cette diﬃculté à concevoir la déformation de la glace est principalement due à la distinction formelle entre un corps solide et un corps fluide, qui repose sur des considérations mécaniques de comportement en traction et cisaillement nécessaire-ment associées à une échelle de temps “humaine”.
En 1920, Weinberg et Somigliana développent la théorie mathématique du fluage (déformation au cours du temps d’un matériau soumis à une charge fixe) où la glace est considérée comme un fluide newtonien très visqueux. Cependant, un certain nombre d’anomalies, pointées par des mesures in-situ (notamment une grande sen-sibilité à une faible variation de contraintes où bien une variabilité importante des vitesses d’écoulement à l’échelle de la journée) soulevaient de nombreuses questions. Le métallurgiste Orowan et le physicien des solides Perutz incorporent en 1948 le concept de plasticité dans la description de la dynamique de la glace. Classique en physique des solides, ce comportement stipule que le corps ne se déforme qu’au delà d’une contrainte seuil à partir de laquelle il se comporte comme un fluide parfait (non visqueux). Cependant, des objections sont émises sur la validité d’une des-cription purement plastique du mouvement de la glace ; selon le type de glaciers, le comportement semblerait plutôt visqueux ou plutôt plastique.
Le débat prendra finalement fin avec John Glen qui, grâce à un protocole ex-périmental rigoureux, observe un comportement intermédiaire entre fluide visqueux newtonien et solide plastique. Dans le célèbre article [59], Glen énonce la loi reliant contrainte et déformation sous la forme d’une loi de puissance. Depuis 1955, la des-cription de la dynamique de la glace par les équations de la mécanique des fluides couplée à la loi de comportement en loi de puissance énoncée par Glen n’a jamais été fondamentalement remise en cause bien que de nombreuses autres propriétés rhéologiques du matériau glace aient été considérées dans les descriptions. La théo-rie du glissement basal est énoncée par Weertman [172] qui propose une description du glissement sur un socle rigide, description encore largement dominante dans la glaciologie moderne.

Rappels de mécanique des milieux continus

Les écoulements de fluide considérés ici sont décrits en utilisant les concepts de la mécanique des milieux continus. L’un des objectifs basique de la mécanique des milieux continus est de décrire comment un matériau continu soumis à une force se déforme et comment la réponse à cette contrainte peut être décrite par une relation constitutive liant la contrainte à la déformation.

Déformations

On se place en description lagrangienne, c’est à dire que l’on considère l’existence d’une fonction ϕ qui à une configuration initiale Ω0 du milieu continu associe l’état courant Ω(t). Autrement dit :
ϕ : Ω0 → Ω(t)
X ϕ(X, t) = x(t)
où X représente les coordonnées initiales de Lagrange d’une particule P du mi-lieu continu et x(t) représente la position de la particule à l’instant t. Dans ce cas, on notera d(X, t) = ϕ(X, t) − X le vecteur déplacement et v(x, t) = ∂x(t) la vitesse ∂t d’une particule.
Un milieu continu subit des déformations si les distances relatives entre les parti-cules varient au cours du temps. Ainsi, considérant un vecteur dX = X′ −X représen-tant la distance entre deux particules P et P ′ à l’instant initial, on s’intéresse alors à la position x′ de la particule P ′ à l’instant t. On obtient, à l’aide d’un développement de Taylor : x′ = ϕ(X ′, t) = ϕ(X + dX , t) = ϕ(X , t) + ∇X ϕ(X , t)dX + |dX |ε(dX )
On définit ainsi la matrice gradient de la déformation F : F = ∇X ϕ(X , t)
À l’aide de F , on peut définir le tenseur des dilatations C (ou tenseur de Cauchy-Green à droite) : C=FTF et, par suite, le tenseur des déformations E (ou tenseur de Green-Lagrange) s’exprime ainsi : E = 12 (C − I d)
La nullité du tenseur des déformations décrit un mouvement de corps rigide.
Dans un cadre de description fluide, on se place sous l’hypothèse des petites déformations soit : d ≪ 1. On a donc : 1 )T (∇X d + 1 (∇X dT + ∇X d) + 1 ∇X dT ∇X d E (d) = (∇X d + I d I d ) − I d = 2 2 2
Au premier ordre, on a donc : E (d) ≃ 1 (∇X d + ∇X dT ). On introduit alors le tenseur des déformations linéarisés ε : ε(d) = 12 (∇X d + ∇X dT )
À noter que les composantes diagonales εii représentent une approximation de l’allongement unitaire dans la direction de base ei tandis que les composantes hors-diagonale εij représentent une approximation du demi-angle de glissement entre les directions de base ei et ej (cf. e.g. [90]). La variation temporelle ε˙ de ε(d) définie par ε˙ = ∂∂tε engendre le tenseur des vitesses de déformation (ou tenseur des taux de déformation) D(v) : ε˙ = D(v) = 12 (∇X v + ∇X vT )

Contraintes : la méthode de Cauchy

On considère que les eﬀorts que subit un sous-domaine matériel ω(t) se classent en deux familles ; d’une part, les eﬀorts volumiques (ou eﬀorts de champ) s’appli-quant à tout point P du domaine (e.g. la gravité) et d’autre part, les eﬀorts sur-faciques (ou eﬀorts de contact) qui s’appliquent en tout point de la frontière ∂ω de ω.
Les eﬀorts surfaciques peuvent prendre la forme d’une action normale et/ou tan-gentielle. On les représente donc par une densité surfacique d’eﬀort s(x, ∂ω, t). Le vecteur s est appelé vecteur contrainte.
Considérer la dépendance de s à toutes les caractéristiques de la frontière du domaine ∂ω est ambitieux. On fera donc appel, de manière classique, au postulat de Cauchy sur la dépendance des eﬀorts de contact vis-à-vis des caractéristiques de la surface : La force de contact F au point x de la frontière de ω ne dépend que de la normale n à ∂ω au point x soit : F = F (n).
Ce postulat définit le caractère local des forces de contact.
De ce postulat découle le théorème de Cauchy (voir par exemple [90]) qui dit que l’application de R3 dans R3 qui à n associe F (n) est linéaire. On note σ une telle application soit :
[σ]n = F (n)
σ s’appelle le tenseur des contraintes de Cauchy. En 3D, il s’écrit :
σ11 σ12 σ13
σ = σ21 σ22 σ23 σ31 σ32 σ33
On le réprésente classiquement sur un cube élémentaire (voir figure 2.1). Dans ce cas, les composantes diagonales de σ sont perpendiculaires aux facettes du cube et sont les contraintes de compression (ou de traction suivant le signe) aussi appelées contraintes normales.
Les composantes hors-diagonale sont tangentielles à la facette considérée et cor-respondent aux contraintes de cisaillement.
Le tenseur des contraintes étant symétrique, il est diagonalisable : il possède trois valeurs propres réelles σi appelées contraintes principales et trois vecteurs propres ni, orthogonaux deux à deux, appelés directions principales de contraintes : σ = σini ⊗ ni, avec σ • ni = σini i=1

Fluide visqueux

Par définition, un fluide est un milieu qui ne peut pas résister aux contraintes de cisaillement qu’on lui applique c’est à dire que sa résistance ne peut empêcher la déformation et le mouvement.
Le tenseur des contraintes, comme tout tenseur, se décompose en une partie sphérique isotrope (ou contrainte moyenne) correspondant à la pression isostatique, et une partie déviatoire non isotrope contribuant aux contraintes tangentielles et correspondant à l’écart de la contrainte par rapport à la pression isostatique : σ = −pI d + S (2.1)
L’existence du déviateur des contraintes S est entièrement due à l’existence d’un mouvement du fluide (cf. e.g. [14]).
Sous l’hypothèse des petites déformations, on peut considérer que le déviateur des contraintes S est une fonction linéaire des composantes du gradient de vitesse. Cela caractérise la linéarité locale de la réponse du milieu continu. Autrement dit : Sij = Aijkl Dkl (2.2) où Dkl représente les coeﬃcients du tenseur des taux de déformations et Aijkl , les coeﬃcients d’un tenseur d’ordre 4, symétrique. Le tenseur A représente alors le tenseur de viscosité. Si l’on considère un fluide isotrope, on peut écrire la relation suivante σ = (−p + λ tr(D))I d + 2ηD
Dans cette relation, η est appelé le coeﬃcient de viscosité apparente (shear visco-sity) et (λ + 23η ) la viscosité de volume (bulk viscosity). La contrainte s’exprimant en PA et les taux de déformations en S−1. On en déduit que les coeﬃcients de viscosité η et λ sont homogène à des PA • S.

Lois de conservation

Conservation de la masse Un écoulement est dit incompressible quand la densité d’un élément fluide n’est pas modifiée par les variations de pression au cours du mouvement. Cela revient à avoir une variation de la densité en fonction du temps nulle. L’équation de conservation de la masse qui en découle s’écrit alors : ∇ • u = 0 (2.3)
Dans le cas incompressible, cela implique : tr(D) = 0, ce qui nous amène à la relation constitutive suivante : S=2ηD (2.4)
Dans le cas d’un fluide idéal (i.e. incompressible et non visqueux), la contrainte est simplement déterminée par la pression isostatique : σ = −pI d
Dans le cas où la relation linéaire (2.2) entre le déviateur des contraintes et les taux de déformations est valable, η est alors une constante de proportionnalité entre contrainte et déformation et on parle d’un fluide newtonien.
Conservation de la quantité de mouvement La loi de conservation de la quantité de mouvement s’obtient en appliquant la relation fondamentale de la dynamique à un volume de fluide, i.e. en écrivant l’égalité entre la variation temporelle de la quantité de mouvement et l’ensemble des forces (de volumes et de surface) exercées sur ce volume. Plus précisément, pour un fluide de masse volumique ρ se déplaçant à une vitesse u, soumis à des contraintes σ décrite par (2.1) et des forces de volume f , on obtient : ρ ∂u + (u • ∇)u = ρf − ∇p + ∇ • S (2.5) ∂t
L’introduction de la relation (2.4) décrivant les contraintes déviatoires pour un fluide visqueux isotrope et incompressible dans la relation (2.5) amène les équa-tions de Navier-Stokes. On introduit par la suite la notion de loi de comportement amenant une modification mineure dans la relation (2.4) en prenant en compte la dépendance expérimentalement observée de la viscosité η au taux de cisaillement. L’équation résultante est très proche des équations de Navier-Stokes classique et décrit la dynamique d’un fluide newtonien généralisé ou quasi-newtonien.
On fera par la suite une approximation des équations de Navier-Stokes associée au écoulements à bas Reynolds menant au problème de Stokes. Cette approximation est introduite dans le cadre de la formulation instationnaire Arbitraire Lagrange Euler présentée en section 2.6.2.

Loi de puissance, loi de Glen et modélisation thermique

Loi de comportement

Dans le cas d’un corps déformable, les lois de conservation du problème se ré-sument aux bilans locaux de masse (équation de continuité) et de quantité de mou-vement (équation d’équilibre) soit une équation scalaire et une équation vectorielle. D’autre part, les inconnues du problème sont la transformation ϕ(X , t) et le tenseur des contraintes σ = (σij ) symétrique ce qui représente 9 fonctions scalaires incon-nues. Dans le cadre d’un problème stationnaire dont on néglige la thermodynamique, l’équation de continuité est automatiquement vérifié et on obtient 3 équations sca-laires pour 9 inconnues. En l’état, le système est donc mal posé. Il reste 6 équations à trouver correspondant aux expressions des 6 composantes des contraintes σ. De telles équations ne sont pas universelles et prennent en compte la spécificité du ma-tériau considéré. La relation (2.4) est la loi de comportement décrivant un fluide newtonien de viscosité constante η.
Une loi de comportement est l’écriture mathématique de certaines hypothèses faites sur le comportement mécanique du matériau considéré. Il en découle une nécessité d’invariance de la loi de comportement dans le cas d’un changement de repère et/ou d’un changement de référentiel. Cette nécessité peut être vue comme un postulat ou admise intuitivement et est dénommée principe d’objectivité matérielle (voir e.g. [6]).
Invariants d’un tenseur On appelle invariants d’un tenseur des fonctions à valeurs réelles des composantes d’un tenseur qui ne dépendent pas du choix de la base de dé-composition. De tels invariants existent pour tout endomorphisme et correspondent à des fonctions des valeurs propres de l’endomorphisme considéré. Dans le cas d’un tenseur symétrique, on est assuré que ses valeurs propres sont réelles.
Dans le cas du tenseur des taux de déformations (symétrique) D, si l’on note ses trois valeurs propres λi, 1 ≤ i ≤ 3, voici les expressions classiquement utilisées :
I1 = λ1 + λ2 + λ3 = tr(D)
I2 = λ1λ2 + λ2λ3 + λ1λ3 = 12 ((tr(D))2 − tr(D2))
I3 = λ1λ2λ3 = det(D)
Ces trois invariants sont appelés invariants principaux de D.
Fluides non newtoniens Comme vu précédemment, un fluide newtonien est carac-térisé par l’hypothèse de linéarité et d’isotropie du déviateur des contraintes par rapport aux composantes du gradient des vitesses. Cependant, on observe expéri-mentalement que de nombreux fluides voient leur viscosité apparente η évoluer avec l’augmentation des contraintes de cisaillement.
Des fluides dont la viscosité apparente diminue avec l’augmentation de la contrainte sont qualifiés de rhéofluidifiants (shear-thinning) ou pseudo-plastiques. Pour décrire ce type de comportement et prendre en compte les eﬀets rhéofluidifiants, on ne fait plus l’hypothèse d’une relation linéaire entre déviateur des contraintes et vitesses de déformations comme dans la relation (2.2). Une généralisation de cette descrip-tion consiste à choisir une loi de comportement de la forme S = g(D). Le principe d’objectivité matérielle appliqué à cette relation impose que le matériau considéré soit isotrope. L’isotropie de la fonction g et la symétrie du tenseur D amènent, par un théorème d’analyse tensorielle, à décrire la viscosité η comme une fonction des invariants du tenseur des vitesses de déformation D.

Table des matières

Introduction générale
1 Introduction
1.1 Historique et généralités sur la modélisation de l’écoulement des glaces terrestres
1.2 Problématique d’échelle
1.3 L’élévation du niveau des mers
1.4 Synthèse des problématiques et annonce du plan
2 Modélisation directe des écoulements de glace
2.1 La glace : histoire d’un fluide rampant
2.2 Rappels de mécanique des milieux continus
2.2.1 Déformations
2.2.2 Contraintes : la méthode de Cauchy
2.2.3 Fluide visqueux
2.3 Lois de conservation
2.4 Loi de puissance, loi de Glen et modélisation thermique
2.4.1 Loi de comportement
2.4.2 Modélisation thermique
2.5 Modélisation du socle : revue des lois de glissement
2.5.1 Glissement sur un socle rigide : la théorie de Weertman
2.5.2 Glissement sur un socle déformable
2.5.3 Bilan
2.6 Modèle d’écoulement de fluide géophysique quasi-newtonien à surface libre
2.6.1 Le modèle d’écoulement
2.6.2 Formulation du problème instationnaire
3 Problème inverse – Approche variationnelle
3.1 La méthode de l’état adjoint
3.1.1 Position du problème
3.1.2 Fonction coût
3.1.3 Modèle linéaire tangent
3.1.4 Modèle adjoint
3.1.5 Analyse de sensibilité locale
3.2 L’assimilation variationnelle de données
3.2.1 Problème de contrôle
3.2.2 Système d’optimalité
3.2.3 Méthodes locales et régularisation
3.3 Implémentation du modèle adjoint
3.3.1 Les différentes stratégies
3.3.2 Différentiation automatique
3.3.2.1 Principes
3.3.2.2 Validation du code adjoint
3.3.3 Résumé de l’algorithme
3.4 Utilisation et rôle des méthodes inverses en glaciologie
3.4.1 Les données en glaciologie
3.4.2 Des données au modèle et du modèle aux données
3.5 Le logiciel DassFlow-Ice
3.5.1 Structure du code de calcul
3.5.2 Une première étude : le glacier Variegated
4 Solveur éléments finis quatre champs pour des écoulements quasinewtonien
4.1 Résumé en français
4.2 Introduction
4.3 The Fluid Model
4.4 Basic Numerical Procedures
4.4.1 Mixed Weak Formulation
4.4.2 Classical Solvers and Limitations
4.5 Four-Field Saddle Point Problem
4.5.1 Four-Field Saddle-Point Formulation
4.5.2 Existence and Uniqueness of the Saddle-Point
4.5.3 Characterization of the saddle-point
4.5.4 LA and LA algorithms
4.5.5 Finite element discretization
4.6 Performances of LA and LA Algorithms
4.6.1 Convergence curves
4.6.2 Comparative performances
4.6.3 Summary
4.7 Local sensitivity analysis
4.7.1 Observations and cost function
4.7.2 Adjoint model
4.8 Numerical experiments
4.8.1 Analytical Poiseuille-like solution
4.8.2 The viscoplastic steady wave
4.8.2.1 Description
4.8.2.2 Numerical results
4.8.3 Geophysical test case : the Mertz glacier
4.8.3.1 Description
4.8.3.2 Numerical results
4.9 Concluding remarks
5 Analyse de sensibilité variationnelle et rhéométrie inverse
5.1 Résumé en français
5.2 Introduction
5.3 Continuous and Discrete Forward Model
5.4 Numerical Inverse Tools
5.4.1 Basic Principles of a Variational Data Assimilation Process
5.4.2 Cost Function in Presence of Observations
5.4.3 Adjoint Model and Gradient Computation
5.4.3.1 Local Sensitivity Analysis
5.4.3.2 Data Assimilation and Twin Experiments
5.5 Analytical Solution and Related Sensitivities
5.6 Are Surface Velocity Observations Reliable ?
5.7 Virtual Rheometry : Power-Law Exponent
5.7.1 Flow Description
5.7.2 Numerical Results
5.8 Virtual Rheometry : Consistency η0 (Temperature Depen- dency)
5.8.1 Description of the Thermal Dependent Flow
5.8.2 Sensitivity With Respect to η0
5.8.2.1 No-Slip Case
5.8.2.2 Sliding Case
5.8.2.3 Discussion
5.8.3 Identification of η0 Based on Surface Velocity Observations
5.8.3.1 Naive Approach
5.8.3.2 Approach Based on Prior Physical Considerations
5.9 Conclusion
6 Sur la précision du gradient et la nécessité d’un adjoint exact
6.1 Introduction
6.2 Adjoint exact, méthode de l’accumulation retour et ap-
proximation « auto-adjointe »
6.2.1 Précision de l’approximation « auto-adjointe »
6.2.2 Troncature de l’accumulation retour
6.3 Identifiabilité du coefficient de friction
6.3.1 Écoulement quasi-uniforme
6.3.2 Écoulement sur topographie réelle du glacier Mertz
6.4 Densité des données
6.5 Conclusion
7 Conclusion générale et perspectives
Bibliographie
A Adjoint of a linear solver
A.1 The direct routine
A.2 The linear tangent routine
A.3 The adjoint routines generated
A.4 The adjoint of the linear system
B Rapport d’analyse numérique pour le logiciel DassFlow-3D
B.1 Weak formulation
B.2 Finite element formulation
B.2.1 Galerkin approximation
B.2.2 The Taylor-Hood (Lagrangian P2/P1) finite element
B.2.3 The friction condition
B.3 Validation
B.3.1 Validation of the bidimensional code
B.3.1.1 Analytical test-case
B.3.1.2 ISMIP test-cases
B.3.2 Explicit solution in 3D
B.3.3 Validation of the friction condition