Une méthode des éléments finis utilisant une variable auxiliaire bornée pour résoudre l'équation de Richards

Informations fondamentales

• ID de l'article : 2510.13012
• Titre : A finite element method using a bounded auxiliary variable for solving the Richards equation
• Auteurs : Abderrahmane Benfanich (Université d'Ottawa), Yves Bourgault (Université d'Ottawa), Abdelaziz Beljadid (Université d'Ottawa & Université Polytechnique Mohammed VI)
• Classification : math.NA cs.NA physics.comp-ph
• Date de publication : 14 octobre 2025 (prépublication arXiv)
• Lien de l'article : https://arxiv.org/abs/2510.13012

Résumé

L'équation de Richards est une équation elliptique-parabolique non linéaire largement utilisée pour modéliser l'infiltration dans les milieux poreux. Cet article développe une méthode des éléments finis pour résoudre l'équation de Richards en introduisant une nouvelle variable auxiliaire bornée qui élimine les termes non bornés de la formulation faible. La méthode utilise un schéma semi-implicite pour la discrétisation temporelle et la méthode de Newton pour résoudre le système non linéaire. Cette approche élimine le besoin de techniques de régularisation et présente des avantages dans le traitement des régions sèches et complètement saturées. Dans la technique proposée, une méthode de décomposition de domaine de Schwarz non-chevauchante est utilisée pour modéliser l'infiltration dans les sols stratifiés. Des expériences numériques valident l'efficacité de la méthode proposée, incluant la modélisation de l'écoulement dans des fibres complètement sèches, les cas avec régions complètement saturées et sèches, ainsi que les problèmes d'infiltration dans les sols stratifiés.

Contexte et motivation de la recherche

Définition du problème

L'équation de Richards décrit l'écoulement non saturé dans les milieux poreux et constitue un modèle fondamental en hydrogéologie, génie environnemental et science des sols. L'équation existe sous plusieurs formes :

1. Forme en charge de pression (Ψ) : Applicable aux écoulements complètement saturés et partiellement non saturés, mais numériquement instable lorsque S→0
2. Forme en saturation (S) : Meilleure performance dans les sols secs, mais difficultés numériques dans la région saturée (S=1)
3. Forme mixte (S,Ψ) : Améliore les propriétés de conservation de masse, mais augmente la complexité

Limitations des méthodes existantes

1. Méthodes de commutation de variables principales : Peuvent produire des solutions non physiques aux interfaces saturé/non saturé
2. Techniques de régularisation : Nécessitent un choix supplémentaire de paramètres de modélisation
3. Méthodes des éléments finis traditionnelles : Présentent des difficultés numériques dans le traitement des régions complètement sèches (S=0) et complètement saturées (S=1)

Motivation de la recherche

Développer une méthode numérique unifiée capable de :

• Éviter la commutation de variables et la régularisation
• Traiter de manière stable les régions sèches et saturées
• Maintenir de bonnes propriétés de conservation de masse
• S'appliquer aux milieux poreux hétérogènes

Contributions principales

1. Proposition d'une nouvelle forme u : Introduction d'une variable auxiliaire bornée u qui élimine les termes non bornés de la formulation faible
2. Développement d'un schéma numérique stable : Combinaison de la discrétisation temporelle semi-implicite et de l'itération de Newton
3. Implémentation d'une méthode de décomposition de domaine : Traitement des problèmes de sols stratifiés et autres milieux hétérogènes
4. Validation de l'efficacité de la méthode : Démonstration de la stabilité et de la précision par plusieurs expériences numériques

Détails de la méthode

Formulation du problème

Résoudre l'équation de Richards : $\frac{\partial\theta}{\partial t} + \nabla \cdot q = 0$

où θ est la teneur en eau et q est le flux d'eau selon la loi de Darcy-Buckingham : $q = -K_s(x)K_r(x,S)\nabla(\Psi(x,S) + z)$

Construction de la nouvelle variable u

Pour tous les modèles, la dérivée de la fonction de Leverett peut s'écrire : $J'(S) = C(x)S^{-a}(1-S^c)^{-b}$

Définition de la variable auxiliaire : $u = \int_0^S (1-s^c)^{-b} ds$

Représentation utilisant la fonction bêta incomplète : $u = \frac{1}{c}B(S^c, \frac{1}{c}, 1-b)$

Équation de Richards sous la forme u

Dans un milieu homogène, la forme u est : $\phi \frac{\partial S(u)}{\partial t} - \nabla \cdot \left(K_sK_r(S(u))\left(h_{cap}CS(u)^{-a}\nabla u + e_z\right)\right) = 0$

Formulation variationnelle

Trouver u ∈ L²(I;H¹(Ω)) tel que : $\int_\Omega \phi \frac{\partial S(u)}{\partial t} v dx + \int_\Omega K_sK_r(S(u))\left(h_{cap}S(u)^{-a}\nabla u + e_z\right) \cdot \nabla v dx + \text{termes de frontière} = 0$

Discrétisation temporelle

Schéma d'Euler semi-implicite : $\int_\Omega \phi \frac{S(u_h^{n+1}) - S(u_h^n)}{\Delta t} v_h dx + \int_\Omega K_sK_r(S(u_h^n))\left(h_{cap}S(u_h^n)^{-a}\nabla u_h^{n+1} + e_z\right) \cdot \nabla v_h dx = 0$

Itération de Newton

Pour les termes non linéaires S(u^{n+1,m+1}), utilisation d'une expansion de Taylor au premier ordre : $S(u^{n+1,m+1}) = S(u^{n+1,m}) + \frac{\partial S}{\partial u}(u^{n+1,m})(u^{n+1,m+1} - u^{n+1,m}) + O((u^{n+1,m+1} - u^{n+1,m})^2)$

Méthode de décomposition de domaine

Pour les sols stratifiés, utilisation de la méthode de Schwarz non-chevauchante :

• Résolution indépendante de l'équation de Richards dans chaque sous-domaine
• Couplage par conditions d'interface : continuité de la charge de pression et conservation du flux
• Résolution itérative utilisant des conditions de transmission de type Robin

Configuration expérimentale

Implémentation numérique

• Plateforme logicielle : FreeFEM++ et FENICSx
• Solveurs linéaires : UMFPACK, PETSc, MUMPS
• Maillages : Maillages triangulaires uniformes, éléments finis linéaires (k=1)

Indicateurs d'évaluation

Erreur L² : $E_{h,\Delta t}(T)^2 = \|u(\cdot,T) - u_{h,\Delta t}(\cdot,T)\|_{L^2(\Omega)}^2$

Ordres de convergence :

• Ordre de convergence spatiale : $p = \log_2(E_{h,\Delta t}(T)/E_{h/2,\Delta t}(T))$
• Ordre de convergence temporelle : $q = \log_2(E_{h,\Delta t}(T)/E_{h,\Delta t/2}(T))$

Cas de test

1. Fibre 1D complètement non saturée : Vérification de la stabilité en conditions sèches
2. Région mixte 1D saturée-non saturée : Test de la capacité à traiter les valeurs initiales discontinues
3. Vérification par solution fabriquée : Calcul des ordres de convergence
4. Écoulement 2D complètement saturé : Comparaison avec le logiciel Hydrus
5. Problème de sol stratifié : Validation de la méthode de décomposition de domaine
6. Inclusion hétérogène dans un milieu : Test du traitement de géométries complexes

Résultats expérimentaux

Résultats principaux

Test de fibre 1D

• Maillage : 500 points de maillage uniformes
• Pas de temps : Δt = 1s
• Résultats : Les résultats de la forme u et de la forme S sont complètement identiques et en bon accord avec les données expérimentales

Région mixte saturée-non saturée

• Maillage : 5000 points de maillage uniformes
• Pas de temps : Δt = 10⁻⁵ heures
• Convergence de Newton : Moyenne de 1-5 itérations
• Stabilité : La solution conserve la monotonie, sans oscillations numériques

Analyse de convergence par solution fabriquée

Convergence spatiale (Δt = 5×10⁻⁵ fixé) :

Convergence temporelle (h = 1×10⁻³ fixé) :

Comparaison d'écoulement 2D saturé

• Maillage : 95,142 éléments triangulaires, 47,972 sommets
• Pas de temps : Δt = 0.001 heures
• Résultats : Accord élevé avec les résultats du logiciel Hydrus

Test de sol stratifié

• Maillage : 89,306 éléments triangulaires, 45,064 sommets
• Pas de temps : Δt = 0.05 heures
• Convergence : Tolérance 10⁻³, maximum 10 itérations, λ = 25
• Résultats : Accord avec les résultats de la littérature 20

Comparaison des méthodes

Dans le test d'inclusion hétérogène :

• Forme u : Convergence stable même lorsque ε=0, moyenne de 4 itérations
• Forme Ψ : Nécessite >1000 itérations lorsque ε=10⁻⁵, résidu de tolérance toujours 0.04

Travaux connexes

Principales directions de recherche

1. Méthodes de commutation de variables : Diersch & Perrochet (1999), Forsyth et al. (1995)
2. Techniques de régularisation : Schweizer (2007), Pop & Schweizer (2011)
3. Formes mixtes : Celia et al. (1990), Maina & Ackerer (2017)
4. Décomposition de domaine : Lions (1990), Manzini & Ferraris (2004)

Avantages de cet article

1. Traitement unifié : Pas besoin de commutation de variables ou de régularisation
2. Stabilité numérique : Stable sur toute la plage S ∈ 0,1
3. Cohérence physique : Maintient la continuité de la charge de pression
4. Efficacité de calcul : Convergence rapide de Newton

Conclusions et discussion

Conclusions principales

1. Validité de la forme u : Élimination réussie des termes non bornés des formes traditionnelles
2. Stabilité numérique : Maintien de la stabilité dans les régions complètement sèches et saturées
3. Performance de convergence : Convergence d'ordre deux en espace, d'ordre un en temps
4. Praticité : Applicable aux problèmes complexes de milieux hétérogènes

Limitations

1. Dépendance du modèle : Nécessite que $\lim_{S\to 0} K(S)S^{-a} < \infty$
2. Résolution non linéaire : Nécessite l'itération de Newton dans la région saturée
3. Traitement d'interface : Nécessite une technique de décomposition de domaine pour les milieux hétérogènes
4. Sensibilité aux paramètres : Le paramètre Robin λ dans la décomposition de domaine nécessite un ajustement

Directions futures

1. Analyse théorique : Preuve rigoureuse d'analyse numérique et de convergence
2. Couplage multi-physique : Extension aux transports de soluté et de chaleur
3. Techniques adaptatives : Développement de maillages adaptatifs et de pas de temps adaptatifs
4. Calcul parallèle : Implémentation parallèle à grande échelle

Évaluation approfondie

Points forts

1. Innovation forte : Construction ingénieuse de la variable u avec fondement théorique solide
2. Stabilité numérique : Évite les difficultés numériques des méthodes traditionnelles
3. Expériences suffisantes : Couvre des cas allant de 1D à 2D, d'homogène à hétérogène
4. Utilité pratique : La comparaison avec le logiciel commercial Hydrus renforce la crédibilité

Insuffisances

1. Analyse théorique insuffisante : Manque de preuves rigoureuses de convergence et de stabilité
2. Complexité de calcul : L'analyse du coût de calcul de l'itération de Newton n'est pas assez détaillée
3. Choix de paramètres : Le choix du paramètre Robin λ manque de guidance théorique
4. Vérification 3D : Absence de vérification sur des problèmes tridimensionnels

Impact

1. Valeur académique : Offre une nouvelle approche pour la résolution numérique de l'équation de Richards
2. Valeur pratique : Perspectives d'application largement applicables en hydrogéologie et science des sols
3. Reproductibilité : Fournit des détails d'implémentation détaillés, facilitant la reproduction

Scénarios d'application

1. Écoulement dans les milieux poreux : Écoulement des eaux souterraines, transport de l'eau dans les sols
2. Applications d'ingénierie : Infiltration à travers les barrages, analyse de stabilité des pentes
3. Problèmes environnementaux : Migration des polluants, restauration des eaux souterraines
4. Science des matériaux : Infiltration de fluides dans les matériaux fibreux

Références

Références clés

1. Richards, L.A. (1931). Capillary conduction of liquids through porous mediums. Physics 1(5), 318-333.
2. Celia, M.A., et al. (1990). A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26(7), 1483-1496.
3. van Genuchten, M.T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44(5), 892-898.
4. Manzini, G., Ferraris, S. (2004). Mass-conservative finite volume methods on 2-d unstructured grids for the Richards' equation. Adv. Water Resour. 27, 1199-1215.

Évaluation globale : Cet article constitue une contribution innovante importante dans le domaine de la résolution numérique de l'équation de Richards. En introduisant une variable auxiliaire bornée u, il résout ingénieusement les difficultés numériques des méthodes traditionnelles dans le traitement des régions sèches et saturées. La méthode possède une base théorique solide et une valeur pratique considérable, et les expériences numériques valident pleinement l'efficacité de la méthode. Bien que l'analyse théorique mérite d'être renforcée, il s'agit globalement d'un travail de recherche de haute qualité.