Applications à l’érosion
Les méthodes numériques implémentées tout au long de cette thèse sont assemblées pour simuler l’érosion interfaciale d’un sol sous l’effet d’un écoulement fluide incompressible. La modélisation décrite précédemment sera appliquée à l’érosion d’un sol cohésif par des écoulements à faible nombre de Reynolds ou modéré. Rappelons que dans le modèle numérique que nous avons présenté dans le premier chapitre, les calculs se déroulent en deux étapes : 1. L’étape de l’écoulement dans laquelle une méthode de pénalisation est utilisée pour calculer les équations de Navier-Stokes autour d’un obstacle avec une méthode aux domaines fictifs pour s’affranchir d’un maillage dépendant de la géométrie. 2. L’étape de l’érosion, où l’interface eau/sol est décrite par une fonction Level Set couplée à une loi d’érosion empirique à seuil. Cette loi est une relation linéaire entre la vitesse d’érosion et la contrainte de cisaillement exercée par le fluide sur le solide. L’objectif ici est de simuler le déplacement de l’interface fluide/solide durant l’érosion. L’évolution de cette interface est pilotée par la vitesse d’érosion qui dépend de la contrainte de cisaillement tangentiel du fluide. Étant donné que l’on s’attend à une discontinuité du gradient de la vitesse de l’écoulement à travers l’interface, des problèmes se poseront pour le calcul de la contrainte de cisaillement. La difficulté est de calculer d’une manière consistante la contrainte de ce problème. Nous mettrons en évidence cette difficulté et nous proposerons par la suite des solutions pour déterminer la contrainte qui doit être estimée correctement au niveau de l’interface. Cette difficulté est rencontrée pour les problèmes à frontière libre résolus sur grille fixe comme les problèmes de Stefan [25,45] et de Hele-Shaw . Enfin, la validation du modèle numérique et la pertinence des solutions proposées seront confirmées par plusieurs résultats de simulation.
Équations du modèle
Pour commencer, nous réécrivons d’abord les différentes équations qui décrivent les deux étapes énoncées dans le modèle, à savoir : l’écoulement incompressible et l’érosion interfaciale. Les milieux en jeux sont le fluide et le solide, qui occupent respectivement deux domaines notés Ωf et Ωs et séparés par une interface Γ = Ωf ∩ Ωf , isocontour 0 d’une Level Set φ En choisissant une approche Eulérienne sur maillage fixe avec des domaines fictifs, toutes les équations sont discrétisées dans le domaine global Ω = Ωf ∪ Ωf . Ainsi, le déroulement des calculs peut être effectué selon le schéma suivant :
Le calcul de l’écoulement : Le système d’équations de
Navier-Stokes décrivant l’écoulement incompressible que nous utilisons est le suivant (cf. chapitre 4) : ρ (∂tu − u ∧ rotu) − µ∆u + µ Ks H(φ)(u − us) = −∇p + f dans ]0, T[×Ω, div(u)=0 dans ]0, T[×Ω, (+C.L.) sur ]0, T[×∂Ω, (5.1) où u est le champ de vitesse défini sur Ω, tout comme la pression p, la densité ρ et la viscosité dynamique µ. Les conditions aux limites (C.L.) permettent de fermer le système et elles seront définies plus tard selon le cas étudié. La fonction Level Set est définie sur Ω telle que : φ < 0 dans le fluide Ωf , φ = 0 sur l’interface Γ, φ > 0 dans le solide Ωs. (5.2) Le paramètre de pénalisation Ks représentant une perméabilité est choisi tel que Ks ‘ 1. Ceci permet d’imposer la vitesse us dans le sol Ωs, caractérisé par φ > 0, en utilisant la fonction Heaviside : H(x) = ! 1 si x > 0 0 si x ≤ 0.
Le calcul de la contrainte
Le calcul de l’érosion commence par l’estimation de la composante tangentielle du vecteur contrainte du fluide à l’interface Γ : −→τ = 2µ[D(u) · −→n − ( −→n · D(u) · −→n ) −→n ], (5.4) où le tenseur D(u) représente le gradient symétrique de la vitesse d’écoulement : D(u) = 1 2 (∇u + ∇Tu), (5.5) et −→n la normale orientée dans le sens « fluide » vers « solide » et naturellement donnée par : −→n = ∇φ #∇φ# .
Le calcul de la vitesse d’érosion
L’érosion est pilotée par la célérité normale de l’interface vΓ. La loi d’érosion reliant cette vitesse à la composante tangentielle à l’interface Γ est la suivante : vΓ = ! Kd(τ − τc) si τ > τc, 0 sinon. (5.7) Où le coefficient Kd est un paramètre décrivant la cinétique de l’érosion et le seuil τc est la contrainte critique. La vitesse représentant l’érosion doit être étendu dans tout le domaine Ω pour faire évoluer toutes les lignes de niveau de la Level Set φ de la même manière. Pour cela, on utilise la technique d’extension des vitesses étudiée et validée dans le chapitre 2. Ainsi, la vitesse étendue vext Γ est la solution de l’équation de transport suivante : ! ∂tvext Γ + sgn(φ) ∇φ %∇φ% · ∇vext Γ = 0, vext Γ (t = 0, x) = vΓ. (5.8) Enfin, l’évolution de l’interface Γ est décrite par l’équation de transport de la Level Set φ : ! ∂tφ + vext Γ ∇φ %∇φ% · ∇φ = 0, φ(t = 0, x) = φ0. (5.9) Le calcul de l’écoulement incompressible, l’extension de la vitesse et le transport de la Level Set ont été étudiés et testés séparément. Le processus complet de l’érosion nécessite l’assemblage de ces calculs. L’estimation de la contrainte représente le chainon important reliant le calcul de l’écoulement et de l’érosion. Comme nous l’avons évoqué dans l’introduction de ce chapitre, avoir une bonne estimation pour la contrainte pose une grande difficulté étant donné que nous utilisons une approche eulérienne sur maillage fixe avec pénalisation. Puis nous mettrons en évidence ce problème, dans le cas HET, avant de proposer des solutions.
Gestion du problème à frontière libre
Le problème qui nous concerne est classiquement rencontré dans les problèmes à frontières libres, comme par exemple ceux posés par les écoulements dans des géométries élastiques qui obligent à résoudre des équations elliptiques sur des domaines en mouvement. Faire intervenir des maillages body-fitted [93] qui collent à la géométrie oblige à remailler les domaines au cours des itérations d’évolution en temps. Une solution très coûteuse en remaillage. En utilisant des maillages fixes, comme dans notre cas, les équations elliptiques peuvent être résolues sur un domaine approché. L’erreur de la solution est alors inconsistante pour les gradients au voisinage de l’interface sur le bord du domaine approché. On retrouve par exemple ce type d’inconsistance dans le problème de Hele-Shaw [7, 22] et dans le problème de Stefan [25, 45] décrivant la fonte des glaces par le transfert de chaleur dans un milieu homogène qui subit un changement de phase. Ces problèmes ont été résolus par exemple à l’aide des méthodes IIM (Immersed Interface Method) proposées par R.J. LeVeque et A.L. Li [62]. Dans ce cas, la continuité des propriétés à l’interface de discontinuité est assurée grâce à des prolongements de Taylor de la solution de chaque coté de l’interface. Mais pour garder la méthode de pénalisation sur un maillage fixe, des traitements spécifiques doivent être apportés