<figure><img src="image_1.png" /></figure>

HAL
open science

# Construction and analysis of numerical schemes for models from physics

Marianne Bessemoulin-Chatard

► To cite this version:

Marianne Bessemoulin-Chatard. Construction and analysis of numerical schemes for models from physics. Numerical Analysis [math.NA]. Nantes Université, 2023. tel-04187120

HAL Id: tel-04187120

https://hal.science/tel-04187120

Submitted on 24 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Laboratoire de Mathématiques Jean Leray
UMR CNRS 6629 et Université de Nantes

Habilitation à diriger des recherches
en Mathématiques Appliquées

présentée par

Marianne BESSEMOULIN-CHATARD

Construction et analyse de schémas
numériques pour des modèles issus de la
physique

Soutenue publiquement le 12 juillet 2023

Après avis de :

<table><tr><td>Boris</td><td>ANDREIANOV</td><td>Université de Tours</td></tr><tr><td>Pauline</td><td>LAFITTE</td><td>CentraleSupélec</td></tr><tr><td>Lorenzo</td><td>PARESCHI</td><td>University of Ferrara</td></tr></table>

Devant le jury composé de :

<table><tr><td>Boris</td><td>ANDREIANOV</td><td>Université de Tours</td></tr><tr><td>Christophe</td><td>BERTHON</td><td>Nantes Université</td></tr><tr><td>Clément</td><td>CANCÈS</td><td>Centre Inria de l'Université de Lille</td></tr><tr><td>Martin</td><td>GANDER</td><td>Université de Genève</td></tr><tr><td>Raphaèle</td><td>HERBIN</td><td>Aix-Marseille Université</td></tr><tr><td>Pauline</td><td>LAFITTE</td><td>CentraleSupélec</td></tr><tr><td>Lorenzo</td><td>PARESCHI</td><td>University of Ferrara</td></tr><tr><td>Marie-Hélène</td><td>VIGNAL</td><td>Université Toulouse 3 - Paul Sabatier</td></tr></table>
<figure><img src="image_3.png" /></figure>
Remerciements

Pour commencer, je tiens à adresser des remerciements tout particuliers à ma directrice et mon directeur de thèse, Claire Chainais-Hillairet et Francis Filbet. Ils ont guidé mes premiers pas dans la recherche, et j'ai pu bénéficier jusqu'à aujourd'hui de leurs conseils et de leur soutien. Je prends toujours beaucoup de plaisir à échanger avec eux, et à collaborer encore aujourd'hui sur des sujets passionnants.

Je remercie chaleureusement Boris Andreianov, Pauline Lafitte et Lorenzo Pareschi, qui m'ont fait l'honneur d'accepter d'être les rapporteurs de ce manuscrit. Leurs remarques judicieuses ainsi que leur avis ont beaucoup d'importance pour moi.

J'adresse également mes plus sincères remerciements à Christophe Berthon, Clément Cancès, Martin Gander, Raphaèle Herbin et Marie-Hélène Vignal qui ont accepté de faire partie de mon jury.

Je tiens aussi à remercier toutes les personnes avec qui j'ai eu la chance de discuter et de travailler depuis mes débuts dans le monde de la recherche. Le travail en collaboration a une très grande importance à mes yeux, c'est un des aspects que j'apprécie le plus dans mon métier de chercheuse. J'adresse en particulier mes remerciements à Maxime Herda, Giulia Lissoni, Hélène Mathis, Thomas Rey et Mazen Saad.

Je me dois également de mentionner ici et de remercier très chaleureusement les mem-
bres du LMJL, et en particulier de l'équipe MACS. Depuis mon arrivée à Nantes après ma
thèse, j'ai pu trouver une ambiance et un cadre de travail propices à mon développement
scientifique et personnel. Un grand merci également à toute l'équipe administrative et
technique qui se démène au quotidien pour nous simplifier la vie malgré des contraintes
de plus en plus complexes.

Je remercie aussi ici toutes les personnes qui font partie de ma vie non mathématique mais ont un impact important de par le soutien et la richesse qu'elles m'apportent au quotidien. Je pense notamment à la team EPSL, notamment les Reines.

Enfin, je termine par les incontournables, qui m'ont accordé un soutien indéfectible au cours de toutes ces années : ma famille. Sans eux, leur patience et leur confiance, rien ne serait possible. Mention spéciale à Stéphane, Léo et Louis, qui sont des sources de bonheur et de joie au quotidien.
<figure><img src="image_4.png" /></figure>
Introduction vii

1 Présentation des travaux vii

2 Liste des publications viii

I Schémas volumes finis pour des modèles de semi-conducteurs 1

1 Comportement en temps long d'une approximation du modèle de dérive-diffusion 3

1.1 Le modèle de dérive-diffusion 4

1.2 Un schéma de Scharfetter-Gummel généralisé 7

1.3 Inégalité d'entropie-production discrète 10

1.4 Existence et bornes $L^\infty$ uniformes 10

1.5 Convergence exponentielle vers l'équilibre thermique 14

1.6 Perspectives 15

2 Schémas numériques pour le modèle de transport d'énergie 19

2.1 Structure entropique discrète 21

2.2 Analyse numérique du schéma 25

2.3 Extension au cadre DDFV 28

2.4 Perspectives 29

II Schémas préservant la limite diffusive de systèmes hyperboliques avec termes sources 31

3 Taux de convergence vers la limite diffusive pour une discrétisation du p-système 33

3.1 Le cadre semi-discret en espace 36
# 1 Présentation des travaux

Je présente dans ce manuscrit une grande partie des travaux effectués depuis ma thèse. Mes recherches concernent le développement et l'analyse de schémas numériques de type volumes finis discrétisant des équations aux dérivées partielles, qui apparaissent notamment dans des modèles issus de la physique et de la biologie. Je m'intéresse plus particulièrement à la préservation au niveau discret de propriétés du modèle continu : comportements asymptotiques (en temps long, à la limite de diffusion) et conservation de quantités physiques entre autres. Dans tous les travaux présentés ici, une attention particulière est accordée à la discrétisation, qui est construite de sorte à pouvoir adapter au cadre discret les outils d'analyse du cadre continu, et ainsi proposer une étude mathématique rigoureuse des méthodes numériques obtenues.

La première partie de ce manuscrit concerne la construction et l'étude de schémas numériques pour des modèles de semi-conducteurs, et s'inscrit dans la continuité de mes travaux de thèse. Dans la vaste hiérarchie de modèles décrivant les semi-conducteurs [Jün11], je m'intéresse à l'approximation de deux modèles macroscopiques : le système de dérive-diffusion et le modèle de transport d'énergie.

Dans le chapitre 1, j'étudie le comportement en temps long d'une discrétisation du système de dérive-diffusion. Le schéma considéré est de type volumes finis, avec flux de Scharfetter-Gummel généralisés, permettant la prise en compte à la fois de lois de pression linéaires et non linéaires. La convergence en temps long des solutions numériques vers une approximation de l'équilibre thermique à taux exponentiel est démontrée en contrôlant l'entropie relative par la production d'entropie. Ce résultat est obtenu sous l'hypothèse de bornes $L^∞$ uniformes en temps sur les densités approchées de porteurs de charge, que nous prouvons rigoureusement dans le cas où la diffusion est linéaire en adaptant au cadre discret une technique d'itération de Moser.

Dans le chapitre 2, je présente la construction et l'analyse de schémas volumes finis pour le modèle de transport d'énergie. L'idée est de construire un schéma de sorte à pouvoir en obtenir une reformulation équivalente en variables duales entropiques, permettant de démontrer la décroissance d'une entropie discrète et le contrôle de la production d'en-
tropie discrète correspondante. À partir de là, un certain nombre de propriétés peuvent être analysées : estimations *a priori*, existence d'une solution approchée, comportement en temps long. Nous considérons tout d'abord un schéma de type TPFA (*Two Points Flux Approximation*) puis une généralisation au cadre DDFV (*Discrete Duality Finite Volume*) permettant de considérer des maillages plus généraux.

La deuxième partie est dédiée à l'étude de la convergence vers la limite de diffusion de schémas numériques pour des systèmes hyperboliques avec termes sources.

Le chapitre 3 est consacré à l'étude de convergence à la limite de diffusion d'un schéma préservant l'asymptotique pour le *p*-système avec amortissement. La nouveauté principale de ce travail est l'établissement d'un taux explicite de convergence, obtenu par une estimation d'erreur entre les solutions approchées du système hyperbolique et la limite diffusive approchée grâce à une version discrète de la méthode d'entropie relative [LT13].

Le modèle considéré dans le chapitre 4, le système de Saint-Venant avec friction de Manning, fait intervenir une difficulté supplémentaire. En effet, le terme source étant quadra-tique, l'équation de diffusion limite implique un opérateur non linéaire de type *p*-laplacien.

Nous proposons deux discrétisations du problème hyperbolique, de type Godunov à deux états intermédiaires, asymptotiquement consistantes avec la limite de diffusion. Cette propriété est établie de manière formelle à l'aide de développements de Chapman-Enskog de la solution approchée.

Dans la troisième partie, je me concentre sur l'analyse numérique de schémas pour des équations cinétiques unidimensionnelles relativement simples.

Dans le chapitre 5, je m'intéresse à la discrétisation d'équations cinétiques linéaires sans champ électromagnétique, où l'opérateur de collision est de type Fokker-Planck ou BGK linéaire. Le schéma volumes finis proposé permet d'établir des estimations *a priori* dont nous déduisons d'une part la préservation de la limite de diffusion au niveau discret, et d'autre part la convergence en temps long, à taux exponentiel, vers l'équilibre. Il s'agit d'une adaptation au cadre discret de la méthode d'hypocoercivité $L^2$ [DMS15].

Enfin, le chapitre 6 est consacré à l'étude de stabilité et de convergence de discrétisations conservatives du système de Vlasov-Poisson écrit sous forme d'un système hyperbolique à l'aide de polynômes de Hermite. L'idée principale est de se placer dans un cadre fonctionnel adapté aux variations de la fonction de distribution au cours du temps. Nous introduisons pour cela un espace $L^2$ à poids dépendant du temps, et étudions les propriétés du schéma obtenu en considérant une discrétisation en espace de type Galerkin discontinue.

## 2 Liste des publications

### 2.1 Travaux présentés dans le mémoire

#### Schémas volumes finis pour des modèles de semi-conducteurs

* (Avec C. Chainais-Hillairet) *Exponential decay of a finite volume scheme to the thermal equilibrium for drift-diffusion systems*, J. Numer. Math., **25** (2017), no. 3,
147–168.

• (Avec C. Chainais-Hillairet et A. Jüngel) *Uniform L∞ Estimates for Approximate Solutions of the Bipolar Drift-Diffusion System*, In International Conference on Finite Volumes for Complex Applications, Springer (2017), 381–389.

• (Avec C. Chainais-Hillairet) *Uniform-in-time bounds for approximate solutions of the drift-diffusion system*, Numerische Mathematik, **141** (2019), no. 4, 881–916.

• (Avec C. Chainais-Hillairet et H. Mathis) *Numerical schemes for semiconductors energy-transport models*, In Finite volumes for complex applications IX, Springer, Cham (2020), 75–90.

• (Avec C. Chainais-Hillairet et H. Mathis) *Analysis of numerical schemes for semiconductor energy-transport models*, Appl. Numer. Math., **168** (2021), 143–169.

• (Avec G. Lissoni et H. Mathis) *Numerical analysis of DDFV schemes for semiconductors energy-transport models*, Comput. Appl. Math., **41** (2022), no. 15.

## Schémas préservant la limite diffusive de systèmes hyperboliques avec termes sources

• (Avec C. Berthon et H. Mathis) *Numerical convergence rate for a diffusive limit of hyperbolic systems: p-system with damping*, SMAI J. Comput. Math., **2** (2016), 99–119.

• (Avec C. Berthon et S. Bulteau) *Convergence rate of an Asymptotic Preserving scheme for the diffusive limit of the p-system with damping*, Communications in Mathematical Sciences, **17** (2019), no. 6, 1459–1486.

• (Avec C. Berthon et S. Bulteau) *An Asymptotic Preserving scheme for the shallow-water equations with Manning friction using viscous correction of the HLL scheme* (2019), soumis.

• (Avec M. Badsi, C. Berthon et S. Bulteau) *A fully well-balanced and asymptotic preserving scheme for the shallow-water equations with a generalized Manning friction source term*, Calcolo, **58** (2021), no. 41.

## Analyse numérique pour des problèmes issus de la théorie cinétique

• (Avec M. Herda et T. Rey) *Hypocoercivity and diffusion limit of a finite volume scheme for linear kinetic equations*, Math. Comp., **89** (2020), no. 323, 1093–1133.

• (Avec F. Filbet) *On the stability of conservative discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system*, J. Comput. Phys., **451** (2022), no. 28.

• (Avec F. Filbet) *On the convergence of discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system* (2022), soumis.
## 2.2 Autres travaux non présentés dans le mémoire

• (Avec M. Saad) *Monotone Combined Finite Volume-Finite Element Scheme for a Bone Healing Model*, In Finite Volumes for Complex Applications VII, Springer (2014), 497-505.

• (Avec M. Saad) *Preserving monotony of combined edge finite volume-finite element scheme for a bone healing model on general mesh*, J. Comput. Appl. Math., **309** (2017), 287-311.

• (Avec C. Berthon, A. Crestetto et F. Foucher) *A Riemann solution approximation based on the zero diffusion-dispersion limit of Dafermos reformulation type problem*, Calcolo, **56** (2019), no. 3.

• (Avec E. H. Quenjel, M. Saad et M. Ghilani) *Convergence of a positive non-linear DDFV scheme for degenerate parabolic equations*, Calcolo, **57** (2020), no. 2.

## 2.3 Travaux de thèse

• (Avec C. Chainais-Hillairet et F. Filbet) *On discrete functional inequalities for some finite volume schemes*, IMA Journal of Numerical Analysis, **35** (2015), no. 3, 1125-1149.

• (Avec C. Chainais-Hillairet et M.-H. Vignal) *Study of a finite volume scheme for the drift-diffusion system. Asymptotic behavior in the quasi-neutral limit*, SIAM Journal on Numerical Analysis, **52** (2014), no. 4, 1666-1691.

• (Avec A. Jüngel) *A finite volume scheme for a Keller-Segel model with additional cross-diffusion*, IMA Journal of Numerical Analysis, **34** (2014), no. 1, 96-122.

• (Avec F. Filbet) *A finite volume scheme for nonlinear degenerate parabolic equations*, SIAM Journal on Scientific Computing, **34** (2012), no. 5, B559-B583.

• *A finite volume scheme for convection-diffusion equations with nonlinear diffusion derived from the Scharfetter-Gummel scheme*, Numerische Mathematik, **121** (2012), no. 4, 637-670.

• *Asymptotic behavior of the Scharfetter-Gummel scheme for the drift-diffusion model*, In Finite Volumes for Complex Applications 6, Springer Verlag (2011).
# PREMIÈRE PARTIE

## SCHÉMAS VOLUMES FINIS POUR DES MODÈLES DE SEMI-CONDUCTEURS
<figure><img src="image_5.png" /></figure>
# CHAPITRE 1

## Comportement en temps long d'une approximation du modèle de dérive-diffusion

Dans ce chapitre, je présente des résultats concernant le comportement en temps long d'une discrétisation du système de dérive-diffusion pour les semi-conducteurs. Dans le cadre continu, la preuve de la convergence en temps long à taux exponentiel vers l'équilibre thermique a notamment été obtenue dans [GG86, GG89]. La technique utilisée s'inscrit dans le cadre des méthodes d'entropie-production [ACD$^{+}$04]. Ces méthodes consistent à étudier une fonction de Liapounov $\mathbb{E}$, appelée entropie, et sa production $\mathbb{I}$ qui est positive. Ces deux quantités sont reliées entre elles par une inégalité du type

$$ (1.1) \qquad 0 \le \mathbb{E}(t) + \int_0^t \mathbb{I}(s) \, ds \le \mathbb{E}(0), $$

qui permet de prouver la convergence en entropie relative des solutions du système de dérive-diffusion vers l'équilibre thermique. Si de plus l'entropie relative est contrôlée par la production d'entropie dans une relation de la forme

$$ (1.2) \qquad \mathbb{E}(t) \le C\mathbb{I}(t), $$

un taux de convergence exponentiel est obtenu. Cette méthode générique peut s'appliquer dans de nombreux contextes, par exemple à des équations de type Fokker-Planck [AMTU01] ou à des équations paraboliques dégénérées [CJM$^{+}$01].

Dans le cas du système de dérive-diffusion, la fonction de Liapounov considérée est l'énergie électrochimique. L'étude de cette quantité permet de prouver la convergence exponentielle vers l'équilibre thermique. Mon objectif est d'adapter la méthode d'entropie-production au niveau discret pour étudier le comportement en temps long des solutions approchées obtenues à partir d'un schéma numérique pour le modèle de dérive-diffusion.
Les travaux présentés dans ce chapitre sont dans la continuité naturelle de mes travaux de thèse. En effet, j'avais d'abord étudié le comportement en temps long du schéma de Scharfetter-Gummel pour le modèle de dérive-diffusion isotherme, et obtenu la convergence en temps long de la solution numérique vers une approximation de l'équilibre thermique, en établissant un analogue discret de (1.1). Cependant, cette première contribution [16] comporte deux limitations importantes. D'une part, en l'absence d'une version discrète de (1.2), le caractère exponentiel de la convergence n'est pas prouvé pour la solution numérique. D'autre part, le schéma de Scharfetter-Gummel classique ne s'applique que dans le cas de la statistique de Boltzmann, exploitant la dépendance exponentielle des densités de porteurs de charge en les potentiels chimiques, et correspondant à une diffusion linéaire. J'ai proposé au cours de ma thèse une généralisation du schéma de Scharfetter-Gummel utilisant une moyenne appropriée de la diffusion non linéaire, qui garantit la consistance thermodynamique [2]. L'idée était donc d'étendre les résultats obtenus dans le cas d'une diffusion linéaire [16] au cas non linéaire plus général en utilisant ce schéma généralisé.

Dans l'article [3] écrit en collaboration avec Claire Chainais-Hillairet, nous abordons ces deux questions. Plus précisément, nous démontrons des analogues discrets des inégalités (1.1) et (1.2), et obtenons ainsi la convergence exponentielle (Théorème 1.3), à la fois dans le cas isotherme (linéaire), mais également dans le cas non linéaire.

Ce résultat est établi sous l'hypothèse assez forte que les densités discrètes vérifient des bornes supérieures et inférieures strictement positives uniformes en temps. Nous avons par la suite obtenu ces bornes dans le cas isotherme [6, 4] (Théorème 1.1), en adaptant au cadre discret une technique d'itération de Moser développée dans le cadre continu dans [GG86, GG89].

Dans la section 1.1, j'introduis le modèle de dérive-diffusion et les hypothèses considérées. La section 1.2 est consacrée à la définition du schéma numérique. La section 1.3 est dédiée à l'analogue discret de l'inégalité d'entropie-production (1.1), qui est fondamental pour établir les résultats présentés dans la suite. Dans les sections 1.4 et 1.5, je présente les principaux résultats que j'ai obtenus concernant l'existence de solutions numériques vérifiant des bornes uniformes en temps, et convergeant à taux exponentiel vers une approximation de l'équilibre thermique.

## 1.1 Le modèle de dérive-diffusion

Le modèle de dérive-diffusion est constitué de deux équations de continuité pour les densités d'électrons N et de trous P, et d'une équation de Poisson pour le potentiel électrique Ψ. Étant donné Ω un ouvert borné de $\mathbb{R}^d$ correspondant à la géométrie du semi-conducteur considéré, et $T > 0$ un temps final, le système s'écrit : pour tous $(t, x) \in (0, T) \times \Omega$,

$$ (1.3) \qquad \partial_t N + \text{div}(-\nabla r(N) + N \nabla \Psi) = -R(N, P), $$

$$ (1.4) \qquad \partial_t P + \text{div}(-\nabla r(P) - P \nabla \Psi) = -R(N, P), $$

$$ (1.5) \qquad -\lambda^2 \Delta \Psi = P - N + C. $$
La fonction donnée $C(x)$ est le profil de dopage, et $\lambda$ est la longueur de Debye adimensionnée. La définition de la fonction $r$ apparaissant dans le terme de diffusion dépend de la statistique choisie pour décrire la relation entre les densités de porteurs de charge et les potentiels chimiques. Je considère dans tout ce chapitre une fonction $r$ de la forme [Jü94, Jü95]

$$r(s) = s^{\alpha}, \quad \alpha \geq 1.$$

Le cas linéaire ($\alpha = 1$) est le modèle isotherme, correspondant au choix de la statistique de Boltzmann. Dans le cas non linéaire, ce choix de fonction puissance correspond à une approximation de la distribution de Fermi-Dirac dans la limite de haute densité.

Le terme $R(N, P)$ est le taux de recombinaison-génération. Dans le cas isotherme ($r = Id$), ce terme peut être supposé de la forme

$$(1.6) \qquad R(N, P) = R_0(N, P)(NP - 1),$$

cette expression générale incluant en particulier les termes de Shockley-Read-Hall et Auger [MRS90].

Le système (1.3)-(1.5) est complété avec des conditions initiales $N_0, P_0$, et des conditions au bord mixtes de type Dirichlet Neumann, supposées indépendantes du temps :

$$(1.7) \qquad N = N^D, \quad P = P^D, \quad \Psi = \Psi^D \quad \text{sur } \Gamma^D,$$

$$(1.8) \qquad \nabla r(N) \cdot \nu = \nabla r(P) \cdot \nu = \nabla \Psi \cdot \nu = 0 \quad \text{sur } \Gamma^N,$$

où $\partial\Omega = \Gamma^D \cup \Gamma^N$ et $\nu$ est le vecteur unitaire sortant normal à $\partial\Omega$.

Tous les résultats présentés dans ce chapitre sont valables sous les hypothèses suivantes.

**Hypothèses H.1.** • Le domaine $\Omega$ est un ouvert borné de $\mathbb{R}^d$ ($d \geq 1$) et $\partial\Omega = \Gamma^D \cup \Gamma^N$, avec $\Gamma^D \cap \Gamma^N = \emptyset$ et $m(\Gamma^D) > 0$.

• Le dopage $C$ est dans $L^\infty(\Omega)$.

• Les conditions au bord $N^D, P^D$ et $\Psi^D$ sont les traces de fonctions définies sur $\Omega$ tout entier, encore notées $N^D, P^D$ et $\Psi^D$, et telles que

$$(1.9) \qquad N^D, P^D \in L^\infty(\Omega) \cap H^1(\Omega), \Psi^D \in H^1(\Omega).$$

• Les conditions initiales $N_0, P_0$ sont supposées dans $L^\infty(\Omega)$.

• Il existe $M, m > 0$ telles que

$$(1.10) \qquad m \leq N_0, P_0, N^D, P^D \leq M \text{ p.p. dans } \Omega.$$

Nous allons par la suite considérer deux cas, correspondant chacun à l'un des deux jeux d'hypothèses suivants : le cas isotherme avec terme de recombinaison-génération (Hypothèses H.2) et le cas non linéaire sans terme de recombinaison génération (Hypothèses H.3).
**Hypothèses H.2.** • La fonction r est l'identité.

* Le terme de recombinaison-génération est de la forme (1.6), avec $R_0$ continue, positive, vérifiant

$$ (1.11) \qquad \exists \bar{R} > 0 \text{ tel que } 0 \le R_0(N, P) \le \bar{R} (1 + |N| + |P|) \quad \forall N, P \in \mathbb{R}. $$

* Les conditions au bord satisfont $N^D P^D = 1$.

**Hypothèses H.3.** • La fonction r est dans $C^1(\mathbb{R})$, et satisfait

$$ r(0) = r'(0) = 0, \quad r'(s) \ge c_0 s^{\alpha-1}, \quad \alpha > 1. $$

* Le taux de recombinaison-génération R est supposé nul.

Je rappelle pour finir la définition de l'équilibre thermique. Il s'agit d'un état stationnaire particulier pour lequel les courants d'électrons et de trous s'annulent :

$$ -\nabla r(N) + N \nabla \Psi = -\nabla r(P) - P \nabla \Psi = 0. $$

Pour définir cet équilibre, introduisons la fonction enthalpie

$$ h(s) = \int_{1}^{s} \frac{r'(\tau)}{\tau} d\tau, $$

et l'inverse généralisée de $h$

$$ g(s) = \begin{cases} h^{-1}(s) & \text{if } h(0^{+}) < s < \infty, \\ 0 & \text{if } s \le h(0^{+}), \end{cases} $$

en supposant que $h(+\infty) = +\infty$. Le cas isotherme correspond simplement à $h = \log$ et $g = \exp$.

Si les conditions de Dirichlet vérifient la condition de compatibilité

$$ (1.12) \qquad h(N^D) - \Psi^D = \alpha_N \text{ et } h(P^D) + \Psi^D = \alpha_P \text{ sur } \Gamma^D, $$

alors l'équilibre thermique est défini dans $\Omega$ par

$$ (1.13) \qquad -\lambda^2 \Delta \Psi^{eq} = g(\alpha_P - \Psi^{eq}) - g(\alpha_N + \Psi^{eq}) + C, $$

$$ (1.14) \qquad N^{eq} = g(\alpha_N + \Psi^{eq}), $$

$$ (1.15) \qquad P^{eq} = g(\alpha_P - \Psi^{eq}), $$

avec les conditions au bord (1.7)–(1.8). L'existence et l'unicité de cet équilibre sont démontrées dans [MRS90] pour le cas linéaire, et dans [MU93] pour le cas non linéaire.

Comme rappelé dans l'introduction de ce chapitre, la preuve de convergence en temps long vers l'équilibre thermique est fondée sur une méthode d'entropie-production. Ici la fonctionnelle de Liapounov considérée est l'énergie relative, déviation de l'énergie totale
(somme des énergies internes pour les densités d'électrons et de trous et de l'énergie due au potentiel électrique) par rapport à l'équilibre thermique :

$$ (1.16) \quad \mathbb{E}(t) = \int_{\Omega} \left( H(N(t)) - H(N^{eq}) - h(N^{eq})(N(t) - N^{eq}) \right. \\ \left. \qquad + H(P(t)) - H(P^{eq}) - h(P^{eq})(P(t) - P^{eq}) \right. \\ \left. \qquad + \frac{\lambda^2}{2} |\nabla(\Psi(t) - \Psi^{eq})|^2 \right) d\mathbf{x}, $$

avec $H$ une primitive de $h$, et la production d'entropie correspondante est donnée par

$$ (1.17) \quad \mathbb{I}(t) = \int_{\Omega} \left( N|\nabla(h(N) - \Psi)|^2 + P|\nabla(h(P) + \Psi)|^2 \right) d\mathbf{x} \\ + \int_{\Omega} R(N, P) (h(N) + h(P) - h(N^{eq}) - h(P^{eq})) d\mathbf{x}. $$

## 1.2 Un schéma de Scharfetter-Gummel généralisé

Je présente maintenant les schémas numériques considérés pour le système de dérive-diffusion (1.3)-(1.8) et pour l'équilibre thermique (1.14)-(1.13). Au niveau numérique, il est essentiel de considérer des méthodes qui préservent les principales propriétés qualitatives du système continu, telles que par exemple la positivité des densités et la convergence vers l'équilibre thermique. Dans le cas d'une diffusion linéaire, le schéma de Scharfetter-Gummel [Il'69, SG69] est particulièrement bien adapté. Dans le cas non linéaire, plusieurs généralisations de ce schéma ont été proposées [Jün95, PCC89], mais ne s'avèrent pas appropriées pour étudier le comportement en temps long. Plus récemment, une autre généralisation a été construite dans l'esprit du schéma de Scharfetter-Gummel original [EFG06], mais apparaît compliquée à mettre en œuvre, car conduit à résoudre un problème non linéaire à chaque interface pour calculer les flux numériques. Nous avons donc plutôt considéré la généralisation que j'ai proposée dans un de mes travaux de thèse [2], dont une interprétation alternative ainsi qu'une extension à des modèles de semi-conducteurs organiques ont été développées dans [KRF+15].

Je me place dans le cadre d'une discrétisation de type volume fini classique en espace [EGH00]. Le maillage $\mathcal{M} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$ du domaine $\Omega$ est donné par une famille $\mathcal{T}$ de volumes de contrôle, une famille $\mathcal{E}$ d'arêtes (ou faces) et une famille $\mathcal{P}$ de points $(x_K)_{K \in \mathcal{T}}$. Le maillage est supposé admissible au sens de [EGH00, Definition 9.1]. Dans l'ensemble des arêtes, nous distinguons les arêtes intérieures $\sigma = K|L \in \mathcal{E}_{int}$ des arêtes de bord $\sigma \in \mathcal{E}_{ext}$. Parmi ces arêtes de bord, nous distinguons également le bord Dirichlet $\Gamma^D$ du bord Neumann $\Gamma^N : \mathcal{E}_{ext} = \mathcal{E}_{ext}^D \cup \mathcal{E}_{ext}^N$. Pour un volume de contrôle $K \in \mathcal{T}$, $\mathcal{E}_K$ est l'ensemble de ses arêtes, qui est également scindé en $\mathcal{E}_K = \mathcal{E}_{K,int} \cup \mathcal{E}_{K,ext}^D \cup \mathcal{E}_{K,ext}^N$. Pour toute arête $\sigma \in \mathcal{E}$, nous définissons $d_\sigma = d(x_K, x_L)$ si $\sigma = K|L \in \mathcal{E}_{int}$ et $d_\sigma = d(x_K, \sigma)$ si $\sigma \in \mathcal{E}_{K,ext}$. Le coefficient de transmissibilité est alors défini par $\tau_\sigma = m(\sigma)/d_\sigma$.
Je suppose dans tout ce chapitre que le maillage vérifie les conditions de régularité et de non dégênerescence suivantes :

$$
(1.18) \qquad \exists \xi > 0 \text{ tel que } d(x_K, \sigma) \ge \xi d_\sigma, \quad \forall K \in \mathcal{T}, \forall \sigma \in \mathcal{E}_K,
$$

$$
(1.19) \qquad \exists \tau_m > 0 \text{ tel que } \tau_\sigma \ge \tau_m, \quad \forall \sigma \in \mathcal{E}.
$$

Concernant la discrétisation en temps, soit $0 < \Delta t \le 1$ le pas de temps. Nous posons
$N_T = E(T/\Delta t)$ et $t^n = n\Delta t$ pour tout $0 \le n \le N_T$. La taille de la discrétisation espace-
temps est définie par

$$
(1.20) \qquad \delta = \max \left( \Delta t, \max_{K \in \mathcal{T}} \operatorname{diam}(K) \right).
$$

Un schéma volumes finis approchant une équation d'inconnue $u$ fournit un vecteur $u_{\mathcal{T}} = (u_K)_{K \in \mathcal{T}} \in \mathbb{R}^{\theta}$ (avec $\theta = \text{Card}(\mathcal{T})$) qui est identifié avec la fonction constante par morceaux

$$
u_{\mathcal{T}} = \sum_{K \in \mathcal{T}} u_K \mathbf{1}_K.
$$

Puisque des conditions de Dirichlet sont imposées sur une partie du bord, nous avons
également besoin de définir $u_{\mathcal{E}^D} = (u_\sigma)_{\sigma \in \mathcal{E}_{\text{ext}}^D}$. Nous notons alors $X(\mathcal{M})$ l'ensemble des
fonctions discrètes $u_{\mathcal{M}} = (u_{\mathcal{T}}, u_{\mathcal{E}^D})$, et pour tout $u_{\mathcal{M}} \in X(\mathcal{M})$, nous définissons pour tout
$K \in \mathcal{T}$ et tout $\sigma \in \mathcal{E}_K$ :

$$
u_{K,\sigma} = \begin{cases} u_L & \text{si } \sigma = K \mid L \in \mathcal{E}_{K,\text{int}} \\ u_\sigma & \text{si } \sigma \in \mathcal{E}_{K,\text{ext}}^D \\ u_K & \text{si } \sigma \in \mathcal{E}_{K,\text{ext}}^N \end{cases}
$$

et

$$
D_{K,\sigma}u = u_{K,\sigma} - u_K, \quad D_{\sigma}u = |D_{K,\sigma}u|.
$$

Pour finir, j'introduis la norme $L^2$ et la semi-norme $H^1$ sur $X(\mathcal{M})$ : pour tout $u_{\mathcal{M}} \in X(\mathcal{M})$,

$$
|u_{\mathcal{M}}|_{1,\mathcal{M}}^2 = \sum_{\sigma \in \mathcal{E}} \tau_{\sigma} (D_{\sigma} u)^2,
$$

$$
\|u_{\mathcal{M}}\|_{0,\mathcal{M}}^2 = \|u_{\mathcal{T}}\|_{0}^2 = \sum_{K \in \mathcal{T}} m(K) u_K^2,
$$

où $\| \cdot \|_0$ est la norme $L^2$ usuelle pour les fonctions constantes par morceaux.

Le schéma pour le système de dérive-diffusion (1.3)–(1.8). Les conditions initiales et conditions au bord sont discrétisées classiquement par

$$
(1.21) \qquad \left(N_K^0, P_K^0\right) = \frac{1}{m(K)} \int_K \left(N_0(x), P_0(x)\right) \mathrm{d}x, \quad \forall K \in \mathcal{T},
$$

$$
(1.22) \qquad \left(N_{\sigma}^{D}, P_{\sigma}^{D}, \Psi_{\sigma}^{D}\right) = \frac{1}{m(\sigma)} \int_{\sigma} \left(N^{D}(\gamma), P^{D}(\gamma), \Psi^{D}(\gamma)\right) d\gamma, \quad \forall \sigma \in \mathcal{E}_{\text{ext}}^{D},
$$
et les valeurs au bord sont définies par

$$
(1.23) \qquad N_{\sigma}^{n} = N_{\sigma}^{D}, \quad P_{\sigma}^{n} = P_{\sigma}^{D}, \quad \Psi_{\sigma}^{n} = \Psi_{\sigma}^{D}, \quad \forall \sigma \in \mathcal{E}_{\text{ext}}^{D}, \quad \forall n \geq 0.
$$

Je considère un schéma complètement implicite en temps et volumes finis en espace :
$\forall K \in \mathcal{T}, \forall n \ge 0,$

$$
(1.24) \qquad \mathrm{m}(K) \frac{N_K^{n+1} - N_K^n}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \mathcal{F}_{K,\sigma}^{n+1} = -\mathrm{m}(K) R(N_K^{n+1}, P_K^{n+1}),
$$

$$
(1.25) \qquad \mathrm{m}(K) \frac{P_K^{n+1} - P_K^n}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \mathcal{G}_{K,\sigma}^{n+1} = -\mathrm{m}(K) R(N_K^{n+1}, P_K^{n+1}),
$$

$$
(1.26) \qquad -\lambda^2 \sum_{\sigma \in \mathcal{E}_K} \tau_\sigma D_{K,\sigma} \Psi^n = \mathrm{m}(K) (P_K^n - N_K^n + C_K).
$$

Les flux numériques sont de type Scharfetter-Gummel généralisés, comme introduits dans
un de mes travaux de thèse [2] :

$$
(1.27) \qquad \mathcal{F}_{K,\sigma}^{n+1} = \tau_\sigma dr(N_K^{n+1}, N_{K,\sigma}^{n+1}) \left[ \begin{aligned} & B \left( \frac{-D_{K,\sigma} \Psi^{n+1}}{dr(N_K^{n+1}, N_{K,\sigma}^{n+1})} \right) N_K^{n+1} \\ & - B \left( \frac{D_{K,\sigma} \Psi^{n+1}}{dr(N_K^{n+1}, N_{K,\sigma}^{n+1})} \right) N_{K,\sigma}^{n+1} \end{aligned} \right],
$$

$$
(1.28) \qquad \mathcal{G}_{K,\sigma}^{n+1} = \tau_\sigma dr(P_K^{n+1}, P_{K,\sigma}^{n+1}) \left[ \begin{aligned} & B \left( \frac{D_{K,\sigma} \Psi^{n+1}}{dr(P_K^{n+1}, P_{K,\sigma}^{n+1})} \right) P_K^{n+1} \\ & - B \left( \frac{-D_{K,\sigma} \Psi^{n+1}}{dr(P_K^{n+1}, P_{K,\sigma}^{n+1})} \right) P_{K,\sigma}^{n+1} \end{aligned} \right],
$$

avec

$$
(1.29) \quad dr(a,b) = \begin{cases} \frac{h(b) - h(a)}{\log(b) - \log(a)} & \text{si } a, b > 0, \ a \neq b, \\ r'\left(\frac{a+b}{2}\right) & \text{sinon,} \end{cases}
$$

$B$ étant la fonction de Bernoulli donnée par

$$
B(0) = 1 \text{ et } B(x) = \frac{x}{\exp(x) - 1} \quad \forall x \neq 0.
$$

**Schéma pour l'équilibre thermique (1.14)–(1.13).** L'approximation ($N_{\mathcal{T}}^{eq}$, $P_{\mathcal{T}}^{eq}$, $\Psi_{\mathcal{T}}^{eq}$) de l'équilibre thermique est obtenue par le schéma volumes finis suivant : $\forall K \in \mathcal{T}$,

$$
(1.30) \qquad -\lambda^2 \sum_{\sigma \in \mathcal{E}_K} \tau_\sigma D_{K,\sigma} \Psi^{eq} = m(K) (g(\alpha_P - \Psi_K^{eq}) - g(\alpha_N + \Psi_K^{eq}) + C_K),
$$

$$
(1.31) \qquad N_K^{eq} = g(\alpha_N + \Psi_K^{eq}),
$$

$$
(1.32) \qquad P_K^{eq} = g(\alpha_P - \Psi_K^{eq}).
$$
## 1.3 Inégalité d'entropie-production discrète

Je présente maintenant l'analogue discret de (1.1), qui est fondamental pour la suite. Pour cela, définissons la fonctionnelle d'entropie relative discrète

$$ (1.33) \qquad \begin{aligned} \mathbb{E}^n = \sum_{K \in \mathcal{T}} & \mathrm{m}(K) (H(N_K^n) - H(N_K^{eq}) - h(N_K^{eq})(N_K^n - N_K^{eq}) \\ & + H(P_K^n) - H(P_K^{eq}) - h(P_K^{eq})(P_K^n - P_K^{eq})) \\ & + \frac{\lambda^2}{2} |\Psi_{\mathcal{M}}^n - \Psi_{\mathcal{M}}^{eq}|_{1,\mathcal{M}}^2, \end{aligned} $$

et la production d'entropie discrète correspondante

$$ (1.34) \qquad \begin{aligned} \mathbb{I}^n = \sum_{\sigma \in \mathcal{E}} \tau_{\sigma} \left[ \min(N_K^n, N_{K,\sigma}^n) (D_{\sigma}(h(N^n) - \Psi^n))^2 \right. \\ \left. \qquad + \min(P_K^n, P_{K,\sigma}^n) (D_{\sigma}(h(P^n) + \Psi^n))^2 \right] \\ + \sum_{K \in \mathcal{T}} \mathrm{m}(K) R(N_K^n, P_K^n) [h(N_K^n) + h(P_K^n) - h(N_K^{eq}) - h(P_K^{eq})]. \end{aligned} $$

Grâce aux propriétés des flux de Scharfetter-Gummel (généralisés), la proposition suivante est établie.

**Proposition 1.1.** Supposons les Hypothèses H.1 et la condition de compatibilité (1.12) vérifiées. Supposons de plus que les Hypothèses H.2 ou H.3 sont satisfaites.

Soit $\mathcal{M} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$ un maillage admissible de $\Omega$. S'il existe une solution $(N_{\mathcal{T}}^n, P_{\mathcal{T}}^n, \Psi_{\mathcal{T}}^n)_{n \ge 0}$ au schéma (1.24)-(1.28), alors l'inégalité suivante est vérifiée pour $n \ge 0$ :

$$ (1.35) \qquad 0 \le \mathbb{E}^{n+1} + \Delta t \mathbb{I}^{n+1} \le \mathbb{E}^n. $$

## 1.4 Existence et bornes $L^\infty$ uniformes

Dans cette section, je présente les résultats obtenus concernant l'existence d'une solution au schéma pour le modèle de dérive-diffusion (1.24)-(1.28) et les bornes $L^\infty$ uniformes en temps obtenues sur les densités approchées. Je distingue deux cas.

* Le cas isotherme avec terme de recombinaison-génération (Hypothèses H.2), pour lequel un résultat complet est démontré, sans restriction sur le dopage C.

* Le cas non linéaire sans terme de recombinaison-génération (Hypothèses H.3), pour lequel un résultat complet est obtenu uniquement dans le cas d'un dopage nul. Dans le cas d'un dopage général, le résultat d'existence est conditionné à une restriction sur le pas de temps, et les bornes $L^\infty$ sur les densités approchées dépendent du temps et dégénèrent quand $t \to +\infty$.
### 1.4.1 Le cas isotherme avec terme de recombinaison-génération

**Théorème 1.1.** Sous les Hypothèses H.1 et H.2, supposons de plus la condition de compatibilité (1.12). Soit $\mathcal{M} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$ un maillage admissible vérifiant (1.18) et (1.19).

Alors il existe une solution $(N_T^n, P_T^n, \Psi_T^n)_{n \ge 0}$ au schéma numérique (1.24)–(1.28), telle que les densités approchées vérifient des bornes inférieures et supérieures uniformes en temps : il existe $D, E > 0$ telles que

$$0 < D \le N_K^n, P_K^n \le E, \quad \forall K \in \mathcal{T}, \quad \forall n \ge 0.$$

Les constantes $D$ et $E$ sont indépendantes des paramètres de discrétisation, et dépendent uniquement des données (conditions initiales et au bord, dopage, longueur de Debye et paramètres $\zeta$ et $\xi$ apparaissant dans (1.18) et (1.19)).

La preuve de ce résultat se fait en deux temps. La première étape consiste à démontrer l'existence d'une solution au schéma telle que les densités sont strictement positives et vérifient une borne supérieure uniforme en temps. La borne inférieure uniforme est ensuite prouvée dans un deuxième temps. Je donne maintenant les principales étapes de cette démonstration.

La preuve de l'existence d'une solution vérifiant la positivité et une borne supérieure uniforme repose sur un argument de degré topologique. Plus précisément, étant donnés $(N_K^0, P_K^0)_{K \in \mathcal{T}}$ et $(N_\sigma^n, P_\sigma^n)_{\sigma \in \mathcal{E}_{ext}^D, n \ge 0}$, nous définissons l'application $\mathcal{H}$ de la manière suivante

$$\mathcal{H}: \quad \begin{array}{ccc} \mathbb{R}^{\theta \times N_T} \times \mathbb{R}^{\theta \times N_T} \times [0, 1] & \to & \mathbb{R}^{\theta \times N_T} \times \mathbb{R}^{\theta \times N_T} \\ ((N_K^n, P_K^n)_{K \in \mathcal{T}, 1 \le n \le N_T}, \kappa) & \mapsto & (((\mathcal{H}_\kappa N)^n_K, (\mathcal{H}_\kappa P)^n_K)_{K \in \mathcal{T}, 1 \le n \le N_T}), \end{array}$$

où

$$(\mathcal{H}_\kappa N)^n_K = \mathrm{m}(K) \frac{N_K^{n+1} - N_K^n}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \hat{\mathcal{F}}_{K,\sigma}^{n+1} + \kappa \mathrm{m}(K) R(N_K^{n+1,+}, P_K^{n+1,+}),$$

$$(\mathcal{H}_\kappa P)^n_K = \mathrm{m}(K) \frac{P_K^{n+1} - P_K^n}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \hat{\mathcal{G}}_{K,\sigma}^{n+1} + \kappa \mathrm{m}(K) R(N_K^{n+1,+}, P_K^{n+1,+}),$$

la notation $x^+ = \max(0, x)$ désignant la partie positive de $x$. Les flux $\hat{\mathcal{F}}_{K,\sigma}^{n+1}$ et $\hat{\mathcal{G}}_{K,\sigma}^{n+1}$ sont définis par

$$\hat{\mathcal{F}}_{K,\sigma}^{n+1} = \tau_{\sigma} \left[ B \left( -D_{K,\sigma} \Psi^{n+1} \right) N_{K}^{n+1,+} - B \left( D_{K,\sigma} \Psi^{n+1} \right) N_{K,\sigma}^{n+1,+} \right],$$

$$\hat{\mathcal{G}}_{K,\sigma}^{n+1} = \tau_{\sigma} \left[ B \left( D_{K,\sigma} \Psi^{n+1} \right) P_{K}^{n+1,+} - B \left( -D_{K,\sigma} \Psi^{n+1} \right) P_{K,\sigma}^{n+1,+} \right],$$

avec $\Psi_T^{n+1}$ satisfaisant

$$(1.36) \quad \left\{ \begin{array}{ll} -\lambda^2 \sum_{\sigma \in \mathcal{E}_K} \tau_{\sigma} D_{K,\sigma} \Psi^{n+1} = \kappa \mathrm{m}(K) (P_K^{n+1} - N_K^{n+1} + C_K) & \forall K \in \mathcal{T}, \\ \Psi_{\sigma}^n = \Psi_{\sigma}^D & \forall \sigma \in \mathcal{E}_{\mathrm{ext}}^D. \end{array} \right.$$
Notons $(S_κ)$ le schéma associé à $\mathcal{H}(\cdot, \kappa)$, c'est-à-dire que $(N_T^n, P_T^n, \Psi_T^n)_{n=1, \dots, N_T}$ est une solution de $(S_κ)$ si $\mathcal{H}((N_K^n, P_K^n)_{K \in \mathcal{T}}, 1 \le n \le N_T}, \kappa) = 0$ et $\Psi_T^n$ vérifie (1.36).

En supposant l'existence d'une solution au schéma $(S_κ)$, nous établissons que les densités correspondantes $N_T^n, P_T^n$ sont strictement positives, et majorées par une borne $E$ uniforme en $\kappa \in [0, 1]$.

Soit $\mathcal{K} = [0, E+1]^{2\theta \times N_T}$ sous-ensemble compact de $\mathbb{R}^{2\theta \times N_T}$. L'application $\mathcal{H}$ est alors uniformément continue sur $\mathcal{K} \times [0, 1]$, et le système non linéaire

$$ \mathcal{H}((N_K^n, P_K^n)_{K \in \mathcal{T}, 1 \le n \le N_T}, \kappa) = 0 $$

ne peut admettre de solution sur $\partial\mathcal{K}$. On en déduit donc que le degré topologique correspondant $\delta(\mathcal{H}, \mathcal{K}, \kappa)$ est constant par rapport à $\kappa$.

Pour terminer cette première étape, il suffit de voir que quand $\kappa = 0$, le problème se réduit à un système linéaire inversible, admettant donc une unique solution, ce qui donne $\delta(\mathcal{H}, \mathcal{K}, 1) = \delta(\mathcal{H}, \mathcal{K}, 0) = 1$, c'est-à-dire que $(S_1)$ admet une solution. Or, du fait de la positivité des densités, nous constatons qu'une solution de $(S_1)$ est en fait une solution du schéma (1.24)–(1.28).

La preuve de la borne supérieure uniforme, tout comme celle de la borne inférieure uniforme, est une adaptation au niveau discret d'une technique d'itération de Moser. Donnons quelques détails sur cette méthode, par exemple pour prouver la borne supérieure.

Afin de ne pas avoir de problème dans le traitement des termes de bord, nous montrons en fait que les quantités

$$ N_{M,K}^n := (N_K^n - M)^+, \quad P_{M,K}^n := (P_K^n - M)^+ $$

sont uniformément bornées, $M$ étant la borne supérieure des conditions au bord. Pour $k \ge 0$, définissons $W_k^n := \|N_{M,T}^n\|_{2k}^{2k} + \|P_{M,T}^n\|_{2k}^{2k}$. En établissant une relation de récurrence entre $W_k^n$ et $W_{k-1}^n$, il est alors possible de prouver que $W_k^n \le c^{2k}$, $c > 0$ étant une constante indépendante de $n$ et de $k$. En pratique, ce résultat utilise une version de l'inégalité d'entropie-production discrète (1.35) adaptée au schéma $(S_κ)$.

Une fois prouvée l'existence d'une solution strictement positive et admettant une borne supérieure uniforme, il reste à établir une borne inférieure uniforme en temps. Ce résultat se démontre dans le même esprit que la borne supérieure, en appliquant une technique d'itération de Moser pour obtenir une borne $L^\infty$ uniforme des quantités

$$ -(\log N_T^n + \bar{m})^{-}, \quad -(\log P_T^n + \bar{m})^{-}, $$

où $x^- = \min(0, x)$ est la partie négative, et $\bar{m} := \max(-\log m, \log E)$, $E$ étant la borne supérieure uniforme déjà établie. Le résultat attendu est ainsi démontré.

Remarquons qu'en général, le principe du maximum n'est pas vérifié, c'est-à-dire que la borne supérieure des densités peut être plus grande que la borne supérieure $M$ des données,
<figure><img src="image_6.png" /><figcaption>FIGURE 1.1 : Évolution des densités approchées dans un cas où $m = 0.5$ et $M = 1$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 1.2 : Bornes supérieure et inférieure en fonction de $\|C\|_\infty$, avec des dopages constants $C = k, k = 0, \dots, 30$.</figcaption></figure>

et de même pour la borne inférieure qui peut être plus petite que $m$ celle des données. Nous pouvons le constater sur la Figure 1.1, qui correspond au cas 1 présenté dans [4] avec un dopage constant $C = 20$, et pour lequel $m = 0.5$ et $M = 1$. L'influence du dopage sur les bornes supérieure et inférieure des densités est illustrée Figure 1.2. D'autres simulations, illustrant également l'influence de la longueur de Debye $\lambda$, sont présentées dans [4]. Du fait de la technique de preuve mise en œuvre, nous ne disposons pas de résultats donnant une dépendance explicite des bornes $D$ et $E$ en fonction des paramètres.

### 1.4.2 Le cas non linéaire sans terme de recombinaison-génération

**Théorème 1.2.** Supposons les Hypothèses H.1 et H.3 vérifiées. Soit $\mathcal{M} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$ un maillage admissible.

Supposons de plus que le pas de temps vérifie

$$ (1.37) \qquad \Delta t \le \frac{\lambda^2}{\|C\|_{\infty}}. $$

Alors il existe une solution $(N_T^n, P_T^n, \Psi_T^n)_{n \ge 0}$ au schéma numérique (1.24)–(1.28), telle que les densités approchées vérifient des bornes inférieures et supérieures : pour tous $K \in$
$T, 0 \le n \le T/\Delta t,$

$$m \exp\left(-\frac{\|C\|_{\infty}}{\lambda^2}T\right) \le m_n \le N_K^n, P_K^n \le M_n \le M \exp\left(\frac{\|C\|_{\infty}}{\lambda^2}T\right),$$

où

$$m_n := m \left( 1 + \frac{\Delta t}{\lambda^2} \|C\|_{\infty} \right)^{-n}, \quad M_n := M \left( 1 - \frac{\Delta t}{\lambda^2} \|C\|_{\infty} \right)^{-n}.$$

En particulier, si $C = 0$, les densités vérifient le principe du maximum :

$$m \le N_K^n, P_K^n \le M, \quad \forall K \in \mathcal{T}, \forall n \ge 0.$$

Ce résultat se démontre en appliquant le théorème de Brouwer à un problème basé sur une linéarisation du schéma (1.24)–(1.28). Soulignons que dans le cas d’un dopage nul $C = 0$, nous obtenons des bornes supérieure et inférieure uniformes en temps pour les densités approchées (qui sont en fait les bornes sur les données initiales et au bord), sans restriction sur le pas de temps.

## 1.5 Convergence exponentielle vers l'équilibre thermique

Je présente maintenant le résultat de convergence vers l'équilibre au niveau discret, sous l'hypothèse que les densités approchées vérifient des bornes uniformes en temps.

**Théorème 1.3.** Supposons les Hypothèses H.1 vérifiées. Soit $\mathcal{M} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$ un maillage admissible vérifiant (1.18). Supposons que les conditions au bord vérifient la condition de compatibilité (1.12). Enfin, plaçons-nous sous l'un des deux jeux d'hypothèses H.2 ou H.3. Si le cas H.3 est considéré, supposons de plus le dopage C nul.

Alors il existe une constante $\beta > 0$ dépendant uniquement des bornes supérieure et inférieure des densités, de $\lambda$, $r$ et du paramètre de régularité du maillage $\xi$, mais pas de la taille du maillage, telle que

$$\mathbb{E}^n \le e^{-\beta t^n} \mathbb{E}^0 \quad \forall n \ge 0.$$

Il existe de plus une constante $c > 0$ dépendant uniquement de la borne inférieure des densités, de $r$ et $\xi$ telle que

$$\|N_{\mathcal{T}}^n - N_{\mathcal{T}}^{eq}\|_0^2 + \|P_{\mathcal{T}}^n - P_{\mathcal{T}}^{eq}\|_0^2 + \|\Psi_{\mathcal{T}}^n - \Psi_{\mathcal{T}}^{eq}\|_0^2 \le c e^{-\beta t^n} \mathbb{E}^0.$$

La preuve de ce théorème est une adaptation au niveau discret de la démarche présentée dans l'introduction de ce chapitre. La première étape est d'établir l'analogue discret de (1.1), ce qui est l'objet de la Proposition 1.1. La seconde étape consiste à obtenir l'analogue discret de (1.2), ce qui se fait grâce à une quantité intermédiaire $\mathbb{F}^n$ comme énoncé dans la proposition suivante.
<figure><img src="image_8.png" /><figcaption>(a) $C = 0, R = 0$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) $C = 0, R = R_{SRH}$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c) $C = 0, R = R_{AU}$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(d) $C \neq 0, R = 0$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(e) $C \neq 0, R = R_{SRH}$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(f) $C \neq 0, R = R_{AU}$</figcaption></figure>

FIGURE 1.3 : Évolution de $\mathbb{E}^n$ dans le cas linéaire avec $C = 0$ et $C \neq 0$, pour différents termes de recombinaison-génération.

**Proposition 1.2.** Sous les hypothèses du théorème 1.3, il existe des constantes $C_{EF}$, $C_{EI} > 0$ dépendant uniquement de $D$, $E$, $\lambda$, $r$ et $\xi$ telles que

$$ \mathbb{E}^n \leq C_{EF} \mathbb{F}^n \leq C_{EI} \mathbb{I}^n, \quad \forall n \geq 0, $$

où

$$ \mathbb{F}^n = \|N_T^n - N_T^{eq}\|_0^2 + \|P_T^n - P_T^{eq}\|_0^2 + \frac{\lambda^2}{2} |\Psi_M^n - \Psi_M^{eq}|_{1,M}^2. $$

C'est dans la preuve de ce résultat que les bornes uniformes en temps sur les densités, établies dans le théorème 1.1 pour le cas isotherme et dans le théorème 1.2 pour le cas non linéaire avec dopage nul, interviennent de manière cruciale. La combinaison des inégalités obtenues dans les propositions 1.1 et 1.2 permet d'obtenir la convergence exponentielle vers 0 de l'entropie relative discrète $\mathbb{E}^n$ quand $n$ tend vers $+\infty$, et la convergence en norme $L^2$ s'en déduit grâce à une inégalité de Poincaré discrète [5, Theorem 4.3].

La figure 1.3 fournit une illustration de la convergence exponentielle vers 0 de l'énergie relative pour une diode à jonction PN, dans le cas linéaire, pour différents choix de dopage et de terme de recombinaison-génération (Auger et Schockley-Read-Hall). Les données utilisées ainsi que d'autres illustrations numériques sont présentées dans [3].

## 1.6 Perspectives

**Modèles avec diffusion non linéaire.** Le résultat de convergence exponentielle des solutions approchées vers l'équilibre thermique a été rigoureusement démontré sous les
<figure><img src="image_14.png" /><figcaption>FIGURE 1.4 : Évolution de $\mathbb{E}^n$ dans le cas non linéaire avec $r(s) = s^{5/3}$.</figcaption></figure>

hypothèses H.2, c'est-à-dire dans le cas d'une diffusion linéaire, avec dopage et terme de recombinaison-génération. Dans le cas non linéaire correspondant aux hypothèses H.3, la preuve complète est établie uniquement dans le cas sans dopage. Il manque l'obtention de bornes $L^\infty$ uniformes en temps sur les densités approchées pour avoir le résultat dans le cadre général, qui semble par ailleurs vrai au vu des résultats numériques obtenus sur la figure 1.4. Le résultat partiel obtenu dans le théorème 1.2 n'est pas optimal puisqu'en pratique la condition (1.37) sur le pas de temps ne semble pas nécessaire dans la mise en œuvre numérique, et des bornes uniformes en temps semblent exister pour les densités approchées, même dans le cas général. Un travail futur pourra donc consister à étendre notre résultat à ce cadre non linéaire.

**Autres discrétisations volumes finis.** Au regard des résultats présentés dans [CHH20], une extension des résultats de ce chapitre à une classe plus large de flux monotones à deux points semble possible. En effet, dans l'article [CHH20], un résultat de convergence exponentielle fondé sur des estimations d'entropie-production discrètes est établi pour des équations de type Fokker-Planck et milieux poreux, pour une famille de B-schémas incluant les flux upwind, centrés et Scharfetter-Gummel. Le cas du système de dérive-diffusion n'est pour l'instant pas traité mais pourrait être une extension naturelle des résultats que je présente dans ce chapitre.

**Modèles de semi-conducteurs organiques.** Depuis quelques années, l'utilisation croissante de semi-conducteurs organiques a conduit à de nombreux travaux sur les modèles associés, qui peuvent être vus comme des généralisations des modèles de semi-conducteurs classiques. Il pourrait être intéressant d'étendre le type d'étude menée dans ce chapitre à ces nouveaux modèles. Citons par exemple le système de dérive-diffusion non isotherme, correspondant au système (1.3)-(1.5) couplé avec l'équation de la chaleur, pour lequel des extensions du schéma de Scharfetter-Gummel généralisé (1.24)-(1.26) sont étudiées [Kan20]. Nous pouvons également mentionner les travaux récents sur un modèle de cellules photovoltaïques à pérovskite, consistant en trois équations de dérive-diffusion (une
espèce supplémentaire de porteurs de charge est considérée pour prendre en compte la migration des ions) couplées avec l'équation de Poisson [ACHFH]. L'analyse du schéma volumes finis proposé pour ce modèle repose sur une inégalité d'entropie-dissipation, et les résultats numériques font apparaître une convergence exponentielle vers l'équilibre thermique associé.
<figure><img src="image_15.png" /></figure>
Schémas numériques pour le modèle de transport d'énergie

Ce chapitre est consacré à la construction et l'analyse de schémas numériques pour des modèles de transport d'énergie unipolaires. Contrairement au modèle de dérive-diffusion considéré dans le chapitre 1, les modèles de transport d'énergie pour les semi-conducteurs prennent en compte la température, et permettent par exemple d'observer l'effet d'électrons chauds.

Le système considéré est constitué de deux équations de continuité sur la densité d'électrons, notée dans ce chapitre $\rho_1$, et la densité d'énergie interne $\rho_2$, couplées avec l'équation de Poisson sur le potentiel électrique $\Psi$. En se plaçant dans le cadre adopté dans [DGJ97], ces densités sont considérées comme des fonctions des variables entropiques $\mathbf{u} = (u_1, u_2)$, où

$$u_1 = \frac{\mu}{T}, \quad u_2 = -\frac{1}{T},$$

$\mu$ étant le potentiel chimique et $T$ la température. Étant donné un ouvert borné $\Omega \subset \mathbb{R}^d$, $d \ge 1$, et $T_{max} > 0$ un temps final, le système de transport d'énergie s'écrit alors : pour tous $(t, x) \in (0, T_{max}) \times \Omega$,

$$(2.1) \qquad \partial_t \rho_1(\mathbf{u}) + \text{div} J_1(\mathbf{u}) = 0,$$

$$(2.2) \qquad \partial_t \rho_2(\mathbf{u}) + \text{div} J_2(\mathbf{u}) = \nabla \Psi \cdot J_1(\mathbf{u}) + W(\mathbf{u}),$$

$$(2.3) \qquad -\lambda^2 \Delta \Psi = C(x) - \rho_1(\mathbf{u}),$$

où les densités de courant d'électrons $J_1$ et d'énergie $J_2$ sont données par :

$$J_1(\mathbf{u}) = -L_{11}(\mathbf{u})(\nabla u_1 + u_2 \nabla \Psi) - L_{12}(\mathbf{u}) \nabla u_2,$$

$$J_2(\mathbf{u}) = -L_{21}(\mathbf{u})(\nabla u_1 + u_2 \nabla \Psi) - L_{22}(\mathbf{u}) \nabla u_2.$$

Les coefficients $(L_{ij}(\mathbf{u}))_{1 \le i,j \le 2}$ forment une matrice $\mathbb{L}(\mathbf{u})$ supposée symétrique et unifor-mément définie positive. Comme dans le chapitre précédent, $C(x)$ est le profil de dopage,
et $\lambda$ la longueur de Debye adimensionnée. Ce système (2.1)–(2.3) est complété par des conditions initiales $\mathbf{u}_0 = (u_{1,0}, u_{2,0})$ et des conditions au bord mixtes de type Dirichlet Neumann, supposées indépendantes du temps :

$$ 
\begin{gathered}
u_1 = u_1^D, \quad u_2 = u_2^D, \quad \Psi = \Psi^D \quad \text{sur } \Gamma^D, \\
J_1 \cdot \nu = J_2 \cdot \nu = \nabla \Psi \cdot \nu = 0 \quad \text{sur } \Gamma^N,
\end{gathered}
$$

où $\partial\Omega = \Gamma^D \cup \Gamma^N$ et $\nu$ est le vecteur unitaire sortant normal à $\partial\Omega$.

Notons que selon les définitions choisies pour $L(\mathbf{u})$ et $W(\mathbf{u})$, ce cadre général inclut les modèles de Chen et de Lyumkis [JØ1].

Dans [DGJ97], l'étude de ce système d'EDP est fondée sur un changement de variables permettant de symétriser le problème :

$$ (2.4) \qquad w_1 = u_1 + u_2 \Psi, \quad w_2 = u_2. $$

En appliquant ce changement de variables, le système (2.1)–(2.3) est équivalent à

$$ (2.5) \qquad \partial_t b_1(\mathbf{w}, \Psi) + \operatorname{div} I_1(\mathbf{w}, \Psi) = 0, $$

$$ (2.6) \qquad \partial_t b_2(\mathbf{w}, \Psi) + \operatorname{div} I_2(\mathbf{w}, \Psi) = \tilde{W}(\mathbf{w}) - \partial_t \Psi b_1(\mathbf{w}, V), $$

$$ (2.7) \qquad -\lambda^2 \Delta \Psi = C - b_1(\mathbf{w}, \Psi), $$

où $b_1(\mathbf{w}, \Psi) = \rho_1(\mathbf{u})$, $b_2(\mathbf{w}, \Psi) = \rho_2(\mathbf{u}) - \Psi \rho_1(\mathbf{u})$ et $\tilde{W}(\mathbf{w}) = W(\mathbf{u})$. Les courants symétrisés s'écrivent alors :

$$ (2.8) \qquad \begin{aligned} I_1(\mathbf{w}, \Psi) &= -D_{11}(\mathbf{w}, \Psi)\nabla w_1 - D_{12}(\mathbf{w}, \Psi)\nabla w_2, \\ I_2(\mathbf{w}, \Psi) &= -D_{21}(\mathbf{w}, \Psi)\nabla w_1 - D_{22}(\mathbf{w}, \Psi)\nabla w_2, \end{aligned} $$

où la nouvelle matrice de diffusion $D(\mathbf{w}, \Psi) = (D_{ij}(\mathbf{w}, \Psi))_{1\le i,j\le 2}$ est définie par

$$ D(\mathbf{w}, \Psi) = \mathbb{P}(\Psi)^T L(\mathbf{u}) \mathbb{P}(\Psi), \quad \text{avec } \mathbb{P}(\Psi) = \begin{pmatrix} 1 & -\Psi \\ 0 & 1 \end{pmatrix}. $$

À l'aide de ce modèle dual (2.5)–(2.7), une propriété d'entropie-production est établie, sur laquelle repose l'analyse du modèle (2.1)–(2.3) proposée dans [DGJ97].

Plus précisément, faisons les hypothèses suivantes, qui seront valables pour le reste du chapitre.

**Hypothèses H.4.** • La condition au bord $u_2^D$ est constante sur $\Gamma^D$ et $u_2^D < 0$.

• Les conditions au bord sont à l'équilibre thermique : $\nabla w_1^D = \nabla w_2^D = 0$.

• Le terme source $W(\mathbf{u})$ vérifie : pour tous $\mathbf{u} \in \mathbb{R}^2$ et $u_2^D < 0$,

$$ W(\mathbf{u})(u_2 - u_2^D) \le 0. $$
* La fonction $\rho = (\rho_1, \rho_2)$ est telle que

$$ - \rho \in W^{1,\infty}(\mathbb{R}^2, \mathbb{R}^2), $$

$$ - \text{il existe } C_0 > 0 \text{ tel que pour tous } \mathbf{u}, \mathbf{v} \in \mathbb{R}^2, $$

$$ (\rho(\mathbf{u}) - \rho(\mathbf{v})) \cdot (\mathbf{u} - \mathbf{v}) \geq C_0 |\mathbf{u} - \mathbf{v}|^2, $$

$$ - \text{il existe } \chi \in C^1(\mathbb{R}^2, \mathbb{R}) \text{ convexe telle que } \rho(\mathbf{u}) = \nabla_{\mathbf{u}}\chi(\mathbf{u}). $$

Définissions l'entropie par

$$ (2.9) \quad S(t) = \int_{\Omega} \left[ \rho(\mathbf{u}) \cdot (\mathbf{u} - \mathbf{u}^D) - (\chi(\mathbf{u}) - \chi(\mathbf{u}^D)) \right] dx - \frac{\lambda^2}{2} u_2^D \int_{\Omega} \left| \nabla(\Psi - \Psi^D) \right|^2 dx. $$

Il est prouvé dans [DGJ97] que sous les hypothèses H.4, la propriété d'entropie-production suivante est vérifiée :

$$ (2.10) \quad \frac{d}{dt} S(t) = - \int_{\Omega} \nabla \mathbf{w}^T \mathbb{D} \nabla \mathbf{w} + \int_{\Omega} W(\mathbf{u})(u_2 - u_2^D) \leq 0. $$

À partir de cette inégalité, des estimations *a priori* sur $\mathbf{u}$, $\mathbf{w}$ et $\Psi$ peuvent être démontrées, permettant notamment d'établir l'existence d'une solution au système (2.1)–(2.3).

Un certain nombre de schémas numériques pour les modèles de transport d'énergie ont été proposés dans la littérature [DJP00, Fou02, HJP03, Rom07, CH09], principalement dans le cas stationnaire, mais sans analyse de convergence. L'objectif des travaux présentés dans ce chapitre est de construire des schémas numériques pour le modèle de transport d'énergie de telle sorte que la méthode d'entropie décrite brièvement ci-dessus soit adaptable au cadre discret.

Dans la section 2.1, je présente le schéma volumes finis de type TPFA (Two Points Flux Approximation) que nous avons proposé avec Claire Chainais-Hillairet et Hélène Mathis [7, 8]. Plus précisément, nous construisons des schémas pour le système en variables primales (2.1)–(2.3) et pour le système en variables duales (2.5)–(2.7) de telle sorte qu'ils soient équivalents, ce qui permet d'établir l'analogue discret de l'inégalité d'entropie (2.10). À partir de cette propriété, l'analyse numérique du schéma peut être conduite dans la section 2.2 : estimations *a priori*, existence d'une solution au schéma et comportement en temps long. Dans la section 2.3, je présente finalement l'extension de ces résultats au cadre DDFV (Discrete Duality Finite Volume), obtenue en collaboration avec Giulia Lissoni et Hélène Mathis [12].

## 2.1 Structure entropique discrète

Je considère une discrétisation de type volume fini classique, en utilisant les notations introduites dans la section 1.2 du chapitre 1. Comme dans le chapitre 1, nous supposons la condition de régularité suivante vérifiée :

$$ (2.11) \quad \exists \xi > 0 \text{ tel que } d(x_K, \sigma) \geq \xi d_\sigma, \quad \forall K \in \mathcal{T}, \quad \forall \sigma \in \mathcal{E}_K. $$
Les conditions initiales et au bord sont discrétisées classiquement :

$$ 
\begin{aligned}
(u_{1,K}^0, u_{2,K}^0) &= \frac{1}{m(K)} \int_K (u_{1,0}(x), u_{2,0}(x)) \, dx, \quad \forall K \in \mathcal{T}, \\
(u_{1,\sigma}^D, u_{2,\sigma}^D, \Psi_{\sigma}^D) &= \frac{1}{m(\sigma)} \int_{\sigma} (u_1^D(\gamma), u_2^D(\gamma), \Psi^D(\gamma)) \, d\gamma, \quad \forall \sigma \in \mathcal{E}_{\text{ext}}^D,
\end{aligned}
$$

et les valeurs au bord sont définies par

$$ u_{1,\sigma}^n = u_{1,\sigma}^D, \quad u_{2,\sigma}^n = u_{2,\sigma}^D, \quad \Psi_{\sigma}^n = \Psi_{\sigma}^D, \quad \forall \sigma \in \mathcal{E}_{\text{ext}}^D, \quad \forall n \ge 0. $$

La difficulté principale pour définir un schéma TPFA pour le système (2.1)–(2.3) est la discrétisation du terme d'effet Joule $\nabla\Psi \cdot J_1$. En effet, ce terme contient le produit de deux gradients, et dans le cadre TPFA seule une approximation unidimensionnelle du gradient selon la normale à l'interface est accessible, ne permettant pas d'obtenir directement une discrétisation consistante pour ce produit. Pour contourner cet obstacle et faciliter l'obtention du schéma équivalent pour le système en variables duales (2.5)–(2.7), le terme d'effet Joule est reformulé comme suit [CCC19] :

$$ (2.12) \qquad \nabla\Psi \cdot J_1 = \operatorname{div}(\Psi J_1) - \Psi \operatorname{div}(J_1). $$

Nous considérons le schéma volumes finis implicite en temps suivant pour le système en variables primales (2.1)–(2.3) : pour tous $n \ge 0$ et $K \in \mathcal{T}$,

$$ (2.13) \qquad m(K) \frac{\rho_{1,K}^{n+1} - \rho_{1,K}^n}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \mathcal{F}_{1,K,\sigma}^{n+1} = 0, $$

$$ (2.14) \qquad m(K) \frac{\rho_{2,K}^{n+1} - \rho_{2,K}^n}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \mathcal{F}_{2,K,\sigma}^{n+1} = \sum_{\sigma \in \mathcal{E}_K} \Psi_{\sigma}^{n+1} \mathcal{F}_{1,K,\sigma}^{n+1} - \Psi_{K}^{n+1} \sum_{\sigma \in \mathcal{E}_K} \mathcal{F}_{1,K,\sigma}^{n+1} + m(K) W_K^{n+1}, $$

$$ (2.15) \qquad -\lambda^2 \sum_{\sigma \in \mathcal{E}_K} \tau_{\sigma} D_{K,\sigma} \Psi^{n+1} = m(K)(C_K - \rho_{1,K}^{n+1}), $$

où

$$ \rho_{i,K}^{n+1} = \rho_i(\mathbf{u}_K^{n+1}), \quad i = 1, 2, \quad \text{et} \quad W_K^{n+1} = W(\mathbf{u}_K^{n+1}), \quad \text{pour tous } K \in \mathcal{T}. $$

Les flux numériques sont définis par

$$ 
\begin{aligned}
\mathcal{F}_{1,K,\sigma}^{n+1} &= -\tau_{\sigma}(L_{11,\sigma}^n(D_{K,\sigma}u_1^{n+1} + u_{2,\sigma}^{n+1}D_{K,\sigma}\Psi^{n+1}) + L_{12,\sigma}^n D_{K,\sigma}u_2^{n+1}), \\
\mathcal{F}_{2,K,\sigma}^{n+1} &= -\tau_{\sigma}(L_{12,\sigma}^n(D_{K,\sigma}u_1^{n+1} + u_{2,\sigma}^{n+1}D_{K,\sigma}\Psi^{n+1}) + L_{22,\sigma}^n D_{K,\sigma}u_2^{n+1}).
\end{aligned}
$$

Les coefficients $L_{ij,\sigma}^n$ sont des approximations des coefficients de la matrice $\mathbb{L}$ à l'interface $\sigma$. La matrice discrète $\mathbb{L}_{\sigma}^n = (L_{ij,\sigma}^n)_{1 \le i,j \le 2}$ est donnée par

$$ \mathbb{L}_{\sigma}^n = \mathbb{L} \left( \frac{\mathbf{u}_K^n + \mathbf{u}_{K,\sigma}^n}{2} \right), \quad \text{pour tous } K \in \mathcal{T}, \sigma \in \mathcal{E}_K. $$
Pour achever la construction de ce schéma, il reste à définir les approximations aux interfaces $\Psi_{\sigma}^{n+1}$ et $u_{2,\sigma}^{n+1}$, ce qui est fait dans la proposition 2.1 de sorte à obtenir l'équivalence avec un schéma pour le système en variables duales.

En appliquant le changement de variables (2.4), le schéma numérique suivant est établi :

$$ (2.17) \qquad m(K) \frac{b_{1,K}^{n+1} - b_{1,K}^{n}}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \mathcal{G}_{1,K,\sigma}^{n+1} = 0, $$

$$ (2.18) \qquad m(K) \frac{b_{2,K}^{n+1} - b_{2,K}^{n}}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \mathcal{G}_{2,K,\sigma}^{n+1} = m(K) \tilde{W}_K^{n+1} - m(K) \frac{\Psi_K^{n+1} - \Psi_K^n}{\Delta t} b_{1,K}^n, $$

$$ (2.19) \qquad -\lambda^2 \sum_{\sigma \in \mathcal{E}_K} \tau_{\sigma} D_{K,\sigma} \Psi^{n+1} = m(K)(C_K - b_{1,K}^{n+1}), $$

où $\tilde{W}_K^{n+1} = W_K^{n+1} = \tilde{W}(w_K^{n+1}, \Psi_K^{n+1})$ et les flux numériques sont donnés par :

$$ (2.20) \qquad \begin{aligned} \mathcal{G}_{1,K,\sigma}^{n+1} &= \mathcal{F}_{1,K,\sigma}^{n+1}, \\ \mathcal{G}_{2,K,\sigma}^{n+1} &= \mathcal{F}_{2,K,\sigma}^{n+1} - \Psi_{\sigma}^{n+1} \mathcal{F}_{1,K,\sigma}^{n+1}. \end{aligned} $$

Le point crucial est de s'assurer que les flux numériques obtenus $\mathcal{G}_{i,K,\sigma}$, $i = 1, 2$ sont effectivement des approximations des courants symétrisés $I_i$ définis par (2.8) pour avoir des schémas pour les systèmes en variables primales et duales, qui sont équivalents entre eux. Cette propriété est obtenue grâce à un choix judicieux des valeurs d'interface $\Psi_{\sigma}^{n+1}$ et $u_{2,\sigma}^{n+1}$, comme le détaille la proposition suivante.

**Proposition 2.1.** Le schéma (2.13)-(2.15) est complété avec les définitions des quantités aux interfaces $(\Psi_{\sigma}^{n+1})_{\sigma \in \mathcal{E}, n \ge 0}$ et $(u_{2,\sigma}^{n+1})_{\sigma \in \mathcal{E}, n \ge 0}$. Nous distinguons deux cas :

* Cas 1 : schéma centré. Pour tous $\sigma \in \mathcal{E}$ et $n \ge 0$,

$$ (2.21) \qquad u_{2,\sigma}^{n+1} = \frac{u_{2,K}^{n+1} + u_{2,K,\sigma}^{n+1}}{2} \quad \text{et} \quad \Psi_{\sigma}^{n+1} = \frac{\Psi_K^{n+1} + \Psi_{K,\sigma}^{n+1}}{2}. $$

* Cas 2 : schéma upwind. Pour tous $\sigma \in \mathcal{E}$ et $n \ge 0$,

$$ (2.22) \qquad u_{2,\sigma}^{n+1} = \begin{cases} u_{2,K,\sigma}^{n+1}, & \text{si } D_{K,\sigma} \Psi^{n+1} > 0, \\ u_{2,K}^{n+1}, & \text{si } D_{K,\sigma} \Psi^{n+1} \le 0, \end{cases} \quad \text{et} \quad \Psi_{\sigma}^{n+1} = \min(\Psi_K^{n+1}, \Psi_{K,\sigma}^{n+1}). $$

Dans ces deux cas, le schéma (2.13)-(2.15) en variables primales est équivalent au schéma (2.17)-(2.19) où les flux numériques sont définis par

$$ (2.23) \qquad \mathcal{G}_{i,K,\sigma}^{n+1} = -\tau_{\sigma} (\mathbb{D}_{i1,\sigma}^{*} D_{K,\sigma} w_{1}^{n+1} + \mathbb{D}_{i2,\sigma}^{*} D_{K,\sigma} w_{2}^{n+1}), \quad i = 1, 2, $$

les coefficients $(\mathbb{D}_{ij,\sigma}^{*})_{1 \le i,j \le 2}$ étant ceux de la matrice suivante :

$$ (2.24) \qquad \mathbb{D}_{\sigma}^{*} = (\mathbb{P}_{\sigma}^{n+1})^{T} \mathbb{L}_{\sigma}^{n} \mathbb{P}_{\sigma}^{n+1} \quad \text{avec} \quad \mathbb{P}_{\sigma}^{n+1} = \begin{pmatrix} 1 & -\Psi_{\sigma}^{n+1} \\ 0 & 1 \end{pmatrix}. $$
<figure><img src="image_16.png" /><figcaption>FIGURE 2.1 : Nombre d'itérations de la méthode de Newton au cours du temps pour différents choix de discrétisation.</figcaption></figure>

Nous avons donc à notre disposition plusieurs discrétisations possibles du système de transport d'énergie, et donc des choix à opérer entre les options suivantes :

* variables primales ou duales,

* discrétisation centrée ou upwind.

La discrétisation temporelle étant implicite, tous ces schémas conduisent en pratique à résoudre un système non linéaire à chaque pas de temps, ce qui est réalisé par une méthode de Newton. Afin de choisir la discrétisation la plus performante, nous avons donc réalisé des comparaisons entre les différents choix, pour différents cas tests unidimensionnels. Une illustration en est donnée sur la figure 2.1, et une étude numérique plus détaillée est proposée dans [8].

Nous avons constaté que les variables primales semblent plus adaptées numériquement, ce qui peut peut-être s'expliquer par le manque de régularité des fonctions impliquées dans la méthode de Newton pour le schéma en variables duales. Concernant le choix entre discrétisation centrée ou upwind, les performances sont comparables dans le cas des variables primales.

D'un point de vue théorique, cette construction de deux schémas équivalents pour les systèmes en variables primales (2.1)–(2.3) et duales (2.5)–(2.7) est essentielle. Définissons l'entropie discrète de la manière suivante :

$$ (2.25) \qquad S^n = \sum_{K \in \mathcal{T}} \mathrm{m}(K) \left[ \rho_K^n \cdot (\mathbf{u}_K^n - \mathbf{u}_K^D) - (\chi(\mathbf{u}_K^n) - \chi(\mathbf{u}_K^D)) \right] \\ - \frac{\lambda^2}{2} u_2^D \sum_{\sigma \in \mathcal{E}} \tau_{\sigma} \left| D_{\sigma}(\Psi^n - \Psi^D) \right|^2. $$
L'équivalence des schémas en variables primales et duales permet d'obtenir un analogue discret de l'inégalité d'entropie-production (2.10).

**Proposition 2.2.** Supposons les hypothèses H.4 vérifiées. Soit $(\mathbf{u}_K^n, \Psi_K^n)_{K \in \mathcal{T}, n \ge 0}$ une solution du schéma en variables primales (2.13)-(2.15) complété par la définition des discrétisations aux interfaces (2.21) ou (2.22).

Alors l'entropie discrète vérifie l'inégalité suivante : pour tous $n \ge 0$,

$$ (2.26) \quad \frac{S^{n+1} - S^n}{\Delta t} \le - \sum_{\sigma \in \mathcal{E}} \tau_{\sigma} (D_{K,\sigma} \mathbf{w}^{n+1})^T \mathbb{D}_{\sigma}^* D_{K,\sigma} \mathbf{w}^{n+1} + \sum_{K \in \mathcal{T}} \mathrm{m}(K) W_K^{n+1} (w_{2,K}^{n+1} - w_{2,K}^D) \le 0. $$

## 2.2 Analyse numérique du schéma

### 2.2.1 Estimations a priori

Les estimations a priori suivantes constituent la première étape dans l'étude des schémas.

**Proposition 2.3.** Supposons les hypothèses H.4 vérifiées. Soit $(\mathbf{u}_K^n, \Psi_K^n)_{K \in \mathcal{T}, n \ge 0}$ une solution du schéma en variables primales (2.13)-(2.15) complété par la définition des discrétisations aux interfaces (2.21) ou (2.22).

Alors il existe des constantes $C_1, C_2 > 0$ dépendant uniquement des données, mais pas des paramètres de discrétisation, telles que

$$ (2.27) \quad \sup_{n=0,...,N_T} \left( \|u_{1,T}^n - u_{1,T}^D\|_2 + \|u_{2,T}^n - u_{2,T}^D\|_2 + |\Psi_M^{n+1} - \Psi_M^D|^2 \right) \le C_1, $$

$$ (2.28) \quad \sum_{n=0}^{N_T-1} \Delta t \left( |w_{1,M}^{n+1}|_{1,M}^2 + |w_{2,M}^{n+1}|_{1,M}^2 \right) \le C_2. $$

L'estimation (2.27) s'obtient en utilisant l'inégalité d'entropie-production discrète (2.26) et les hypothèses sur $\rho$. Pour établir (2.28), le point crucial est d'avoir une borne $L^\infty$ uniforme sur le potentiel approché $\Psi_T$. Ce résultat est établi dans un cadre général pour les schémas TPFA dans [CCHFG21, Proposition A1]. Il permet de montrer que la matrice $\mathbb{D}_\sigma^*$ est uniformément définie positive, et donc de déduire (2.28) en sommant sur $n$ l'inégalité d'entropie-production (2.26).

### 2.2.2 Existence d'une solution au schéma

En adaptant au cadre discret la preuve proposée dans [DGJ97, Lemma 3.2], fondée sur une application du théorème du point fixe de Leray-Schauder, nous obtenons l'existence de solutions aux schémas en variables primales et duales. La démonstration utilise des arguments similaires à ceux développés dans les preuves des propositions 2.2 et 2.3.
**Théorème 2.1.** Soit $\mathcal{M} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$ un maillage admissible vérifiant la contrainte de régularité (2.11). Sous les hypothèses H.4, il existe une solution $(u_{1,\mathcal{T}}^n, u_{2,\mathcal{T}}^n, \Psi_{\mathcal{T}}^n)_{n\ge 0}$ au schéma en variables primales (2.13)-(2.3). Par équivalence, le schéma en variables duales (2.5)-(2.7) admet donc également une solution $(w_{1,\mathcal{T}}^n, w_{2,\mathcal{T}}^n, \Psi_{\mathcal{T}}^n)_{n\ge 0}$.

### 2.2.3 Comportement en temps long

J'établis maintenant la convergence exponentielle en temps long de la solution approchée du système de transport d'énergie vers une approximation de l'équilibre thermique. Nous nous plaçons dans le cadre considéré dans [DGJ97] pour le modèle continu. Sous les hypothèses H.4, supposons de plus les conditions de compatibilité suivantes sur les données au bord :

$$ (2.29) \qquad W(\mathbf{u}^D) = 0, $$

$$ (2.30) \qquad -\lambda^2 \Delta \Psi^D = C - \rho_1(\mathbf{u}^D). $$

L'équilibre thermique $(\mathbf{w}^{eq}, \Psi^{eq})$ est alors défini comme étant l'unique solution du problème elliptique non linéaire suivant :

$$ (2.31) \qquad -\lambda^2 \Delta \Psi^{eq} = C - b_1(\mathbf{w}^{eq}, \Psi^{eq}), \quad \mathbf{w}^{eq} = \mathbf{w}^D \quad \text{sur } \Omega, $$

avec $\Psi^{eq} = \Psi^D$ sur $\Gamma^D$ et $\nabla\Psi^{eq} \cdot \nu = 0$ sur $\Gamma^N$.

Nous considérons une discrétisation TPFA classique de cet équilibre thermique : pour tout $K \in \mathcal{T}$,

$$ (2.32) \qquad -\lambda^2 \sum_{\sigma \in \mathcal{E}_K} \tau_{\sigma} D_{K,\sigma} \Psi^{eq} = m(K) (C_K - b_1(\mathbf{w}_K^{eq}, \Psi_K^{eq})), \quad \mathbf{w}_K^{eq} = \mathbf{w}_K^D, $$

avec $\Psi_{\sigma}^{eq} = \Psi_{\sigma}^{D}$ pour tous $\sigma \in \mathcal{E}_{ext}^{D}$.

L'étude du comportement en temps long des solutions approchées est fondée, comme dans le chapitre 1, sur une estimation d'entropie relative par rapport à l'équilibre thermique définie par

$$ (2.33) \qquad \begin{aligned} \tilde{S}^n = \sum_{K \in \mathcal{T}} & m(K) [\rho_K^n \cdot (\mathbf{u}_K^n - \mathbf{u}_K^{eq}) - (\chi(\mathbf{u}_K^n) - \chi(\mathbf{u}_K^{eq}))] \\ & - \frac{\lambda^2}{2} u_2^D \sum_{\sigma \in \mathcal{E}} \tau_{\sigma} |D_{\sigma}(\Psi^n - \Psi^{eq})|^2. \end{aligned} $$

En adaptant les preuves des propositions 2.2 et 2.3, nous établissons la double inégalité d'entropie-production discrète suivante :

$$ \tilde{S}^{n+1} - \tilde{S}^n \leq -\Delta t \left( |w_{1,M}^{n+1}|_{1,M}^2 + |w_{2,M}^{n+1}|_{1,M}^2 \right) \leq -\frac{\Delta t}{C_{SI}} \tilde{S}^{n+1}, $$

de laquelle nous déduisons la convergence exponentielle vers l'équilibre thermique quand $t \to +\infty$.
<figure><img src="image_17.png" /><figcaption>FIGURE 2.2 : État stationnaire correspondant à un transistor de type MESFET.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 2.3 : Évolution au cours du temps de l'entropie relative $\tilde{S}^n$ et des normes $\|\mathbf{u}_T^n - \mathbf{u}_T^{eq}\|_2 + \|\Psi_M^n - \Psi_M^{eq}\|_{1,M}$ dans le cas d'un transistor MESFET.</figcaption></figure>

**Théorème 2.2.** Soit $\mathcal{M} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$ un maillage admissible vérifiant la contrainte de régularité (2.11). Sous les hypothèses H.4, (2.29) et (2.30), il existe $\gamma > 0$ dépendant uniquement des données et du paramètre de régularité du maillage $\xi > 0$ tel que

$$ (2.34) \qquad \tilde{S}^n \le \tilde{S}^0 e^{-\gamma t^n}, \quad \forall n \ge 0. $$

De plus, il existe $c > 0$ dépendant uniquement des données et de $\xi$ tel que

$$ (2.35) \qquad \|\mathbf{u}_T^n - \mathbf{u}_T^{eq}\|_2^2 + \|\Psi_M^n - \Psi_M^{eq}\|_{1,M}^2 \le c \tilde{S}^0 e^{-\gamma t^n} \quad \forall n \ge 0. $$

Si d'un point de vue théorique la convergence exponentielle en temps long vers un état stationnaire est démontrée uniquement dans le cas très particulier de l'équilibre thermique, elle semble cependant bien vérifiée dans des cas plus généraux. Par exemple, dans le cas d'un transistor de type MESFET [HJP03], nous observons sur la figure 2.3 une convergence exponentielle des solutions approchées vers l'état stationnaire illustré sur la figure 2.2.
## 2.3 Extension au cadre DDFV

Pour des cas tests physiquement réalistes, la différence d'épaisseur entre différentes parties du domaine de simulation peut nécessiter l'usage de maillages localement raffinés (voir par exemple [HJP04]). Dans ce contexte, il peut s'avérer compliqué de construire un maillage admissible, requis pour appliquer les schémas de type TPFA présentés depuis le début de ce chapitre. Je m'intéresse donc dans cette section à une généralisation au cadre DDFV des idées présentées jusqu'ici dans le cadre TPFA. En effet, les méthodes de type DDFV permettent de considérer des maillages très généraux. Contrairement aux schémas TPFA, elles nécessitent des inconnues à la fois aux centres et aux sommets des volumes de contrôle, ce qui permet de construire une approximation bidimensionnelle des gradients, en dualité discrète avec les opérateurs de divergence correspondants. Ce cadre est donc bien adapté pour transposer au niveau discret les propriétés des opérateurs différentiels continus [ABH07].

Un schéma DDFV pour le modèle de transport d'énergie avait déjà été proposé dans [CH09]. Celui que nous considérons dans [12] est différent, construit de telle sorte que les idées développées dans le cadre TPFA (équivalence entre les discrétisations en variables entropiques primales et duales, inégalité d'entropie-production) puissent être retrouvées dans le cadre DDFV général.

Comme dans la section 2.1, le terme d'effet Joule est réécrit sous la forme (2.12). Ici, nous reformulons également le terme $u_2 \nabla \Psi$ comme

$$ (2.36) \qquad u_2 \nabla \Psi = \operatorname{div}(u_2 \Psi) - \Psi \nabla u_2. $$

Cette reformulation permet de prouver directement l'équivalence entre les schémas en variables primales et duales. En effet, contrairement au cadre continu, la version discrète de (2.36) n'est en générale pas vérifiée. Dans le cadre TPFA développé dans la section 2.1, les approximations de $u_2$ et $\nabla \Psi$ aux interfaces sont justement construites de sorte à obtenir un analogue discret de (2.36). Ici, nous choisissons plutôt d'avoir une unique définition de l'opérateur de reconstruction par diamant (correspondant à une approximation aux interfaces), et de discrétiser directement le terme $u_2 \nabla \Psi$ sous sa forme reformulée (2.36).

En construisant ainsi le schéma DDFV pour le système en variables primales (2.1)–(2.3), nous obtenons dans ce cadre généralisé :

* un schéma DDFV en variables duales équivalent à celui en variables primales (analogue dans le cadre DDFV de la proposition 2.1),

* une inégalité d'entropie-production discrète (analogue de la proposition 2.2)

* une estimation a priori $L^\infty$ en temps, $L^2$ en espace pour les approximations de $u - u^D$ et $\nabla(\Psi - \Psi^D)$ (analogue de (2.27)),

* l'existence de solutions aux schémas en variables primales et duales (analogue du théorème 2.1).
## 2.4 Perspectives

**Convergence des schémas TPFA.** À partir des estimations *a priori* énoncées dans la proposition 2.3, on déduit par inégalité de Poincaré discrète les bornes suivantes :

$$ \sum_{n=0}^{N_T-1} \Delta t (\|w_{1,M}^{n+1}\|_{1,M}^2 + \|w_{2,M}^{n+1}\|_{1,M}^2) \le c, \quad \sup_{n=0,...,N_T} (\|\Psi_M^{n+1} - \Psi_M^D\|_{1,M}) \le c. $$

Ces estimations sont la première étape pour effectuer une preuve complète de convergence des schémas TPFA pour le modèle de transport d'énergie. La principale difficulté pour montrer rigoureusement ce résultat provient de la non linéarité des densités $\rho_1(\mathbf{u})$, $\rho_2(\mathbf{u})$. Dans ce contexte, il semble nécessaire d'utiliser le résultat général donné dans [ACM17] pour déduire de la compacité en espace et en temps à partir des estimations *a priori* en espace. Ensuite, l'identification de la limite obtenue comme une solution faible du problème continu devrait s'obtenir de manière assez standard.

**Analyse des schémas DDFV.** Dans le cadre DDFV, nous obtenons uniquement un analogue de l'estimation (2.27) pour les approximations de $\mathbf{u}$ et $\Psi$. À ce stade, nous ne pouvons pas obtenir d'estimation discrète $L^2$ en temps, $H^1$ en espace de $w_M$. En effet, ce résultat nécessite une borne $L^\infty$ sur le potentiel $\Psi_T$, qui est obtenu dans le cadre TPFA en tirant parti de la structure particulière de la matrice associée au schéma. Il n'existe pas de résultat analogue dans le cadre DDFV. Plus généralement, la discrétisation DDFV classique d'une équation elliptique linéaire $-\Delta u = f$ dans le cas de maillages non admissibles est non monotone, dans le sens où en considérant un terme source positif, la solution approchée obtenue peut présenter des valeurs négatives. Dans [CH16], un schéma DDFV monotone non linéaire pour l'équation de diffusion est proposé, mais il présente un défaut de coercivité ne permettant pas de démontrer sa convergence.

À ce stade, l'étude proposée dans le cadre TPFA ne peut donc pas s'adapter au contexte DDFV et l'analyse complète nécessitera la mise en place d'une nouvelle stratégie.
<figure><img src="image_19.png" /></figure>
# DEUXIÈME PARTIE

## SCHÉMAS PRÉSERVANT LA LIMITE DIFFUSIVE DE SYSTÈMES HYPERBOLIQUES AVEC TERMES SOURCES
<figure><img src="image_20.png" /></figure>
Taux de convergence vers la limite diffusive pour une discrétisation du
$p$-système

Ce chapitre est consacré à l'étude de la limite de diffusion d'un schéma numérique pour un système hyperbolique particulier : le p-système avec friction. Ce modèle décrit la dynamique isotherme lagrangienne d'un gaz de covolume $\tau > 0$ et de vitesse $u$ :

$$
(3.1) \qquad \left\{
\begin{array}{l}
\partial_t \tau - \partial_x u = 0, \\
\partial_t u + \partial_x p(\tau) = -\sigma u,
\end{array}
\right.
\qquad (x,t) \in \mathbb{R} \times \mathbb{R}_+.
$$

La loi de pression $p$ est supposée vérifier les hypothèses suivantes, assurant notamment
l'hyperbolicité du système homogène correspondant :

$$
(3.2) \quad \begin{array}{l}
p \in C^2(\mathbb{R}_+^*), \quad p(\tau) > 0, \quad p'(\tau) < 0, \\
\quad \text{si } \tau \ge c > 0 \text{ alors il existe } m \text{ tel que } p(\tau) \ge m > 0 \text{ et } p'(\tau) \le -m < 0.
\end{array}
$$

Typiquement, les fonctions de la forme $p(s) = s^{-\gamma}$, $\gamma \ge 1$ peuvent être considérées.

À ce système est associée une inégalité d'entropie permettant de déterminer les solutions physiquement admissibles :

$$
(3.3) \qquad \partial_t \eta(\tau, u) + \partial_x \psi(\tau, u) \le -\sigma u^2 \le 0,
$$

où le couple entropie-flux $(\eta, \psi)$ est défini par

$$
\eta(\tau, u) = \frac{u^2}{2} - P(\tau), \quad \psi(\tau, u) = u p(\tau),
$$

l'énergie interne $-P(\tau)$ étant une primitive de $-p(\tau)$.
Nous nous intéressons au comportement en temps long des solutions $w = (\tau, u)^T$ de ce système, et plus précisément à leur limite $\bar{w} = (\bar{\tau}, \bar{u})^T$ dans le régime de diffusion, vérifiant le problème parabolique non linéaire suivant (équation des milieux poreux) :

$$ (3.4) \qquad \begin{cases} \partial_t \bar{\tau} + \frac{1}{\sigma} \partial_{xx} p(\bar{\tau}) = 0, \\ \partial_x p(\bar{\tau}) = -\sigma \bar{u}, \end{cases} \quad (x,t) \in \mathbb{R} \times \mathbb{R}_+ $$

Il existe une littérature importante concernant l'étude de la limite de diffusion de systèmes hyperboliques de lois de conservation. Citons par exemple [HL92] et [Nis96] où des taux de convergence en temps sont obtenus par des méthodes d'énergie, avec utilisation d'inégalités de Sobolev. Des taux explicites de convergence en temps ont également été démontrés pour une large classe de systèmes hyperboliques dissipatifs [BHN07]. Ces résultats semblant difficiles à adapter dans un cadre discret, nous considérons plutôt le point de vue adopté par exemple dans [MMS88, Mei10], qui repose sur l'introduction d'un paramètre $\varepsilon$ de mise à l'échelle.

Plus précisément, nous nous intéressons au régime temps long et friction dominante pour le système (3.1). Dans ce cadre de forte friction, il apparaît que la vitesse est de l'ordre de $\varepsilon$. Le système mis à l'échelle que nous considérons s'écrit finalement

$$ (3.5) \qquad \begin{cases} \partial_t \tau^\varepsilon - \partial_x u^\varepsilon = 0, \\ \varepsilon^2 \partial_t u^\varepsilon + \partial_x p(\tau^\varepsilon) = -\sigma u^\varepsilon, \end{cases} \quad (x,t) \in \mathbb{R} \times \mathbb{R}_+ $$

On appelle tout couple de fonctions couple entropie-flux d'entropie muni de l'inégalité d'entropie suivante :

$$ (3.6) \qquad \partial_t \eta^\varepsilon(\tau^\varepsilon, u^\varepsilon) + \partial_x \psi(\tau^\varepsilon, u^\varepsilon) \le -\sigma |u^\varepsilon|^2 \le 0, $$

avec

$$ (3.7) \qquad \eta^\varepsilon(\tau, u) = \varepsilon^2 \frac{u^2}{2} - P(\tau). $$

Les solutions fortes entropiques $w^\varepsilon = (\tau^\varepsilon, u^\varepsilon)^T$ de (3.5) vérifient en particulier l'inégalité d'entropie (3.6) avec $\eta^\varepsilon$ définie par (3.7).

Quand $\varepsilon$ tend vers 0, les solutions $w^\varepsilon = (\tau^\varepsilon, u^\varepsilon)^T$ de (3.5) convergent vers les solutions $\bar{w} = (\bar{\tau}, \bar{u})^T$ de (3.4). Suivant [LT13], ce résultat peut être établi grâce à une méthode d'entropie relative dont je rappelle maintenant les grandes lignes, puisque l'idée est de l'adapter au cadre discret.

L'entropie relative est définie comme étant le terme quadratique du développement de Taylor au voisinage de $\bar{w}$ solution régulière de (3.4) :

$$ (3.8) \qquad \begin{aligned} \eta^\varepsilon(w|\bar{w}) &= \eta^\varepsilon(w) - \eta^\varepsilon(\bar{w}) - \nabla\eta^\varepsilon(\bar{w}) \cdot (w - \bar{w}) \\ &= \frac{\varepsilon^2}{2}(u - \bar{u})^2 - P(\tau|\bar{\tau}), \end{aligned} $$

où $P(\tau|\bar{\tau}) = P(\tau) - P(\bar{\tau}) - p(\bar{\tau})(\tau - \bar{\tau})$.

Cette entropie relative satisfait une loi d'évolution décrite dans le lemme suivant.
**Lemme 3.1.** Soit $w^\epsilon = (\tau^\epsilon, u^\epsilon)^T$ une solution forte entropique de (3.5) et $\bar{w} = (\bar{\tau}, \bar{u})^T$ une solution régulière de (3.4). Alors l'entropie relative vérifie

$$ (3.9) \qquad \partial_t \eta^\epsilon(w^\epsilon|\bar{w}) + \partial_x \psi(w^\epsilon|\bar{w}) = -\sigma(u^\epsilon - \bar{u})^2 + \frac{1}{\sigma} p(\tau^\epsilon|\bar{\tau}) \partial_{xx} p(\bar{\tau}) + \frac{\epsilon^2}{\sigma} (u^\epsilon - \bar{u}) \partial_{xt} p(\bar{\tau}), $$

où le flux relatif est défini par

$$ (3.10) \qquad \psi(w|\bar{w}) = (u - \bar{u})(p(\tau) - p(\bar{\tau})). $$

Introduisons maintenant la quantité suivante, qui va permettre de quantifier l'écart entre $w^\epsilon$ et $\bar{w}$ :

$$ (3.11) \qquad \phi_{\epsilon}(t) = \int_{\mathbb{R}} \eta^{\epsilon}(w^{\epsilon}|\bar{w}) \, dx. $$

En supposant que les conditions initiales des systèmes (3.5) et (3.4) sont telles que les conditions aux limites suivantes sont satisfaites :

$$ (3.12) \qquad \begin{aligned} \lim_{x \to \pm\infty} \tau^{\epsilon}(x,t) &= \lim_{x \to \pm\infty} \bar{\tau}(x,t) = \tau_{\pm} > 0, \\ \lim_{x \to \pm\infty} u^{\epsilon}(x,t) &= \lim_{x \to \pm\infty} \bar{u}(x,t) = 0, \end{aligned} $$

le résultat suivant peut finalement être démontré.

**Théorème 3.1.** Soit $\bar{w}$ une solution régulière de (3.4) telle que $\bar{\tau} \ge c > 0$ sur $Q_T = \mathbb{R} \times [0, T)$. Supposons qu'il existe $K > 0$ telle que

$$ (3.13) \qquad \|\partial_{xx} p(\bar{\tau})\|_{L^{\infty}(Q_T)} \le K, \qquad \|\partial_{xt} p(\bar{\tau})\|_{L^2(Q_T)} \le K. $$

Soit $w^\epsilon$ une solution forte entropique de (3.5) avec $\tau \ge c > 0$ sur $Q_T$, telle que $\phi_\epsilon(0) < +\infty$. Alors

$$ (3.14) \qquad \phi_{\epsilon}(t) \le C(\phi_{\epsilon}(0) + \epsilon^4), \quad t \in [0, T), $$

où C est une constante dépendant de T, $\sigma$ et K. De plus, si $\phi_\epsilon(0) \to 0$ quand $\epsilon \to 0$, alors

$$ (3.15) \qquad \sup_{t \in [0,T]} \phi_{\epsilon}(t) \to 0 \quad \text{quand } \epsilon \to 0. $$

La preuve de ce théorème repose sur l'intégration en espace et en temps de la loi d'évolution de l'entropie relative (3.9), et sur l'utilisation des propriétés de $p(\bar{\tau})$ et ses dérivées, grâce notamment à l'application d'inégalités de Cauchy-Schwarz et Young. La conclusion est obtenue par application du lemme de Grönwall.

Ce résultat (3.14) permet de quantifier la vitesse de convergence vers la limite de diffusion quand $\epsilon$ tend vers 0. Remarquons de plus que du fait de la convexité de l'entropie $\eta^\epsilon$, l'intégrale de l'entropie relative $\phi_\epsilon$ se comporte comme $\|w^\epsilon - \bar{w}\|_{L^2(\mathbb{R})}^2$, permettant de
déduire une convergence en norme $L^2$ à vitesse $\varepsilon^2$.

L'objectif des travaux présentés dans ce chapitre est de retrouver rigoureusement au niveau discret ce taux de convergence, quand $w^\varepsilon$ et $\bar{w}$ sont approchés par des schémas numériques appropriés. En effet, d'un point de vue numérique, il n'est pas immédiat de construire un schéma pour (3.5) restant valide dans la limite $\varepsilon \to 0$. Nous considérons donc une méthode *préservant l'asymptotique (AP)* au sens de Jin [Jin99, JPT98], c'est-à-dire fournissant une discrétisation stable de (3.5) pour toute valeur de $\varepsilon > 0$, et conduisant lorsque $\varepsilon$ tend vers 0, à paramètres de discrétisation fixés, à un schéma consistant avec le problème limite (3.4). Établir l'analogue discret du théorème 3.1 permet en particulier de prouver le caractère AP du schéma numérique considéré.

Dans la section 3.1, je considère le cas semi-discret (continu en temps et discrétisé en espace) que nous avons étudié avec Christophe Berthon et Hélène Mathis [1]. Nous prouvons, pour une discrétisation volumes finis avec flux de type HLL [HLvL83], l'analogue semi-discret du théorème 3.1. Ce résultat est fondé sur un analogue semi-discret de la loi d'évolution de l'entropie relative (3.9), mais nécessite le contrôle de termes additionnels provenant de la viscosité numérique. Dans la section 3.2, je présente l'extension de ces résultats au cadre complètement discret, que nous avons obtenue en collaboration avec Christophe Berthon et Solène Bulteau [15]. L'idée générale de la preuve est la même, mais la discrétisation temporelle fait apparaître beaucoup de termes supplémentaires à contrôler, rendant la preuve assez technique.

## 3.1 Le cadre semi-discret en espace

Considérons un maillage uniforme de l'espace constitué de cellules $(x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}})_{i \in \mathbb{Z}}$ de taille $\Delta x$. La solution de (3.6) est approchée sur chaque cellule $(x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}})$ par une fonction (dépendant du temps) constante en espace $w_i(t) = (\tau_i(t), u_i(t))^T$ (à partir de maintenant, la dépendance en $\varepsilon$ est omise pour alléger les notations).

Définissons maintenant les normes semi-discrètes qui seront utilisées par la suite. Pour $v(t) = (v_i(t))_{i \in \mathbb{Z}}$ une fonction de $t \in [0, T)$, nous notons :

$$
\begin{align}
\|D_x v\|_{L^\infty(Q_T)} &= \sup_{t \in [0,T]} \sup_{i \in \mathbb{Z}} \left| \frac{v_{i+1} - v_i}{\Delta x} \right|, \notag \\
\|\tilde{D}_{xx} v\|_{L^\infty(Q_T)} &= \sup_{t \in [0,T]} \sup_{i \in \mathbb{Z}} \left| \frac{v_{i+2} - 2v_i + v_{i-2}}{(2\Delta x)^2} \right|, \notag \\
\|D_{xx} v\|_{L^\infty(Q_T)} &= \sup_{t \in [0,T]} \sup_{i \in \mathbb{Z}} \left| \frac{v_{i+1} - 2v_i + v_{i-1}}{(\Delta x)^2} \right|, \tag{3.16} \\
\|\tilde{D}_{tx} v\|_{L^2(Q_T)} &= \left( \int_0^T \sum_{i \in \mathbb{Z}} \Delta x \left| \frac{d}{dt} \left( \frac{v_{i+1} - v_{i-1}}{2\Delta x} \right) \right|^2 (s) ds \right)^{1/2}, \notag \\
\|D_{xx} v\|_{L^2(Q_T)} &= \left( \int_0^T \sum_{i \in \mathbb{Z}} \Delta x \left| \frac{v_{i+1} - 2v_i + v_{i-1}}{(\Delta x)^2} \right|^2 (s) ds \right)^{1/2}. \notag
\end{align}
$$
Le schéma numérique semi-discret en espace que nous étudions, basé sur le flux HLL [HLvL83], s'écrit

$$ (3.17) \qquad \begin{cases} \frac{d}{dt} \tau_i = \frac{1}{2\Delta x} (u_{i+1} - u_{i-1}) + \frac{\lambda}{2\Delta x} (\tau_{i+1} - 2\tau_i + \tau_{i-1}), \\ \frac{d}{dt} u_i = \frac{\lambda}{2\Delta x} (u_{i+1} - 2u_i + u_{i-1}) - \frac{1}{2\varepsilon^2 \Delta x} (p(\tau_{i+1}) - p(\tau_{i-1})) - \frac{\sigma}{\varepsilon^2} u_i, \end{cases} $$

avec

$$ (3.18) \qquad \lambda = \sup_{t \in (0,T)} \max_{i \in \mathbb{Z}} (\sqrt{-p'(\tau_i)}). $$

Formellement, en prenant $\varepsilon \to 0$ dans ce schéma (3.17), nous obtenons une discrétisation consistante du problème limite (3.4) :

$$ (3.19) \qquad \begin{cases} \frac{d}{dt} \bar{\tau}_i = \frac{1}{2\Delta x} (\bar{u}_{i+1} - \bar{u}_{i-1}) + \frac{\lambda}{2\Delta x} (\bar{\tau}_{i+1} - 2\bar{\tau}_i + \bar{\tau}_{i-1}), \\ \sigma \bar{u}_i = - \frac{p(\bar{\tau}_{i+1}) - p(\bar{\tau}_{i-1})}{2\Delta x}. \end{cases} $$

En notant $\bar{w}_i(t) = (\bar{\tau}_i(t), \bar{u}_i(t))^T$, l'entropie relative semi-discrète est donnée par

$$ (3.20) \qquad \eta_i^{\varepsilon}(t) = \eta^{\varepsilon}(w_i(t)|\bar{w}_i(t)) = \frac{\varepsilon^2}{2}(u_i(t) - \bar{u}_i(t))^2 - P(\tau_i(t)|\bar{\tau}_i(t)), $$

et l'analogue semi-discret de $\phi_\varepsilon$, noté $\phi_\varepsilon^{\Delta x}$, est donné par

$$ (3.21) \qquad \phi_{\varepsilon}^{\Delta x}(t) = \sum_{i \in \mathbb{Z}} \Delta x \eta_i^{\varepsilon}(t). $$

En supposant que l'analogue de la condition aux limites (3.12) est vérifiée par les solutions approchées :

$$ (3.22) \qquad \begin{aligned} \lim_{i \to \pm\infty} \tau_i &= \lim_{i \to \pm\infty} \bar{\tau}_i = \tau_{\pm}, \\ \lim_{i \to \pm\infty} u_i &= \lim_{i \to \pm\infty} \bar{u}_i = 0, \end{aligned} $$

le résultat suivant peut être établi.

**Théorème 3.2.** *Soit $(\bar{w}_i(t))_i$ solution régulière de (3.19) telle que $\bar{\tau}_i \ge c > 0$ sur $Q_T$. Supposons qu'il existe $K > 0$ telle que*

$$ (3.23) \qquad \|\tilde{D}_{tx}p(\bar{\tau})\|_{L^2(Q_T)} \le K, \qquad \|\tilde{D}_{xx}p(\bar{\tau})\|_{L^\infty(Q_T)} \le K, $$

$$ (3.24) \qquad \|D_{xx}\bar{\tau}\|_{L^\infty(Q_T)} \le K, \qquad \|D_x\bar{\tau}\|_{L^\infty(Q_T)} \le K, \qquad \|D_{xx}\bar{u}\|_{L^2(Q_T)} \le K. $$

*Soit $(w_i(t))_i$ une solution de (3.17) avec $\tau_i \ge c > 0$ sur $Q_T$, telle que $\phi_\varepsilon^{\Delta x}(0) < +\infty$. Alors*

$$ (3.25) \qquad \phi_{\varepsilon}^{\Delta x}(t) \le C(\phi_{\varepsilon}^{\Delta x}(0) + \varepsilon^4), \quad t \in [0, T), $$

où C est une constante dépendant de T, $\sigma$ et K. De plus, si $\phi_\varepsilon^{\Delta x}(0) \to 0$ quand $\varepsilon \to 0$, alors $\sup_{t \in [0,T)} \phi_\varepsilon^{\Delta x}(t) \to 0$ quand $\varepsilon \to 0$.
De manière analogue au cadre continu, cet énoncé fait intervenir des hypothèses sur des normes semi-discrètes de la solution régulière du schéma (3.19) pour le problème parabolique, définies par (3.16). Plus précisément, les hypothèses (3.23) correspondent aux hypothèses (3.13), tandis que les hypothèses supplémentaires (3.24) sont nécessaires pour contrôler les termes additionnels dus à la viscosité numérique.

La preuve de ce théorème repose sur une version semi-discrète de la loi d'évolution de l'entropie relative énoncée dans le lemme suivant.

**Lemme 3.2.** Soit $(w_i)_i = (\tau_i, u_i)_i^T$ une solution de (3.17) et $(\bar{w}_i)_i = (\bar{\tau}_i, \bar{u}_i)_i$ une solution de (3.19). Alors l'entropie relative $\eta_i^\varepsilon$ vérifie

$$ (3.26) \qquad \begin{aligned} \frac{d\eta_i^\varepsilon}{dt} + \frac{1}{\Delta x}(\psi_{i+1/2} - \psi_{i-1/2}) = & -\sigma(u_i - \bar{u}_i)^2 \\ & + \frac{1}{\sigma} \frac{p(\bar{\tau}_{i+2}) - 2p(\bar{\tau}_i) + p(\bar{\tau}_{i-2})}{(2\Delta x)^2} p(\tau_i | \bar{\tau}_i) \\ & + \frac{\varepsilon^2}{\sigma}(u_i - \bar{u}_i) \frac{d}{dt} \left( \frac{p(\bar{\tau}_{i+1}) - p(\bar{\tau}_{i-1})}{2\Delta x} \right) \\ & + R_i^u + R_i^\tau, \end{aligned} $$

où $\psi_{i+1/2}$ correspond à une approximation du flux relatif $\psi$ à l'interface $x_{i+1/2}$ donnée par

$$ (3.27) \qquad \psi_{i+1/2} = \frac{1}{2}(u_i - \bar{u}_i)(p(\tau_{i+1}) - \bar{p}(\bar{\tau}_{i+1})) + \frac{1}{2}(u_{i+1} - \bar{u}_{i+1})(p(\tau_i) - \bar{p}(\bar{\tau}_i)), $$

et les quantités $R_i^u$ et $R_i^\tau$ sont des restes numériques définis par

$$ (3.28) \qquad \begin{aligned} R_i^u &= \frac{\lambda\varepsilon^2}{2\Delta x}(u_i - \bar{u}_i)(u_{i+1} - 2u_i + u_{i-1}), \\ R_i^\tau &= -\frac{\lambda}{2\Delta x}((p(\tau_i) - p(\bar{\tau}_i))(\tau_{i+1} - 2\tau_i + \tau_{i-1}) - (\tau_i - \bar{\tau}_i)p'(\bar{\tau}_i)(\bar{\tau}_{i+1} - 2\bar{\tau}_i + \bar{\tau}_{i-1})). \end{aligned} $$

Comme dans le cas continu, le résultat du théorème 3.2 se démontre en intégrant en temps l'identité (3.26). Les deux termes résiduels $R_i^u$ et $R_i^\tau$ sont alors contrôlés en utilisant les hypothèses additionnelles (3.24), et l'estimation finale est obtenue en utilisant les propriétés sur $p(\bar{\tau}_i)$ et ses dérivées, des inégalités de Cauchy-Schwarz et Young, et finalement le lemme de Grönwall.

## 3.2 Le cadre complètement discret

Pour obtenir un résultat complet, il est nécessaire de considérer également une discrétisation temporelle. Pour cela, notons $\Delta t > 0$ le pas de temps, et $t^n = n \Delta t$ pour $n = 0, \dots, N_T + 1$, avec $T = t^{N_T+1}$ le temps final.

Nous considérons un schéma analogue à celui proposé dans [JPT98], basé sur la reformulation suivante du système (3.5) :

$$ \left\{ \begin{aligned} \partial_t \tau - \partial_x u &= 0, \\ \partial_t u + \partial_x p(\tau) &= -\frac{1}{\varepsilon^2} \left( \sigma u + (1 - \varepsilon^2) \partial_x p(\tau) \right). \end{aligned} \right. $$
La discrétisation est fondée sur un splitting en deux étapes. Dans un premier temps, le système hyperbolique homogène (non raide)

$$
\begin{cases}
\partial_t \tau - \partial_x u = 0, \\
\partial_t u + \partial_x p(\tau) = 0.
\end{cases}
$$

est traité par un schéma HLL explicite :

$$
(3.29a) \qquad \tau_i^{n+\frac{1}{2}} = \tau_i^n - \frac{\Delta t}{\Delta x} \left( \mathcal{F}_{i+\frac{1}{2}}^{\tau} - \mathcal{F}_{i-\frac{1}{2}}^{\tau} \right),
$$

$$
(3.29b) \qquad u_i^{n+\frac{1}{2}} = u_i^n - \frac{\Delta t}{\Delta x} \left( \mathcal{F}_{i+\frac{1}{2}}^u - \mathcal{F}_{i-\frac{1}{2}}^u \right),
$$

les flux numériques étant donnés par

$$
\begin{align*}
\mathcal{F}_{i+\frac{1}{2}}^{\tau} &= \frac{1}{2}(-u_i^n - u_{i+1}^n) - \frac{\lambda}{2}(\tau_{i+1}^n - \tau_i^n), \\
\mathcal{F}_{i+\frac{1}{2}}^{u} &= \frac{1}{2}(p(\tau_i^n) + p(\tau_{i+1}^n)) - \frac{\lambda}{2}(u_{i+1}^n - u_i^n),
\end{align*}
$$

où

$$
(3.30) \qquad \lambda = \sup_{0 \le n \le N_T} \max_{i \in \mathbb{Z}} (\sqrt{-p'(\tau_i)}).
$$

Dans un deuxième temps, le terme source raide est pris en compte de manière implicite :

$$
\tau_i^{n+1} = \tau_i^{n+\frac{1}{2}},
$$
$$
\frac{u_i^{n+1} - u_i^{n+\frac{1}{2}}}{\Delta t} = -\frac{1}{\varepsilon^2} \left( \sigma u_i^{n+1} + (1-\varepsilon^2) \frac{p_{i+\frac{1}{2}}^{n+1} - p_{i-\frac{1}{2}}^{n+1}}{\Delta x} \right),
$$

où $p_{i+1/2}$ sont des discrétisations centrées de $p(\tau)$.

Puisque $\tau_i^{n+1} = \tau_i^{n+\frac{1}{2}}$, $u_i^{n+1}$ peut finalement être calculé explicitement à partir de
$(\tau_i^n, u_i^n)_i$, et l'étape de relaxation s'écrit finalement

$$
(3.31a) \qquad \tau_i^{n+1} = \tau_i^{n+\frac{1}{2}},
$$

$$
(3.31b) \quad u_i^{n+1} = \left( \frac{\varepsilon^2}{\varepsilon^2 + \sigma \Delta t} \right) u_i^{n+\frac{1}{2}} - \Delta t \left( \frac{1-\varepsilon^2}{\Delta t \sigma + \varepsilon^2} \right) \frac{p(\tau_{i+1}^{n+\frac{1}{2}}) - p(\tau_{i-1}^{n+\frac{1}{2}})}{2 \Delta x}.
$$

L'approximation du problème limite (3.4) est obtenue en passant à la limite $\varepsilon \to 0$
dans le schéma (3.29)–(3.31). Nous obtenons

$$
(3.32a) \quad \bar{\tau}_i^{n+1} = \bar{\tau}_i^n + \frac{\Delta t}{2 \Delta x} (\bar{u}_{i+1}^n - \bar{u}_{i-1}^n) - \frac{\lambda \Delta t}{2 \Delta x} (\bar{\tau}_{i+1}^n - 2\bar{\tau}_i^n + \bar{\tau}_{i-1}^n),
$$

$$
(3.32b) \quad \bar{u}_i^{n+1} = - \frac{1}{2 \sigma \Delta x} \left( p(\bar{\tau}_{i+1}^{n+1}) - p(\bar{\tau}_{i-1}^{n+1}) \right),
$$

qui est bien consistant avec le problème parabolique limite puisque le dernier terme de la
première équation est un terme de viscosité numérique.

Le résultat principal est énoncé dans le théorème suivant.
**Théorème 3.3.** Soient $(w_i^n)_{i,n}$ la solution du schéma (3.29)-(3.31) et $(\bar{w}_i^n)_{i,n}$ la solution du schéma limite (3.32). Supposons les hypothèses suivantes satisfaites :

* il existe $K > 0$ tel que

$$ (3.33) \qquad \|D_t \bar{u}^{n+\frac{1}{2}}\|_{L^2}, \quad \|D_{xx} \bar{u}^n\|_{L^2}, \quad \|D_t \bar{\tau}^{n+\frac{1}{2}}\|_{L^\infty}, \quad \|D_x \bar{\tau}^n\|_{L^\infty} \le K, $$

* il existe $L_\tau > 0$ tel que

$$ (3.34) \qquad \frac{1}{L_\tau} \le \tau_i^n, \bar{\tau}_i^n \le L_\tau, $$

* la fonction de pression p vérifie, pour tous $\tau \in [1/L_\tau, L_\tau]$,

$$ (3.35) \qquad \frac{1}{L_p} \le p(\tau), p''(\tau) \le L_p, \quad -L_p \le p'(\tau), p^{(3)}(\tau) \le -\frac{1}{L_p}, $$

* les paramètres de discrétisation vérifient une condition CFL parabolique $\Delta t/\Delta x^2 \le C_p, C_p$ constante dépendant de $K, L_\tau$ et $L_p$,

* $\varepsilon$ est tel que

$$ (3.36) \qquad \varepsilon^2 \le \min \left( \frac{\sigma}{C_p(2 + 15 L_p^2)}, \frac{\sigma}{8\lambda} \Delta x \right). $$

Alors, l'analogue discret de $\phi_\varepsilon$ vérifie

$$ (3.37) \qquad \phi_\varepsilon^{N_T+1} \le C(\phi_\varepsilon^0 + \|u^0 - \bar{u}^0\|_{L^2}^2 + \varepsilon^4), $$

où $C$ est une constante dépendant uniquement de $T$ et des paramètres $\sigma, \lambda, K, L_\tau$ et $L_p$.

Comme dans les cas continu et semi-discret, cet énoncé fait intervenir des hypothèses (3.33) sur l'approximation $\bar{w}_i^n$ de la solution du problème limite, qui permettent de contrôler les termes de viscosité additionnels dus à la discrétisation en espace et en temps. Dans le cas continu et semi-discret, des conditions sont imposées uniquement sur les deux premières dérivées de la pression $p$. Ici, des hypothèses jusqu'à la dérivée d'ordre 3 (3.35) sont requises pour contrôler les termes de reste additionnels issus de la discrétisation temporelle. Notons que ces hypothèses sont vérifiées par la fonction de pression classique $p(\tau) = \tau^{-\gamma}$, $\gamma > 1$. L'hypothèse (3.36) n'est pas restrictive en pratique puisque nous sommes intéressés par la limite $\varepsilon \to 0$. Enfin, la condition CFL parabolique est naturelle, le schéma limite étant une discrétisation explicite en temps d'un problème diffusif.

La preuve du théorème 3.3 suit le même cheminement que celles des cas continu et semi-discret. Une version discrète de la loi d'entropie relative (3.9) est tout d'abord établie. L'idée générale pour contrôler les nombreux termes de reste est de faire apparaître soit des termes en $\varepsilon^4$, soit des termes impliquant uniquement la condition initiale, soit des termes qui vont pouvoir être compensés par les termes négatifs déjà présents.
<figure><img src="image_21.png" /><figcaption>FIGURE 3.1 : Intégrale en espace de l'entropie relative $\phi_{\epsilon}^{N_T+1}$ en fonction de $\epsilon$, en échelle logarithmique.</figcaption></figure>

Nous constatons que le taux de convergence en $\epsilon^4$ établi dans les cadres continu [LT13] et semi-discret est bien retrouvé pour le schéma (3.29)–(3.31). Sur la figure 3.1 est représentée l’intégrale spatiale de l’entropie relative discrète $\phi_{\epsilon}^{N_T+1}$ en fonction de $\epsilon$ en échelle logarithmique, pour deux données initiales, et pour différents nombres $N$ de cellules en espace. Les simulations numériques donnent un taux de convergence en parfait accord avec celui établi dans le théorème 3.3, ce qui laisse à penser que ce résultat est optimal.

## 3.3 Perspectives

Dans l'article [LT13], la méthode d'entropie relative est également appliquée au système d'Euler isentropique et à un système quasi-linéaire de viscoélasticité. Les preuves sont fondées sur les expressions des entropies particulières à chaque système. Cependant, la méthode n'est pas formalisée pour un système général, comme c'est le cas pour l'asymptotique hyperbolique [Tza05]. Pourtant, la structure de relaxation est définie de façon générale dans [BHN07, BLT13]. Ainsi, dans [BLT13], un système hyperbolique non linéaire général avec terme de relaxation est considéré :

$$ \epsilon \partial_t U + \partial_x F(U) = - \frac{R(U)}{\epsilon}. $$

À l'aide d'un développement de Chapman-Enskog, un système effectif d'équations décrivant la limite singulière de diffusion $\epsilon \to 0$ est dérivé. La structure de ce système est discutée, de même que le rôle de l'entropie. Ainsi, il semble qu'en établissant une version discrète de ce formalisme, et donc une structure générale de schémas volumes finis pré servant l'asymptotique, nous pourrions exploiter la méthode d'entropie relative pour étudier la convergence de tels schémas.

Un point de départ pourrait être d'adapter au niveau discret des travaux récents [Bia19, Bia21] dans lesquels la méthode d'entropie relative est utilisée pour étudier la
limite diffusive de systèmes de type BGK à vitesses discrètes.
Schémas consistants à la limite diffusive pour les équations de
Saint-Venant avec friction de Manning

Dans ce chapitre, je considère des discrétisations des équations de Saint-Venant avec fric-
tion de Manning [Man90] consistantes à la limite de diffusion. Ce sytème, qui décrit l’évo-
lution de la hauteur d’eau h et du débit q, s’écrit

$$
(4.1) \qquad \begin{cases} \partial_t h + \partial_x q = 0, \\ \partial_t q + \partial_x \left( \frac{q^2}{h} + \frac{gh^2}{2} \right) = -k |q| q h^{-\eta}, \end{cases} \quad (x,t) \in \mathbb{R} \times \mathbb{R}_+,
$$

g étant la constante universelle de gravitation, et les paramètres de friction apparaissant
dans le terme source de Manning étant η, en général égal à 7/3, et k qui définit l'intensité
de la friction.

Je m’intéresse, comme dans le chapitre précédent, au comportement des solutions de
ce problème en temps long et friction dominante. Pour cela, une mise à l’échelle en temps
et en friction est considérée, en introduisant un paramètre ε > 0 qui sera amené à tendre
vers 0 :

$$
(4.2) \qquad t \leftarrow t/\epsilon, \quad k \leftarrow k/\epsilon^2.
$$

La mise à l'échelle de la friction (en ε²) est différente de celle proposée au chapitre 3 (en
ε), du fait du terme quadratique en q dans le terme source [BLT13]. Après mise à l'échelle,
le système (4.1) se réécrit :

$$
(4.3) \qquad \begin{cases} \varepsilon \partial_t h + \partial_x q = 0, \\ \varepsilon \partial_t q + \partial_x \left( \frac{q^2}{h} + \frac{gh^2}{2} \right) = -\frac{k}{\varepsilon^2} |q| q h^{-\eta}. \end{cases}
$$
Le comportement quand $\varepsilon$ tend vers 0 de $h$ et $q$ peut être étudié formellement en
utilisant un développement de Chapman-Enskog :

$$h = h^0 + \varepsilon h^1 + \dots, \quad q = q^0 + \varepsilon q^1 + \dots$$

La hauteur d'eau $h$ étant ici considérée strictement positive, avec de plus $h^0 > 0$, il est
clair que $q^0$ est nécessairement nul, et finalement le problème limite est donné par

$$(4.4) \qquad \begin{cases} \partial_t h^0 + \partial_x q^1 = 0, \\ \partial_x \left( \frac{g(h^0)^2}{2} \right) = -k |q^1| q^1 (h^0)^{-\eta}. \end{cases}$$

La deuxième équation de ce système conduit à l'expression suivante pour $q^1$ :

$$(4.5) \qquad q^1 = -\operatorname{sign}(\partial_x h^0) \sqrt{\frac{(h^0)^\eta}{k} \left| \partial_x \left( \frac{g(h^0)^2}{2} \right) \right|},$$

ce qui mène finalement à l'équation diffusive non linéaire suivante (de type p-laplacien)
pour $h^0$ :

$$(4.6) \qquad \partial_t h^0 + \partial_x \left( -\operatorname{sign}(\partial_x h^0) \sqrt{\frac{(h^0)^\eta}{k} \left| \partial_x \left( \frac{g(h^0)^2}{2} \right) \right|} \right) = 0.$$

L'objectif des travaux présentés dans ce chapitre, effectués dans le cadre de la thèse
de Solène Bulteau, est de construire des discrétisations de (4.1) asymptotiquement consis-
tantes (AC) avec la limite de diffusion, c'est-à-dire consistantes avec le modèle hyperbo-
lique (4.1) pour tous $\varepsilon > 0$, et convergeant, quand $\varepsilon$ tend vers 0, vers un schéma consistant
avec le modèle limite de diffusion (4.5)-(4.6). Cette propriété de consistance asymptotique
est moins forte que la propriété AP étudiée dans le chapitre 3 puisqu'elle n'impose pas la
stabilité de la méthode indépendamment de la valeur de $\varepsilon > 0$.

Je présente dans ce chapitre deux schémas ayant cette propriété de consistance asymp-
totique. Ces deux discrétisations sont de type Godunov à deux états intermédiaires. Le
cadre général considéré est donné dans la section 4.1.

Dans la section 4.2, je présente une généralisation de la méthode de perturbation
HLL introduite dans [BT11] que nous avons proposée avec Christophe Berthon et Solène
Bulteau [14]. L'idée de cette méthode est d'inclure la discrétisation du terme source dans
le solveur de Riemann approché (de type HLL [HLvL83]). En introduisant dans ce schéma
HLL perturbé un paramètre supplémentaire, la viscosité numérique peut être contrôlée
de sorte à être consistante avec la diffusion du problème limite considéré. Dans l'article
[BT11], le cadre traité est linéaire, dans le sens où le terme source est linéaire et conduit
à une limite de diffusion impliquant un opérateur linéaire (de type laplacien). Dans le cas
qui m'intéresse ici, le système de Saint-Venant avec friction de Manning (4.1), le terme
source est quadratique, conduisant à une limite de diffusion impliquant un opérateur non
linéaire de type p-laplacien. Nous proposons donc une extension adaptée à ce contexte pour
obtenir une discrétisation consistante à la limite de diffusion. Cette consistance est établie formellement à l'aide de développements de Chapman-Enskog de la solution approchée.

Dans la section 4.3, je m'intéresse au schéma *équilibre* de type Godunov proposé dans [MDBCF17], construit initialement dans le but de préserver tous les états stationnaires du système de Saint-Venant avec friction de Manning (4.1). Ces états stationnaires sont gouvernés par un débit uniforme $q_0$, et l'équation suivante sur la hauteur d'eau

$$ (4.7) \qquad \partial_x \left( \frac{q_0^2}{h} + \frac{gh^2}{2} \right) = -kq_0|q_0|h^{-\eta}, $$

qui peut être intégrée et donne une équation non linéaire sur $h$ [MDBCF17]. Avec Solène Bulteau, Mehdi Badsi et Christophe Berthon [13], nous montrons que, à une petite modification près, le schéma *équilibre* proposé est en fait également consistant à la limite de diffusion, tout en continuant de préserver tous les états stationnaires.

## 4.1 Schémas de type Godunov à deux états intermédiaires

Considérons un maillage uniforme en espace, composé de cellules $K_i = (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}})$, $i \in \mathbb{Z}$, de longueur $\Delta x$ et de centre $x_i$. Concernant la discrétisation temporelle, notons $\Delta t$ le pas de temps et $t^n = n\Delta t$ pour $n \in \mathbb{N}$.

Le système (4.1) peut se réécrire sous la forme générique d'un système hyperbolique avec terme source

$$ (4.8) \qquad \partial_t W + \partial_x F(W) = S(W), $$

avec

$$ W = (h, q)^T, \quad F(W) = \left( q, \frac{q^2}{h} + g \frac{h^2}{2} \right)^T, \quad S(W) = (0, -k q |q| h^{-\eta})^T. $$

Notons $W_i^n$ une approximation de la solution $W$ de (4.8) sur la cellule $K_i$, au temps $t^n$. Nous considérons des schémas de type Godunov; au temps $t^{n+1}$, l'approximation est mise à jour comme

$$ (4.9) \qquad W_i^{n+1} = \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} W^\Delta(x, t^n + \Delta t) dx, $$

où

$$ W^\Delta(x, t^n + t) = \tilde{W} \left( \frac{x - x_{i+\frac{1}{2}}}{t}; W_i^n, W_{i+1}^n \right), \quad \text{si } x \in (x_i, x_{i+1}). $$

Les solveurs de Riemann approchés $\tilde{W}$ utilisés dans ce chapitre sont composés de trois ondes discontinues de vitesses $\lambda_L < 0 < \lambda_R$ séparant quatre états constants. Je suppose
dans la suite que $-\lambda_L = \lambda_R = \lambda_{LR} > 0$ pour alléger les notations. Considérons donc

$$
(4.10) \quad \tilde{W} \left( \frac{x}{t}, W_L, W_R \right) = \begin{cases} W_L & \text{si } \frac{x}{t} \le -\lambda_{LR}, \\ W^{L*} & \text{si } -\lambda_{LR} < \frac{x}{t} \le 0, \\ W^{R*} & \text{si } 0 < \frac{x}{t} \le \lambda_{LR}, \\ W_R & \text{si } \lambda_{LR} < \frac{x}{t}. \end{cases}
$$

Les deux états intermédiaires $W^{L*}$ et $W^{R*}$ sont déterminés de telle sorte que la relation
de consistance intégrale suivante soit vérifiée [BC16] :

$$
\begin{equation}
\begin{split}
\frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \tilde{W} \left( \frac{x}{\Delta t}; W_L, W_R \right) dx = \frac{1}{2} (W_L + W_R) - \frac{\Delta t}{\Delta x} (F(W_R) - F(W_L)) \\
\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad + \Delta t \bar{S}(W_L, W_R),
\end{split}
\tag{4.11}
\end{equation}
$$

$\bar{S}(W_L, W_R)$ étant une approximation du terme source. Cette relation (4.11) est valide tant que les vitesses d'ondes exactes de (4.1) appartiennent à $(-\lambda_{LR}, \lambda_{LR})$, conduisant à prendre

$$
(4.12) \qquad \lambda_{LR} = \max \left( \frac{|q_L|}{h_L} + \sqrt{gh_L}, \frac{|q_R|}{h_R} + \sqrt{gh_R} \right),
$$

et à imposer la condition CFL suivante :

$$
(4.13) \qquad \frac{\Delta t}{\Delta x} \max_{i \in \mathbb{Z}} \lambda_{i+\frac{1}{2}} \le \frac{1}{2}.
$$

En utilisant la relation de consistance (4.11) ainsi que la définition du solveur de Riemann
approché, le schéma (4.9) s’écrit finalement

$$
(4.14) \quad W_i^{n+1} = W_i^n + \frac{\Delta t}{\Delta x} \left( \lambda_{i+\frac{1}{2}} (W_{i+\frac{1}{2}}^{L*} - W_i^n) + \lambda_{i-\frac{1}{2}} (W_{i-\frac{1}{2}}^{R*} - W_i^n) \right).
$$

Il reste donc à définir de manière adéquate les deux états intermédiaires $W^{L*}$ et $W^{R*}$. Les deux choix exposés dans les sections 4.2 et 4.3 font intervenir l'état intermédiaire $W^*$ du schéma HLL classique [HLvL83] (ou simplement Rusanov ici puisque $-\lambda_L = \lambda_R = \lambda_{LR} > 0$) pour un système hyperbolique homogène :

$$
(4.15) \qquad W^* = \frac{1}{2} (W_R + W_L) - \frac{1}{2\lambda_{LR}} (F(W_R) - F(W_L)).
$$

4.2 Un schéma basé sur la généralisation du cas linéaire

Je m'intéresse tout d'abord à une généralisation au cadre non linéaire du schéma proposé
dans [BT11]. Pour rappel, ce schéma est fondé sur l'écriture de terme source sous la forme

$$
S(W) = k(R(W) - W) = (k + \bar{k})(\bar{R}(W) - W),
$$
le paramètre $\bar{k}$ étant un degré de liberté permettant d'ajuster la diffusion numérique de sorte à obtenir une discrétisation consistante à la limite de diffusion. Dans le cadre linéaire, les états intermédiaires sont définis par

$$ (4.16) \qquad W^{L*} = \alpha_{LR} W^* + (1 - \alpha_{LR}) R(W_L), $$

$$ (4.17) \qquad W^{R*} = \alpha_{LR} W^* + (1 - \alpha_{LR}) R(W_R), $$

où $W^*$ est l'état intermédiaire HLL classique (4.15) et le paramètre $\alpha_{LR}$ est donné par

$$ (4.18) \qquad \alpha_{LR} = \frac{2\lambda_{LR}}{2\lambda_{LR} + (k + \bar{k}_{LR})\Delta x}. $$

Ce schéma est adopté dans [DMTB15] pour les équations de Saint-Venant avec friction de Manning, dans lequel le régime limite considéré est obtenu par la mise à l'échelle suivante (qui n'est pas consistante avec celle du cadre continu (4.2)) :

$$ (4.19) \qquad \Delta t \leftarrow \Delta t/\epsilon, \quad k \leftarrow k/\epsilon, \quad \bar{k}_{LR} \leftarrow \bar{k}_{LR}/\epsilon. $$

En appliquant aux équations de Saint-Venant le schéma (4.14) avec les états intermédiaires (4.16) et (4.17), et en considérant cette mise à l'échelle (4.19), la discrétisation obtenue suivant [DMTB15] s'écrit

$$ (4.20) \qquad h_i^{n+1} = h_i^n - \frac{\Delta t}{2\epsilon\Delta x} \left( \alpha_{i+\frac{1}{2}}^{\epsilon} (q_{i+1}^n - q_i^n) + \alpha_{i-\frac{1}{2}}^{\epsilon} (q_i^n - q_{i-1}^n) \right) \\ \qquad \qquad + \frac{\Delta t}{2\epsilon\Delta x} \left( \lambda_{i+\frac{1}{2}} \alpha_{i+\frac{1}{2}}^{\epsilon} (h_{i+1}^n - h_i^n) - \lambda_{i-\frac{1}{2}} \alpha_{i-\frac{1}{2}}^{\epsilon} (h_i^n - h_{i-1}^n) \right), $$

$$ (4.21) \qquad q_i^{n+1} = q_i^n - \frac{\Delta t}{2\epsilon\Delta x} \left( \alpha_{i+\frac{1}{2}}^{\epsilon} \left( \frac{(q_{i+1}^n)^2}{h_{i+1}^n} + g \frac{(h_{i+1}^n)^2}{2} - \frac{(q_i^n)^2}{h_i^n} - g \frac{(h_i^n)^2}{2} \right) \right. \\ \qquad \qquad \left. + \alpha_{i-\frac{1}{2}}^{\epsilon} \left( \frac{(q_i^n)^2}{h_i^n} + g \frac{(h_i^n)^2}{2} - \frac{(q_{i-1}^n)^2}{h_{i-1}^n} - g \frac{(h_{i-1}^n)^2}{2} \right) \right) \\ \qquad \qquad + \frac{\Delta t}{2\epsilon\Delta x} \left( \lambda_{i+\frac{1}{2}} \alpha_{i+\frac{1}{2}}^{\epsilon} (q_{i+1}^n - q_i^n) - \lambda_{i-\frac{1}{2}} \alpha_{i-\frac{1}{2}}^{\epsilon} (q_i^n - q_{i-1}^n) \right) \\ \qquad \qquad \left. - \frac{k\Delta t}{2\epsilon^2} \left( \alpha_{i+\frac{1}{2}}^{\epsilon} + \alpha_{i-\frac{1}{2}}^{\epsilon} \right) q_i^n |q_i^n|(h_i^n)^{-\eta}, \right) $$

où le paramètre $\alpha$ défini par (4.18) est lui aussi mis à l'échelle :

$$ (4.22) \qquad \alpha_{i+\frac{1}{2}}^{\epsilon} = \frac{2\lambda_{i+\frac{1}{2}}\epsilon}{2\lambda_{i+\frac{1}{2}}\epsilon + (k + \bar{k}_{i+\frac{1}{2}})\Delta x}. $$

Le paramètre $\bar{k}_{i+\frac{1}{2}}$ est déterminé de telle sorte que lorsque $\epsilon \to 0$, le schéma (4.20)–(4.21) fournisse une approximation consistante de l'équation de diffusion (4.6) :

$$ (4.23) \qquad \bar{k}_{i+\frac{1}{2}} = 2\lambda_{i+\frac{1}{2}}^2 \sqrt{\frac{k|h_{i+1}^n - h_i^n|}{g\Delta x((h_{i+1}^n)^\eta + (h_i^n)^\eta)(h_{i+1}^n + h_i^n)}} - k. $$

Le théorème suivant établit l'échec de ce premier schéma (dans le cas d'une hauteur d'eau non constante).
<figure><img src="image_22.png" /><figcaption>FIGURE 4.1 : Comparaison de la solution approchée obtenue par le schéma (4.20)–(4.22) pour différentes valeurs de ε avec la solution du problème limite (4.5)–(4.6).</figcaption></figure>

**Théorème 4.1.** En développant $h_i^n = h_i^{0,n} + O(\epsilon)$, supposons $h_i^{n,0} > 0$ pour tous $i \in \mathbb{Z}$.
Soit $K_i$ une cellule telle que $h_i^{n,0} \neq h_{i-1}^{n,0}$ et $h_i^{n,0} \neq h_{i+1}^{n,0}$. Alors, quand $\epsilon \to 0$, la première
équation (4.20) du schéma s'écrit

$$
\begin{equation}
\begin{split}
\frac{h_i^{0,n+1} - h_i^{0,n}}{\Delta t} ={}& \frac{1}{\Delta x} \left( \operatorname{sign}(h_{i+1}^{0,n} - h_i^{0,n}) \sqrt{\frac{(h_{i+1}^{0,n})^\eta + (h_i^{0,n})^\eta}{2k} \left| \frac{g}{2}(h_{i+1}^{0,n})^2 - \frac{g}{2}(h_i^{0,n})^2 \right|} \right. \\
& \left. - \operatorname{sign}(h_i^{0,n} - h_{i-1}^{0,n}) \sqrt{\frac{(h_i^{0,n})^\eta + (h_{i-1}^{0,n})^\eta}{2k} \left| \frac{g}{2}(h_i^{0,n})^2 - \frac{g}{2}(h_{i-1}^{0,n})^2 \right|} \right),
\end{split}
\tag{4.24}
\end{equation}
$$

tandis que la seconde équation (4.21) se réduit à une discrétisation de $\partial_x h^0 = 0$.

Comme attendu par construction de $\bar{k}_{i+\frac{1}{2}}$, l'équation limite discrète (4.24) est une discrétisation consistante de l'équation de diffusion (4.6). Cependant, la seconde équation (4.21) du schéma développé dans [DMTB15] ne tend pas vers une discrétisation de la relation (4.5), et donne à la limite $\epsilon \to 0$ une hauteur d'eau constante, ce qui est en contradiction avec le cas spécifique considéré ici. Ce comportement non consistant à la limite de diffusion du schéma proposé dans [DMTB15] est illustré sur la figure 4.1.

Partant de ce constat, nous proposons une méthode générique pour construire des schémas consistants à la limite de diffusion pour des problèmes de type (4.8), dont le régime de diffusion est gouverné par une mise à l'échelle de la forme

$t \leftarrow t/\epsilon, \quad k \leftarrow k/\epsilon^m.$

Il s'agit donc de généraliser le cas $m = 1$ traité dans [BT11]. Pour cela, l'idée est d'intro-
duire, en plus de $\bar{k}$, un second degré de liberté $\delta > 0$ :

$$
\partial_t W + \partial_x F(W) = \delta(k + \bar{k})(\bar{R}(W) - W).
$$
En remplaçant $k + \bar{k}_{LR}$ par $\delta(k + \bar{k}_{LR})$, et en procédant à la mise à l'échelle

$$ \Delta t \leftarrow \Delta t / \varepsilon, \quad k \leftarrow k / \varepsilon^m, \quad \bar{k}_{LR} \leftarrow \bar{k}_{LR} / \varepsilon^m, $$

le paramètre $\delta$ est alors fixé égal à $\varepsilon^{m-1}$.

En appliquant cette procédure pour le cas qui nous intéresse, à savoir les équations de Saint-Venant avec friction de Manning pour lequel $m = 2$, nous obtenons le schéma suivant :

$$ (4.25) \qquad h_i^{n+1} = h_i^n - \frac{\Delta t}{2\varepsilon\Delta x} \left( \alpha_{i+\frac{1}{2}}^{\varepsilon} (q_{i+1}^n - q_i^n) + \alpha_{i-\frac{1}{2}}^{\varepsilon} (q_i^n - q_{i-1}^n) \right) \\ \qquad \qquad + \frac{\Delta t}{2\varepsilon\Delta x} \left( \lambda_{i+\frac{1}{2}} \alpha_{i+\frac{1}{2}}^{\varepsilon} (h_{i+1}^n - h_i^n) - \lambda_{i-\frac{1}{2}} \alpha_{i-\frac{1}{2}}^{\varepsilon} (h_i^n - h_{i-1}^n) \right), $$

$$ (4.26) \qquad q_i^{n+1} = q_i^n - \frac{\Delta t}{2\varepsilon\Delta x} \left( \alpha_{i+\frac{1}{2}}^{\varepsilon} \left( \frac{(q_{i+1}^n)^2}{h_{i+1}^n} + g \frac{(h_{i+1}^n)^2}{2} - \frac{(q_i^n)^2}{h_i^n} - g \frac{(h_i^n)^2}{2} \right) \right. \\ \qquad \qquad \left. + \alpha_{i-\frac{1}{2}}^{\varepsilon} \left( \frac{(q_i^n)^2}{h_i^n} + g \frac{(h_i^n)^2}{2} - \frac{(q_{i-1}^n)^2}{h_{i-1}^n} - g \frac{(h_{i-1}^n)^2}{2} \right) \right) \\ \qquad \qquad + \frac{\Delta t}{2\varepsilon\Delta x} \left( \lambda_{i+\frac{1}{2}} \alpha_{i+\frac{1}{2}}^{\varepsilon} (q_{i+1}^n - q_i^n) - \lambda_{i-\frac{1}{2}} \alpha_{i-\frac{1}{2}}^{\varepsilon} (q_i^n - q_{i-1}^n) \right) \\ \qquad \qquad \left. - \frac{k\Delta t}{2\varepsilon^3} \left( \alpha_{i+\frac{1}{2}}^{\varepsilon} + \alpha_{i-\frac{1}{2}}^{\varepsilon} \right) q_i^n |q_i^n|(h_i^n)^{-\eta}, \right) $$

où $\alpha_{i+\frac{1}{2}}^{\varepsilon}$ est toujours défini par (4.22) et $\bar{k}_{i+\frac{1}{2}}$ par (4.23).

Remarquons que nous retrouvons en fait la discrétisation proposée dans [BLT13] pour le cas non linéaire. La différence entre ce schéma (4.25)-(4.26) et le schéma (4.20)-(4.21) développé dans [DMTB15] réside dans l'introduction du coefficient $\delta$ qui modifie l'ordre en $\varepsilon$ de la discrétisation du terme source (ordre $1/\varepsilon^3$ dans (4.26) contre $1/\varepsilon^2$ dans (4.21)). Cette modification permet d'avoir un schéma asymptotiquement consistant à la limite de diffusion.

**Théorème 4.2.** Le schéma (4.25)-(4.26) vérifie les propriétés suivantes à la limite de diffusion $\varepsilon \to 0$ :

* l'équation (4.25) coïncide avec l'équation de diffusion discrète (4.24) qui est une discrétisation de l'équation de diffusion non linéaire (4.6),

* l'équation (4.26) coïncide avec une discrétisation de l'expression (4.5) de $q^1$.

Ce résultat est établi de manière formelle, à l'aide de développements de Chapman-Enskog des solutions approchées $h_i^n$ et $q_i^n$, et l'identification des termes du même ordre en $\varepsilon$. La figure 4.2 présente les résultats obtenus avec ce schéma corrigé pour le même cas test que celui présenté dans la figure 4.1.
<figure><img src="image_23.png" /><figcaption>FIGURE 4.2 : Comparaison de la solution approchée obtenue par le schéma (4.25)–(4.26) pour différentes valeurs de $\epsilon$ avec la solution du problème limite (4.5)–(4.6).</figcaption></figure>

## 4.3 Un schéma construit pour préserver les équilibres

Je m’intéresse maintenant au schéma équilibre proposé dans [MDBCF17]. Il s’écrit sous la forme (4.14), et les états intermédiaires sont définis par :

$$ (4.27) \qquad \begin{aligned} h^{L*} &= h^* + \frac{k\bar{q}|\bar{q}|\overline{h^{-\eta}}\Delta x}{2\gamma_{LR}}, \\ h^{R*} &= h^* - \frac{k\bar{q}|\bar{q}|\overline{h^{-\eta}}\Delta x}{2\gamma_{LR}}, \\ q^{L*} &= q^{R*} = q^* - \frac{k\bar{q}|\bar{q}|\overline{h^{-\eta}}\Delta x}{2\lambda_{LR}} =: q^{LR*}, \end{aligned} $$

où le paramètre $\gamma_{LR}$ est donné par

$$ (4.28) \qquad \gamma_{LR} = -\frac{\bar{q}^2}{h_L h_R} + \frac{g}{2}(h_R + h_L), $$

et supposé strictement positif.

Ce choix d’états intermédiaires vérifie la condition de consistance intégrale (4.11) quand l’approximation du terme source est donnée par

$$ \bar{S} = \begin{pmatrix} 0 \\ -k\bar{q}|\bar{q}|\overline{h^{-\eta}} \end{pmatrix}. $$

Les approximations $\bar{q}$ et $\overline{h^{-\eta}}$ sont alors définies de telle sorte que tous les états stationnaires soient préservés. Pour cela,

$$ (4.29) \qquad \overline{h^{-\eta}} = \frac{[h^2]}{2} \frac{\eta + 2}{[h^{\eta+2}]} - \frac{\operatorname{sign}(\bar{q})}{k \Delta x} \left( \left[ \frac{1}{h} \right] + \frac{[h^2]}{2} \frac{\eta + 2}{[h^{\eta+2}]} \frac{[h^{\eta-1}]}{\eta - 1} \right), $$
avec la notation $[X] := X_R - X_L$. Quant à $\bar{q} = \bar{q}(q_L, q_R)$, il apparaît dans [MDBCF17] que la seule condition de consistance à assurer est

$$ (4.30) \qquad \bar{q}(q, q) = q \quad \forall q \in \mathbb{R}. $$

L'idée est d'utiliser ce degré de liberté restant sur la définition de $\bar{q}$ pour obtenir un schéma consistant à la limite de diffusion.

En procédant à la mise à l'échelle temps long et friction dominante (4.2), et en utilisant la définition (4.27) des états intermédiaires, le schéma (4.14) s'écrit finalement

$$ (4.31) \qquad \begin{aligned} h_i^{n+1} = h_i^n &+ \frac{\Delta t}{2\epsilon\Delta x} \left( \lambda_{i+\frac{1}{2}} (h_{i+1}^n - h_i^n) - \lambda_{i-\frac{1}{2}} (h_i^n - h_{i-1}^n) \right) \\ &- \frac{\Delta t}{2\epsilon\Delta x} (q_{i+1}^n - q_{i-1}^n) \\ &+ \frac{k\Delta t}{2\epsilon^3} \left( \lambda_{i+\frac{1}{2}} \frac{\bar{q}_{i+\frac{1}{2}} |\bar{q}_{i+\frac{1}{2}}| \overline{h^{-\eta}_{i+\frac{1}{2}}}}{\gamma_{i+\frac{1}{2}}} - \lambda_{i-\frac{1}{2}} \frac{\bar{q}_{i-\frac{1}{2}} |\bar{q}_{i-\frac{1}{2}}| \overline{h^{-\eta}_{i-\frac{1}{2}}}}{\gamma_{i-\frac{1}{2}}} \right), \end{aligned} $$

$$ (4.32) \qquad \begin{aligned} q_i^{n+1} = q_i^n &+ \frac{\Delta t}{2\epsilon\Delta x} \left( \lambda_{i+\frac{1}{2}} (q_{i+1}^n - q_i^n) - \lambda_{i-\frac{1}{2}} (q_i^n - q_{i-1}^n) \right) \\ &- \frac{\Delta t}{2\epsilon\Delta x} \left( \frac{(q_{i+1}^n)^2}{h_{i+1}^n} + g \frac{(h_{i+1}^n)^2}{2} - \frac{(q_{i-1}^n)^2}{h_{i-1}^n} - g \frac{(h_{i-1}^n)^2}{2} \right) \\ &- \frac{k\Delta t}{2\epsilon^3} \left( \bar{q}_{i+\frac{1}{2}} |\bar{q}_{i+\frac{1}{2}}| \overline{h^{-\eta}_{i+\frac{1}{2}}} + \bar{q}_{i-\frac{1}{2}} |\bar{q}_{i-\frac{1}{2}}| \overline{h^{-\eta}_{i-\frac{1}{2}}} \right). \end{aligned} $$

En ajoutant une condition supplémentaire sur la définition de $\bar{q}$, nous établissons que ce schéma (4.31)–(4.32) est consistant avec le problème de diffusion (4.4).

**Théorème 4.3.** *Les vitesses d'ondes* $(\lambda_{i+\frac{1}{2}})_{i \in \mathbb{Z}}$ *sont supposées définies de sorte que*

$$ (4.33) \qquad \lambda_{i+\frac{1}{2}}\gamma_{i-\frac{1}{2}} + \lambda_{i-\frac{1}{2}}\gamma_{i+\frac{1}{2}} \neq 0 \quad \forall i \in \mathbb{Z}. $$

Considérons une approximation $\bar{q}_{i+\frac{1}{2}}$ vérifiant

$$ (4.34) \qquad \bar{q}_{i+\frac{1}{2}} = 0 \text{ si et seulement si } q_i^n = q_{i+1}^n = 0. $$

Alors, quand $\epsilon \to 0$, les discrétisations données par le schéma (4.31)–(4.32) satisfont

$$ h_i^n = h_i^{n(0)} + \mathcal{O}(\epsilon) \quad \text{et} \quad q_i^n = \epsilon q_i^{n(1)} + \mathcal{O}(\epsilon^2) \quad \forall i \in \mathbb{Z}, \forall n \in \mathbb{N}, $$

$$ (4.35) \qquad h_i^{n+1(0)} = h_i^{n(0)} - \frac{\Delta t}{2\Delta x} (q_{i+1}^{n(1)} - q_{i-1}^{n(1)}) + \mathcal{O}(\Delta t \Delta x), $$

$$ (4.36) \qquad \frac{g_2(h_{i+1}^{n(0)})^2 - g_2(h_{i-1}^{n(0)})^2}{2\Delta x} = -\frac{k}{2} \left( \widetilde{q|q|}_{i+\frac{1}{2}}^{(2)} \overline{h^{-\eta}_{i+\frac{1}{2}}}^{(0)} + \widetilde{q|q|}_{i-\frac{1}{2}}^{(2)} \overline{h^{-\eta}_{i-\frac{1}{2}}}^{(0)} \right), $$

où

$$ (4.37) \qquad \overline{h^{-\eta}}_{i+\frac{1}{2}}^{(0)} = \frac{(h_{i+1}^{n(0)})^2 - (h_i^{n(0)})^2}{2} \frac{\eta+2}{(h_{i+1}^{n(0)})^{\eta+2} - (h_i^{n(0)})^{\eta+2}}, $$

$$ (4.38) \qquad \widetilde{q|q|}_{i+\frac{1}{2}}^{(2)} = \text{sign}(q_i^{n(1)} + q_{i+1}^{n(1)}) \frac{1}{2} ((q_i^{n(1)})^2 + (q_{i+1}^{n(1)})^2). $$
L'hypothèse technique (4.33) est en pratique peu restrictive et dans tous les cas tests que nous avons considérés, la définition classique (4.12) des vitesses d'ondes convenait.

Remarquons que l'approximation $\bar{q}$ doit vérifier simultanément la condition de consistance (4.30) et la restriction (4.34). En pratique, la définition suivante peut par exemple être adoptée :

$$ (4.39) \qquad \bar{q} = \mathrm{sign}(q_L + q_R) \sqrt{\frac{1}{2}(q_L^2 + q_R^2)}. $$

Il apparaît en fait que si la restriction (4.34) est nécessaire à la preuve du théorème 4.3, la consistance à la limite de diffusion est obtenue en pratique avec un choix de $\bar{q}$ vérifiant uniquement (4.30). Ce fait est illustré sur la figure 4.3, où nous comparons les résultats obtenus avec la définition (4.39) de $\bar{q}$ et ceux obtenus avec la définition adoptée dans [MDBCF17] :

$$ (4.40) \qquad \bar{q} = \begin{cases} \mathrm{sign}(q_L + q_R) \frac{2|q_L||q_R|}{|q_L| + |q_R|} & \text{si } q_L \neq 0 \text{ et } q_R \neq 0, \\ 0 & \text{sinon.} \end{cases} $$

Notons enfin qu'en pratique, le schéma limite associé à (4.31)–(4.32) ne s'écrit pas de manière explicite : des termes de viscosité apparaissant dans (4.35)–(4.36) dépendent des termes d'ordre 1 et 2 du développement asymptotique de $h$ et $q$, et ne sont pas calculables numériquement.
<figure><img src="image_24.png" /><figcaption>FIGURE 4.3 : Comparaison de la solution approchée obtenue par le schéma (4.31)–(4.32) pour différentes valeurs de $\epsilon$ avec la solution du problème limite (4.5)–(4.6).</figcaption></figure>
<figure><img src="image_25.png" /></figure>
# Troisième partie

## Analyse numérique pour des problèmes issus de la théorie cinétique
<figure><img src="image_26.png" /></figure>
## Limite de diffusion et hypocoercivité discrètes pour un modèle cinétique linéaire

Dans ce chapitre, je présente des résultats concernant l’analyse de propriétés asymptotiques d’un schéma numérique pour des équations cinétiques linéaires unidimensionnelles. L’équation considérée décrit l’évolution d’une distribution de particules $f^\varepsilon(t, x, v)$, pour un temps $t \ge 0$, une position $x \in \mathbb{T}$ et une vitesse $v \in \mathbb{R}$ :

$$ (5.1) \qquad \begin{cases} \varepsilon \frac{\partial f^\varepsilon}{\partial t} + v \frac{\partial f^\varepsilon}{\partial x} = \frac{1}{\varepsilon} Q(f^\varepsilon), \\ f^\varepsilon(0, x, v) = f_0(x, v) \ge 0. \end{cases} $$

Le paramètre $\varepsilon > 0$ est le ratio entre le libre parcours moyen des particules et la longueur caractéristique du problème, correspondant par exemple dans le contexte de la dynamique des gaz raréfiés au nombre de Knudsen adimensionné. Le système est dit en régime *cinétique* quand $\varepsilon \sim 1$, et en régime *diffusif* quand $\varepsilon \ll 1$. L’opérateur $Q$ décrit le processus de collision agissant uniquement sur la variable de vitesse $v$, et préservant la masse. Je me concentre ici sur deux opérateurs de collision particuliers :

* un opérateur de type Fokker-Planck

$$ (5.2) \qquad Q_{FP}(f)(v) = \partial_v (\partial_v f + v f), $$

* un opérateur de type BGK [BGK54] linéaire (ou relaxation)

$$ (5.3) \qquad Q_{BGK}(f)(v) = \rho M(v) - f(v), \qquad \rho = \int_{\mathbb{R}} f(v) \, dv, $$

$M$ étant la Gaussienne centrée réduite

$$ M(v) := \frac{1}{\sqrt{2\pi}} \exp \left( -\frac{|v|^2}{2} \right), $$
de moments $m_k = \int_{\mathbb{R}} |v|^k M(v) \, dv$.

Remarquons que la Gaussienne centrée réduite vérifie

$$ (5.4) \qquad M'(v) = -v M(v) \quad \forall v \in \mathbb{R}, $$

ce qui permet de réécrire l'opérateur de Fokker-Planck sous la forme

$$ (5.5) \qquad Q_{FP}(f) = \partial_v \left( M \partial_v \left( \frac{f}{M} \right) \right). $$

L'objectif des travaux présentés dans ce chapitre est de construire et d'analyser un schéma numérique efficace pour deux asymptotiques de (5.1) : la limite de diffusion $\varepsilon \to 0$ et la limite en temps long $t \to \infty$.

**Limite de diffusion.** Elle permet de capturer des dynamiques macroscopiques non triviales du modèle. Cette limite peut s'obtenir en écrivant les équations sur les premiers moments de $f^\varepsilon$ :

$$ (5.6) \qquad \begin{cases} \partial_t \rho^\varepsilon + \partial_x j^\varepsilon = 0, \\ \varepsilon^2 \partial_t j^\varepsilon + m_2 \partial_x \rho^\varepsilon + \partial_x S^\varepsilon = -j^\varepsilon, \end{cases} $$

où les moments sont définis par

$$ \rho^\varepsilon := \int_{\mathbb{R}} f^\varepsilon \, dv, \quad j^\varepsilon := \frac{1}{\varepsilon} \int_{\mathbb{R}} f^\varepsilon \, v \, dv, \quad S^\varepsilon := \int_{\mathbb{R}} (v^2 - m_2) f^\varepsilon \, dv. $$

Formellement, en prenant $\varepsilon \to 0$ dans l'équation cinétique (5.1) et les équations aux moments (5.6), nous constatons que $f^\varepsilon \to f = \rho M$ ce qui implique $S^\varepsilon \to 0$, et la seconde équation de (5.6) donne $j^\varepsilon \to -m_2 \partial_x \rho$. Ainsi, la densité limite $\rho$ est solution de l'équation de la chaleur

$$ (5.7) \qquad \frac{\partial \rho}{\partial t} - \partial_x (m_2 \partial_x \rho) = 0. $$

Cette limite de diffusion a été étudiée dans de nombreux articles (voir par exemple [BLP79, BSS84, Pou91, DGP00]).

**Hypocoercivité.** L'autre problème classique auquel je m'intéresse dans ce chapitre est le comportement en temps long de l'équation cinétique (5.1). Dans un espace fonctionnel approprié $\mathcal{X}$, le retour exponentiel à l'équilibre peut être démontré : il existe des constantes $\kappa > 0$ et $C \ge 1$ telles que

$$ (5.8) \qquad \|f^\varepsilon(t) - \mu_f M\|_{\mathcal{X}} \le C \|f_0 - \mu_f M\|_{\mathcal{X}} e^{-\kappa t}, $$

où

$$ \mu_f = \iint_{\mathbb{T} \times \mathbb{R}} f_0 \, dx \, dv. $$
Dans le cadre de la mise à l'échelle diffusive que nous considérons ici, les constantes C et $\kappa$ peuvent être choisies indépendamment de $\varepsilon$.

Depuis quelques années, des méthodes systèmatiques et robustes, dites d'hypocoercivité, ont été développées pour démontrer (5.8) (voir par exemple [HN04, MN06, Vil09, DMS15]). Ces méthodes s'inscrivent dans la classe des méthodes d'entropie. En effet, l'idée générale est de construire une fonctionnelle d'entropie modifiée $H(f^\varepsilon)$ dissipée le long des solutions. Le contrôle de l'entropie modifiée par sa dissipation permet alors d'établir la convergence, avec un taux explicite.

Dans ce chapitre, je propose une synthèse de résultats obtenus en collaboration avec Maxime Herda et Thomas Rey [11]. La section 5.1 est consacrée à la construction de schémas volumes finis implicites en temps pour l'équation (5.1) avec opérateurs de Fokker-Planck (5.2) et de BGK linéaire (5.3). Nous proposons l'étude complète de leurs limites de diffusion et de leur comportement en temps long, qui repose notamment sur des estimations uniformes en temps et en $\varepsilon$ présentées dans la section 5.2. La propriété de préservation de l'asymptotique de diffusion $\varepsilon \to 0$ de nos discrétisations est établie dans la section 5.3, et une implémentation adéquate, détaillée dans la section 5.4, permet de prendre $\varepsilon = 0$ dans le schéma. Enfin, nous prouvons dans la section 5.5 l'analogue discret de l'estimation de décroissance exponentielle (5.8), en adaptant la méthode proposée dans le cadre continu dans [DMS15].

## 5.1 Définition du schéma numérique

Considérons un domaine borné symétrique en vitesse $[-v_*, v_*]$, décomposé en $2L$ volumes de contrôle arrangés symétriquement autour de 0

$$V_j := (v_{j-\frac{1}{2}}, v_{j+\frac{1}{2}}), \quad j \in \mathcal{J} := \{-L+1, \dots, L\},$$

avec

$$v_{-L+1/2} = -v_{\star}, \quad v_{1/2} = 0, \quad v_{j+1/2} = -v_{-j+1/2}, \quad \forall j = 0, \dots, L.$$

Pour chaque cellule $V_j$, notons $\Delta v_j$ sa longueur et $v_j$ son milieu. Définissons également le maillage dual composé des cellules

$$V_{j+1/2}^* := (v_j, v_{j+1}), \quad j \in \mathcal{J}^* := \{-L, \dots, L\},$$

avec $v_{-L} := v_{-L+1/2} = -v_{\star}$ et $v_{L+1} := v_{L+1/2} = v_{\star}$.

Concernant l'espace, considérons une discrétisation du tore $\mathbb{T}$ en $N$ sous-intervalles

$$\mathcal{X}_i := (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}), \quad i \in \mathcal{I} := \mathbb{Z}/N\mathbb{Z}$$

de longueur $\Delta x_i$ et centres $x_i$.

Les volumes de contrôle dans l'espace des phases sont donc définis par

$$K_{ij} := \mathcal{X}_i \times V_j, \quad \forall (i, j) \in \mathcal{I} \times \mathcal{J}.$$
Finalement, notons $\Delta t > 0$ le pas de temps et $t^n = n\Delta t$ pour tout $n \ge 0$.

Un point crucial pour établir nos résultats est de considérer une discrétisation adéquate de la Maxwellienne $M(v)$. Dans le cas de l'opérateur BGK (5.3), considérons des valeurs par cellules $M = (M_j)_{j \in \mathcal{J}}$ vérifiant les hypothèses suivantes :

$$ (5.9) \qquad \left\{ \begin{array}{l} M_j > 0, \quad M_j = M_{-j+1}, \quad \forall j = 1, \dots, L; \\ \displaystyle\sum_{j \in \mathcal{J}} M_j \Delta v_j = 1; \\ 0 < \underline{m}_2 \le m_2^{\Delta v} \le \bar{m}_2, \quad m_4^{\Delta v} \le \bar{m}_4, \end{array} \right. $$

où pour $k \in \mathbb{N}$,

$$ m_k^{\Delta v} := \sum_{j \in \mathcal{J}} |v_j|^k M_j \Delta v_j $$

et $\underline{m}_2, \bar{m}_2, \bar{m}_4$ sont des constantes.

Dans le cas de l'opérateur de Fokker-Planck (5.2), nous considérons des valeurs aux interfaces $(M_{j+1/2}^*)_{j \in \mathcal{J}^*} \in \mathbb{R}^{\mathcal{J}^*}$ vérifiant

$$ (5.10) \qquad \left\{ \begin{array}{l} M_{j+1/2}^* = M_{-j+1/2}^*, \quad \forall j \in \mathcal{J}^*; \\ M_{L+1/2}^* = M_{-L+1/2}^* = 0; \\ M_j := \frac{M_{j-1/2}^* - M_{j+1/2}^*}{v_j \Delta v_j} > 0, \quad \forall j \in \mathcal{J}; \\ \displaystyle\sum_{j \in \mathcal{J}} M_j \Delta v_j = 1; \\ 0 < \underline{m}_2 \le m_2^{\Delta v} \le \bar{m}_2, \quad m_4^{\Delta v} \le \bar{m}_4. \end{array} \right. $$

De manière pratique, de telles Maxwelliennes discrètes peuvent facilement être construites. La deuxième hypothèse dans (5.10) dans le cas Fokker-Planck est nécessaire pour des raisons techniques lors des intégrations par parties discrètes. Pour des domaines suffisamment grands, cette hypothèse est pertinente du fait de la décroissance exponentielle de la Gaussienne. La troisième hypothèse dans (5.10) est un analogue discret de la relation (5.4), permettant de réécrire au niveau continu l'opérateur de Fokker-Planck sous la forme (5.5), ce qui nous est essentiel au niveau discret.

Remarquons enfin que les hypothèses de symétrie des vitesses et des Maxwelliennes impliquent

$$ (5.11) \qquad \sum_{j \in \mathcal{J}} \Delta v_j v_j M_j = 0. $$

Présentons maintenant les discrétisations de l'équation (5.1) dans les cas Fokker-Planck et BGK. Nous considérons des schémas volumes finis, avec une discrétisation temporelle de type Euler implicite. Afin d'alléger les notations, j'omets à partir de maintenant les exposants $\epsilon$ dans la discrétisation de $f^\epsilon$ et ses moments, même si ces quantités dépendent bien entendu de $\epsilon$.
**Schéma pour l'équation de Fokker-Planck.** Pour tous $i \in \mathcal{I}$, $j \in \mathcal{J}$, $n \ge 0$,

$$ (5.12) \qquad \varepsilon \Delta x_i \Delta v_j (f_{ij}^{n+1} - f_{ij}^n) + \Delta t \left( \mathcal{F}_{i+\frac{1}{2},j}^{n+1} - \mathcal{F}_{i-\frac{1}{2},j}^{n+1} \right) = \frac{\Delta t}{\varepsilon} \left( \mathcal{G}_{i,j+\frac{1}{2}}^{n+1} - \mathcal{G}_{i,j-\frac{1}{2}}^{n+1} \right). $$

Le transport libre est approché par un flux centré :

$$ (5.13) \qquad \mathcal{F}_{i+\frac{1}{2},j}^{n+1} = v_j \frac{f_{i+1,j}^{n+1} + f_{ij}^{n+1}}{2} \Delta v_j, \quad \forall j \in \mathcal{J}, \quad \forall i \in \mathcal{I}, $$

tandis que la définition des flux en vitesse $\mathcal{G}_{i,j+\frac{1}{2}}$ est basée sur la forme (5.5) de l'opérateur de Fokker-Planck :

$$ (5.14) \qquad \mathcal{G}_{i,j+\frac{1}{2}}^{n+1} = \mathcal{M}_{j+\frac{1}{2}}^* \left( \frac{f_{i,j+1}^{n+1}}{\mathcal{M}_{j+1}} - \frac{f_{ij}^{n+1}}{\mathcal{M}_j} \right) \frac{1}{\Delta v_{j+\frac{1}{2}}} \Delta x_i, \quad \forall j \in \mathcal{J}^* \setminus \{-L, L\}, \quad \forall i \in \mathcal{I}. $$

Nous imposons finalement des conditions au bord de flux nul :

$$ (5.15) \qquad \mathcal{G}_{i,-L+\frac{1}{2}}^{n+1} = \mathcal{G}_{i,L+\frac{1}{2}}^{n+1} = 0, \quad \forall i \in \mathcal{I}. $$

**Schéma pour l'équation de BGK.** Pour tous $i \in \mathcal{I}$, $j \in \mathcal{J}$, $n \ge 0$,

$$ (5.16) \qquad \varepsilon \Delta x_i \Delta v_j (f_{ij}^{n+1} - f_{ij}^n) + \Delta t \left( \mathcal{F}_{i+\frac{1}{2},j}^{n+1} - \mathcal{F}_{i-\frac{1}{2},j}^{n+1} \right) = \frac{\Delta t}{\varepsilon} \Delta x_i \Delta v_j \left( \rho_i^{n+1} \mathcal{M}_j - f_{ij}^{n+1} \right), $$

où le flux $\mathcal{F}_{i+\frac{1}{2},j}$ est défini par (5.13) et la densité discrète est donnée par :

$$ (5.17) \qquad \rho_i^n := \sum_{j \in \mathcal{J}} \Delta v_j f_{ij}^n. $$

## 5.2 Estimations uniformes discrètes

Dans cette section, je présente les estimations uniformes discrètes qui sont le point clé pour étudier à la fois la limite de diffusion du schéma et le retour exponentiel à l'équilibre des solutions discrètes.

Le premier résultat est une version discrète de l'estimation d'entropie $L^2$ à poids $\gamma = 1/M$. La norme $L^2$ discrète à poids $\gamma$ étant définie pour $g = (g_{ij})_{i \in \mathcal{I}, j \in \mathcal{J}}$ par :

$$ \|g\|_{2,\gamma}^2 := \sum_{(i,j) \in \mathcal{I} \times \mathcal{J}} |g_{ij}|^2 \gamma_j \Delta x_i \Delta v_j, $$

le lemme suivant est obtenu.

**Lemme 5.1.** Soit $\mathcal{M}$ une Maxwellienne discrète vérifiant (5.10) (resp. (5.9)) et $(f_{ij}^n)_{i \in \mathcal{I}, j \in \mathcal{J}, n \in \mathbb{N}}$ solution du schéma (5.12) (resp. (5.16)). Alors, pour tous $n \ge 0$,

$$ (5.18) \qquad \frac{\|f^{n+1}\|_{2,\gamma}^2 - \|f^n\|_{2,\gamma}^2}{2\Delta t} + \frac{1}{\varepsilon^2} \|f^{n+1} - \rho^{n+1}\mathcal{M}\|_{2,\gamma}^2 \le 0. $$
En particulier, les estimations suivantes sont vérifiées :

$$
(5.19) \quad \max \left( \sup_{n \ge 0} \|f^n\|_{2,\gamma}^2, \frac{2}{\epsilon^2} \sum_{n=1}^{\infty} \Delta t \|f^n - \rho^n \mathcal{M}\|_{2,\gamma}^2 \right) \le \|f^0\|_{2,\gamma}^2.
$$

Ce résultat est démontré en utilisant la définition de la densité discrète $(\rho_i)_i$ et les hypothèses (5.9) sur la Maxwellienne discrète dans le cas BGK. Dans le cas Fokker-Planck, la discrétisation particulière (5.14) des flux ainsi qu'une inégalité de Poincaré à poids gaussien discrète sont utilisées.

Puisque les schémas (5.12) et (5.16) sont des systèmes linéaires de dimension finie, l'estimation (5.19) permettant en particulier de montrer l'unicité de la solution discrète, elle en assure également l'existence.

Le deuxième lemme concerne des estimations uniformes sur les moments discrets définis par :

$$
(5.20) \quad J_i^n := \frac{1}{\epsilon} \sum_{j \in \mathcal{J}} \Delta v_j v_j f_{ij}^n, \quad S_i^n := \sum_{j \in \mathcal{J}} \Delta v_j (v_j^2 - m_2^{\Delta v}) f_{ij}^n,
$$

Ces estimations sont obtenues grâce à des inégalités de Cauchy-Schwarz et les majorations
(5.19) obtenues dans le lemme précédent.

**Lemme 5.2.** *Sous les hypothèses du lemme 5.1, les moments discrets vérifient les estimations suivantes :*

$$
(5.21) \qquad \|\rho^n\|_2 \le \|f^n\|_{2,\gamma},
$$

$$
(5.22) \qquad \varepsilon \|J^n\|_2 \le (m_2^{\Delta v})^{1/2} \|f^n\|_{2,\gamma},
$$

$$
(5.23) \qquad \varepsilon \|J^n\|_2 \le (m_2^{\Delta v})^{1/2} \|f^n - \rho^n \mathcal{M}\|_{2,\gamma},
$$

$$
(5.24) \quad \|S^n\|_2 \le (m_4^{\Delta v} - (m_2^{\Delta v})^2)^{1/2} \|f^n - \rho^n \mathcal{M}\|_{2,\gamma}.
$$

En particulier, il existe une constante $C > 0$ dépendant uniquement des bornes uniformes sur $m_2^{\Delta v}$ et $m_4^{\Delta v}$ telle que

$$
(5.25) \quad \max \left( \sup_{n \in \mathbb{N}} \|\rho^n\|_2, \varepsilon \sup_{n \in \mathbb{N}} \|J^n\|_2, \left( \sum_{n=0}^{\infty} \Delta t \|J^n\|_2^2 \right)^{1/2}, \frac{1}{\varepsilon} \left( \sum_{n=0}^{\infty} \Delta t \|S^n\|_2^2 \right)^{1/2} \right) \\
\le C \|f^0\|_{2,\gamma}.
$$
## 5.3 Préservation de l'asymptotique de diffusion

Comme dans le cadre continu, l'étude de la limite de diffusion se base sur une version discrète des équations aux moments (5.6), donnée dans le lemme suivant.

**Lemme 5.3.** Soit $\mathcal{M}$ une Maxwellienne discrète vérifiant (5.10) (resp. (5.9)) et $(f_{ij}^n)_{i \in \mathcal{I}, j \in \mathcal{J}, n \in \mathbb{N}}$ solution du schéma (5.12) (resp. (5.16)). Alors les moments discrets vérifient : pour tous $i \in \mathcal{I}$, $n \ge 0$,

$$ (5.26) \qquad \Delta x_i (\rho_i^{n+1} - \rho_i^n) + \Delta t (J_{i+\frac{1}{2}}^{n+1} - J_{i-\frac{1}{2}}^{n+1}) = 0, $$

$$ (5.27) \qquad \varepsilon^2 \Delta x_i (J_i^{n+1} - J_i^n) + \Delta t (S_{i+\frac{1}{2}}^{n+1} - S_{i-\frac{1}{2}}^{n+1}) + \Delta t m_2^{\Delta v} (\rho_{i+\frac{1}{2}}^{n+1} - \rho_{i-\frac{1}{2}}^{n+1}) = - \Delta t \Delta x_i J_i^{n+1}, $$

où les valeurs approchées aux interfaces sont définies par la moyenne

$$ Q_{i+\frac{1}{2}} = \frac{Q_i + Q_{i+1}}{2}, \quad Q = \rho, J, S. $$

À partir de là, la propriété de préservation de l'asymptotique de diffusion $\varepsilon \to 0$ des schémas peut être rigoureusement démontrée.

**Théorème 5.1.** Soit $\mathcal{M}$ une Maxwellienne discrète vérifiant (5.10) (resp. (5.9)), soit $f_{\varepsilon}^n = (f_{ij}^n)_{i \in \mathcal{I}, j \in \mathcal{J}}$ la solution de (5.12) (resp. (5.16)). Alors il existe $\rho^n = (\rho_i^n)_{i \in \mathcal{I}}$ tel que quand $\varepsilon \to 0$,

$$ f_{\varepsilon}^n \to \rho^n \mathcal{M} \quad \text{dans } \mathbb{R}^{\mathcal{I} \times \mathcal{J}}, \quad \text{pour tous } n \ge 1, $$

et la densité limite $\rho^n$ vérifie le schéma aux différences finies suivant pour l'équation de la chaleur

$$ (5.28) \qquad \Delta x_i \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} = \frac{m_2^{\Delta v}}{2} \left( (D_x \rho^{n+1})_{i+1} - (D_x \rho^{n+1})_{i-1} \right), \quad \forall i \in \mathcal{I}, $$

le gradient discret $D_x \rho$ étant défini sur le maillage primal par

$$ (5.29) \qquad (D_x \rho)_i = \frac{\rho_{i+1/2} - \rho_{i-1/2}}{\Delta x_i}. $$

En utilisant les résultats établis dans les lemmes précédents, la preuve de ce théorème est assez directe. En effet, l'estimation (5.25) sur la densité discrète donne l'existence d'une sous-suite convergente $\rho_{\varepsilon_k}^n \to \rho^n$. Grâce à l'estimation (5.19), la suite $(f_{\varepsilon_k}^n)_k$ converge vers $(\rho_i^n M_j)_{i \in \mathcal{I}, j \in \mathcal{J}}$. Le schéma limite (5.28) s'obtient en combinant les équations aux moments (5.26) et (5.27), et en passant à la limite $\varepsilon \to 0$ en utilisant les estimations (5.25). Finalement, le schéma limite ayant une unique solution, la convergence de toute la suite est obtenue.

Cette propriété de préservation de l'asymptotique de diffusion au niveau discret est illustrée dans le cas Fokker-Planck sur la figure 5.1.
<figure><img src="image_27.png" /><figcaption>FIGURE 5.1 : Comparaison de la solution de l'équation de la chaleur obtenue par le schéma (5.28) avec la solution obtenue par le schéma (5.34) pour différents ε, aux temps t = 0.05, 0.1, 0.15 et 10.</figcaption></figure>
## 5.4 Implémentation

En pratique, l'implémentation des schémas écrits sous la forme (5.12) (resp. (5.16)) présente deux problèmes. D'une part, la matrice à inverser est mal conditionnée pour $\varepsilon$ petit. Cela conduit à des erreurs numériques importantes proche du régime diffusif, et détériore donc la propriété de préservation de l'asymptotique. D'autre part, l'évaluation de la différence de quantités qui sont du même ordre, telle que par exemple $\|f^n - \mu_f M\|_{2,M}$ quand $n$ est grand, conduit à une accumulation d'erreurs machine et donc à de mauvaises approximations de ces différences.

Pour résoudre ces problèmes, nous proposons une écriture équivalente des schémas, fondée sur une formulation micro-macro. Ainsi, introduisons les inconnues micro et macro, $(h_{ij}^n)_{i \in \mathcal{I}, j \in \mathcal{J}, n \in \mathbb{N}}$ et $(\lambda_i^n)_{i \in \mathcal{I}, n \in \mathbb{N}}$ respectivement, vérifiant

$$ (5.30) \qquad f_{ij}^n = \mu_f M_j + \lambda_i^n M_j + \varepsilon h_{ij}^n M_j, $$

$$ (5.31) \qquad \lambda_i^n = \rho_i^n - \mu_f = \sum_{j \in \mathcal{I}} f_{ij}^n \Delta v_j - \mu_f. $$

En conséquence immédiate de ces définitions, la relation suivante est vérifiée

$$ (5.32) \qquad \sum_{j \in \mathcal{J}} h_{ij}^n M_j \Delta v_j = 0 \quad \forall i \in \mathcal{I}, \forall n \in \mathbb{N}. $$

Avec ces nouvelles inconnues, le schéma (5.26) donne l'équation d'évolution macroscopique suivante :

$$ (5.33) \qquad \Delta x_i (\lambda_i^{n+1} - \lambda_i^n) + \Delta t \sum_{j \in \mathcal{J}} \Delta v_j v_j M_j \frac{h_{i+1,j}^{n+1} - h_{i-1,j}^{n+1}}{2} = 0, \quad \forall i \in \mathcal{I}. $$

**Cas Fokker-Planck.** En utilisant (5.30) dans le schéma cinétique initial (5.12) et en soustrayant (5.33) multiplié par $-\Delta v_j \varepsilon$, nous obtenons l'équation d'évolution microscopique suivante :

$$ \begin{align} (5.34) \qquad & \varepsilon^2 \Delta x_i \Delta v_j (h_{ij}^{n+1} - h_{ij}^n) + \Delta t \Delta v_j v_j \frac{\lambda_{i+1}^{n+1} - \lambda_{i-1}^{n+1}}{2} \nonumber \\ & + \varepsilon \Delta t \Delta v_j \left( v_j \frac{h_{i+1,j}^{n+1} - h_{i-1,j}^{n+1}}{2} - \sum_{k \in \mathcal{J}} \Delta v_k v_k M_k \frac{h_{i+1,k}^{n+1} - h_{i-1,k}^{n+1}}{2} \right) \nonumber \\ & = \Delta t \Delta x_i \left( \frac{M_{j+1/2}^*}{\Delta v_{j+1/2} M_j} (h_{i,j+1}^{n+1} - h_{ij}^{n+1}) - \frac{M_{j-1/2}^*}{\Delta v_{j-1/2} M_j} (h_{ij}^{n+1} - h_{i,j-1}^{n+1}) \right). \end{align} $$
**Cas BGK.** En utilisant (5.30) dans le schéma cinétique initial (5.16), nous obtenons l'équation d'évolution microscopique suivante :

$$ (5.35) \quad \begin{aligned} & \varepsilon^2 \Delta x_i \Delta v_j (h_{ij}^{n+1} - h_{ij}^n) + \Delta t \Delta v_j v_j \frac{\lambda_{i+1}^{n+1} - \lambda_{i-1}^{n+1}}{2} \\ & \qquad + \varepsilon \Delta t \Delta v_j \left( v_j \frac{h_{i+1,j}^{n+1} - h_{i-1,j}^{n+1}}{2} - \sum_{k \in \mathcal{J}} \Delta v_k v_k \mathcal{M}_k \frac{h_{i+1,k}^{n+1} - h_{i-1,k}^{n+1}}{2} \right) \\ & = -\Delta t \Delta x_i \Delta v_j h_{ij}^{n+1}. \end{aligned} $$

Partant de ces calculs, trois schémas peuvent être définis :

(S<sub>f</sub><sup>ε</sup>) L'inconnue (f<sub>ij</sub><sup>n</sup>)<sub>i,j,n</sub> est solution de (5.12) dans le cas Fokker-Planck et (5.16) dans le cas BGK. Les inconnues micro et macro (h<sub>ij</sub><sup>n</sup>)<sub>i,j,n</sub> et (λ<sub>i</sub><sup>n</sup>)<sub>i,n</sub> sont définies par (5.30) and (5.31).

(S<sub>λ,h</sub><sup>ε</sup>) Les inconnues micro et macro (h<sub>ij</sub><sup>n</sup>)<sub>i,j,n</sub> et (λ<sub>i</sub><sup>n</sup>)<sub>i,n</sub> sont solutions de (5.33)–(5.34) dans le cas Fokker-Planck et (5.33)–(5.35) dans le cas BGK. L'inconnue (f<sub>ij</sub><sup>n</sup>)<sub>i,j,n</sub> est définie par (5.30).

(S̃<sub>λ,h</sub><sup>ε</sup>) Les inconnues micro et macro (h<sub>ij</sub><sup>n</sup>)<sub>i,j,n</sub> et (λ<sub>i</sub><sup>n</sup>)<sub>i,n</sub> sont solutions du système surdéterminé (5.32)–(5.33)–(5.34) dans le cas Fokker-Planck et (5.32)–(5.33)–(5.35) dans le cas BGK. L'inconnue (f<sub>ij</sub><sup>n</sup>)<sub>i,j,n</sub> est définie par (5.30).

Dans le régime cinétique ε > 0, ces trois schémas sont bien posés et équivalents, mais ce n'est plus le cas dans le régime diffusif ε = 0.

**Proposition 5.1.** Les trois versions du schéma vérifient les propriétés suivantes.

(i) Si ε > 0, (S<sub>f</sub><sup>ε</sup>), (S<sub>λ,h</sub><sup>ε</sup>) et (S̃<sub>λ,h</sub><sup>ε</sup>) définissent les mêmes uniques solutions (f<sub>ij</sub><sup>n</sup>)<sub>i,j,n</sub>, (h<sub>ij</sub><sup>n</sup>)<sub>i,j,n</sub> et (λ<sub>i</sub><sup>n</sup>)<sub>i,n</sub>.

(ii) Si ε = 0, l'équivalence des schémas n'est plus vérifiée. Plus précisément, (S<sub>f</sub><sup>0</sup>) est mal posé, (S<sub>λ,h</sub><sup>0</sup>) est bien posé dans le cas BGK et mal posé dans le cas Fokker-Planck, et (S̃<sub>λ,h</sub><sup>0</sup>) est bien posé.

Ce résultat repose sur une étude du noyau de l'opérateur de collision discret associé à chacun des schémas.

Nous déduisons de cette proposition que le conditionnement du système linéaire apparaissant dans (S̃<sub>λ,h</sub><sup>ε</sup>) est uniformément borné par rapport à ε ∈ [0, 1].

## 5.5 Hypocoercivité discrète

Je présente maintenant l'adaptation au cadre discret de la méthode d'hypocoercivité proposée dans [DMS15]. Notre résultat principal, analogue discret de (5.8), est présenté dans le théorème suivant.
**Théorème 5.2.** Soit $\mathcal{M}$ une Maxwellienne discrète vérifiant (5.10) (resp. (5.9)), et sup-
posons le nombre de points $N$ de la discrétisation spatiale impair. Alors il existe des
constantes $C \ge 1$ et $\beta > 0$ telles que pour tous $\varepsilon \in (0,1)$, tous $\Delta t \le \Delta t_{max}$, la solu-
tion $(f^n)_{i \in \mathcal{I}, j \in \mathcal{J}, n \in \mathbb{N}}$ de (5.12) (resp. (5.16)) vérifie

$$
(5.36) \qquad \|f^n - \mu_f \mathcal{M}\|_{2,\gamma} \le C \|f^0 - \mu_f \mathcal{M}\|_{2,\gamma} e^{-\frac{\beta}{2} t^n}.
$$

De plus, les constantes C et β ne dépendent pas des paramètres de discrétisation, et Δt_max peut être choisi de façon arbitraire.

Sans perte de généralité, la masse totale discrète étant préservée par les schémas, il
suffit de démontrer le résultat dans le cas $\mu_f = 0$. La preuve repose sur l'étude d'une
fonctionnelle d'entropie modifiée définie par
(5.37)

$$
\mathbb{H}_{\eta}(f^n) := \frac{1}{2} \|f^n\|_{2,\gamma}^2 + \eta \epsilon^2 \sum_{i \in \mathcal{I}} \Delta x_i J_i^n (D_{x\phi})_i^n + \frac{\eta \epsilon^2}{2} \sum_{i \in \mathcal{I}} \Delta x_i \frac{\left( (D_{x\phi})_i^n - (D_{x\phi})_{i-1}^{n-1} \right)^2}{\Delta t},
$$

$\eta > 0$ étant un paramètre libre à déterminer, et $(\phi_i^n)_{i \in \mathcal{I}}$ vérifiant le problème de Poisson discret suivant

$$
(5.38) \qquad - \frac{(D_x \phi)^n_{i+1} - (D_x \phi)^n_{i-1}}{2} = \Delta x_i \rho_i^n, \quad \forall i \in \mathcal{I},
$$

$$
(5.39) \qquad \sum_{i \in \mathcal{I}} \Delta x_i \phi_i^n = 0.
$$

Les deux premiers termes de $\mathbb{H}_\eta(f^n)$ correspondent à une discrétisation de la fonction-
nelle considérée dans le cadre continu [DMS15], et le troisième terme est nécessaire pour
contrôler des termes de restes induits par la discrétisation temporelle, mais d'ordre $\Delta t$.
Si les paramètres de discrétisation tendent vers 0, l'entropie modifiée continue est donc
retrouvée.

Le problème de Poisson (5.38)–(5.39) est bien posé si et seulement si $N$ est pair. Dans
ce cas-là, des estimations sur $(\phi_i^n)_{i \in \mathcal{I}}$ par les moments discrets de $f$ sont obtenues :

$$
(5.40) \qquad \|D_x \phi^n\|_2 \le C_P \|\rho^n\|_2,
$$

$$
(5.41) \qquad \left\| \frac{(D_x \phi)_i^{n+1} - (D_x \phi)_i^n}{\Delta t} \right\|_2 \le \| J_i^{n+1} \|_2,
$$

$C_P$ étant la constante de l'inégalité de Poincaré discrète sur $\mathbb{T}$, qui peut être prouvée dans
le cas $N$ impair. Cette hypothèse d'imparité résulte de la discrétisation particulière du
gradient (5.29) intervenant ici.

La preuve s'articule classiquement en deux grandes étapes. D'une part, pour $\eta > 0$
suffisamment petit, la racine de l'entropie modifiée $\mathbb{H}_\eta(f^n)$ est équivalente à la norme
$\| \cdot \|_{2,\gamma}$. Ce résultat, énoncé dans le lemme suivant, s'obtient en utilisant les estimations
(5.39) et (5.40) sur $(\phi_i^n)_{i \in \mathcal{I}}$, ainsi que le contrôle des moments discrets (lemme 5.2).
**Lemme 5.4.** Sous les hypothèses du théorème 5.2, en supposant $\varepsilon \le 1$ et $\Delta t \le \Delta t_{max}$, l'estimation suivante est vérifiée pour tous $n \ge 1$
(5.42)

$$ \left( \frac{1}{2} - \eta (m_2^{\Delta v})^{1/2} \right) \|f^n\|_{2,\gamma}^2 \le \mathbb{H}_{\eta}(f^n) \le \left( \frac{1}{2} + \eta (m_2^{\Delta v})^{1/2} C_P + \frac{\eta}{2} m_2^{\Delta v} \Delta t_{max} \right) \|f^n\|_{2,\gamma}^2. $$

D'autre part, une estimation d'entropie-dissipation discrète est démontrée, basée elle aussi sur l'utilisation de l'estimation des moments discrets (lemme 5.2) et les estimations (5.39) et (5.40) sur $(\phi_i^n)_{i \in \mathcal{I}}$, ainsi que sur l'estimation d'entropie $L^2$ à poids $\gamma$ (lemme 5.1). Le même type de calculs que dans le cadre continu sont mis en œuvre pour prouver ce résultats, avec des termes additionnels à contrôler provenant de la discrétisation temporelle.

**Proposition 5.2.** Sous les hypothèses du théorème 5.2, il existe $\eta_2 > 0$ tel que pour tous $\varepsilon \le 1$, $\Delta t \le \Delta t_{max}$ et $\eta \le \eta_2$,

$$ (5.43) \quad \frac{\mathbb{H}_{\eta}(f^{n+1}) - \mathbb{H}_{\eta}(f^n)}{\Delta t} + K(\eta) \left( \|f^{n+1} - \rho^{n+1} \mathcal{M}\|_{2,\gamma}^2 + \|\rho^{n+1}\|_{2}^2 \right) \le 0, \quad \forall n \ge 1, $$
avec $K(\eta) = \frac{1}{2} \min(1 - \eta \bar{m}_2, \eta \underline{m}_2)$.

Le résultat d'hypocoercivité discrète (5.36) s'obtient finalement en combinant (5.42) et (5.43), et en choisissant $\eta$ de manière adéquate.

La figure 5.2 fournit une illustration du taux de décroissance exponentielle observé dans le cas Fokker-Planck, pour différentes valeurs de $\varepsilon$, et met en évidence une borne inférieure uniforme en $\varepsilon$ de ce taux. Le résultat du théorème 5.2 pourrait être affiné en suivant précisément la dépendance en $\varepsilon$ du taux de convergence $\beta$.

Enfin, remarquons que dans le cas BGK, notre résultat s'applique plus largement que dans le cas d'un équilibre gaussien, et nécessite uniquement que les hypothèses (5.9) soient vérifiées. À titre d'exemple, le cas d'un équilibre non gaussien à queue polynomiale

$$ (5.44) \qquad M(v) = \frac{1}{1 + 0.1|v|^6} (\cos(\pi v) + 1.1) $$

est présenté dans la figure 5.3.

## 5.6 Perspectives

**Extension avec un terme de potentiel externe confinant.** Une première extension naturelle consiste à ajouter un terme $\partial_x V \partial_v f^\varepsilon$ dans l'équation (5.1), modélisant le comportement de particules soumises à une force dérivant d'un champ $V$. Le champ $V$ étant supposé connu, le problème reste dans ce cas linéaire. Pour un potentiel confinant $V$ vérifiant des conditions de régularité et de croissance à l'infini, la convergence exponentielle vers l'équilibre [HN04, Vil09, DMS15] ainsi que la limite de diffusion [Pou91] sont étudiées dans le cadre continu.

Partant du schéma que nous avons proposé pour (5.1), l'objectif est d'ajouter une discrétisation appropriée du terme $\partial_x V \partial_v f^\varepsilon$, de sorte à obtenir les estimations nécessaires pour :
<figure><img src="image_28.png" /><figcaption>FIGURE 5.2 : Comparaison des taux de convergence de $\|f^{\varepsilon} - \mu_f M\|_{2,\gamma}$ pour différentes valeurs de $\varepsilon$, dans le cas Fokker-Planck.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>FIGURE 5.3 : Cas BGK non gaussien : distribution à l'équilibre (gauche), et évolution de la norme $L^2$ à poids de la différence entre $f^{\varepsilon}$ et l'équilibre global.</figcaption></figure>
* adapter au cadre discret le résultat d'hypocoercivité $L^2$ démontré dans [DMS15],

* démontrer rigoureusement la préservation de la limite de diffusion par le schéma construit.

Notons que dans un article récent [BF], la préservation de la limite de diffusion et l'hypocoercivité au niveau discret sont établies pour l'équation de Vlasov-Fokker-Planck approchée en vitesse par une méthode spectrale de Hermite et en espace par une méthode volumes finis classique. Mon objectif est d'obtenir ce type de résultats pour une discrétisation de type volumes finis à la fois en espace et en vitesse, ce qui fait apparaître d'autres difficultés.

**Extension au couplage avec l'équation de Poisson.** Une fois ces premiers résultats établis, l'idée est de passer au système de Vlasov-Fokker-Planck-Poisson

$$ (5.45) \qquad \varepsilon \partial_t f^\varepsilon + v \partial_x f^\varepsilon - (\partial_x V + \partial_x \Psi^\varepsilon) \partial_v f^\varepsilon = \frac{1}{\varepsilon} \partial_v (v f^\varepsilon + \partial_v f^\varepsilon), $$

$$ (5.46) \qquad - \partial_{xx} \Psi^\varepsilon = \rho^\varepsilon := \int_{\mathbb{R}} f^\varepsilon dv. $$

Ici, le potentiel auto-consistant $\Psi^\varepsilon$, défini par le couplage avec l'équation de Poisson (5.46), correspond à des forces électrostatiques répulsives et le potentiel extérieur $V$ est toujours supposé confinant. Du fait du couplage avec l'équation de Poisson, le problème n'est plus linéaire, ce qui complique fortement l'étude par rapport à celle menée dans les cas linéaires. Sur le modèle de ce que nous avons fait dans le cas très simplifié (5.1), l'idée serait d'étudier au niveau discret à la fois la limite de diffusion $\varepsilon \to 0$ et la convergence vers l'équilibre quand $t \to \infty$, en adaptant des outils et méthodes du cadre continu. La convergence vers le système de dérive-diffusion-Poisson quand $\varepsilon \to 0$ est par exemple étudiée dans le cas $V = 0$ dans [PS00, Gou05]. Concernant le comportement en temps long, la piste qui semble la plus prometteuse serait d'adapter au cadre discret les travaux récents proposés dans [ADLT21]. Dans cet article, les auteurs considèrent une version linéarisée autour de l'équilibre du système (5.45)-(5.46) et établissent rigoureusement la convergence exponentielle vers l'équilibre de la solution du linéarisé quand $t \to \infty$. À partir de cette étude du linéarisé, il est possible en dimension 1 d'en déduire le résultat pour le problème non linéaire de départ.
Analyse de méthodes spectrales de Hermite pour le système de
Vlasov-Poisson

Ce chapitre est consacré à l'étude de discrétisations du système de Vlasov-Poisson écrit sous la forme d'un système hyperbolique en utilisant des polynômes de Hermite en vitesse. Contrairement au modèle étudié dans le chapitre 5, le modèle de Vlasov-Poisson considéré ici ne prend pas en compte de collisions entre les particules. La seule interaction décrite est la force du champ moyen créée par les effets électrostatiques. Le système de Vlasov-Poisson décrit donc le comportement d'une distribution de particules $f = f(t, x, v)$ et du champ électrique auto-consistant $E = E(t, x)$. Je le considère dans ce chapitre sous sa forme adimensionnée, périodique en espace :

$$
(6.1) \qquad \left\{
\begin{array}{l}
\displaystyle \frac{\partial f}{\partial t} + v \frac{\partial f}{\partial x} + E \frac{\partial f}{\partial v} = 0, \\
\\
\displaystyle \frac{\partial E}{\partial x} = \rho - \rho_0,
\end{array}
\right.
$$

où $\rho_0$ est une constante assurant la quasi-neutralité du plasma :

$$
(6.2) \quad \int_{\mathbb{T}} \rho(t,x)dx = \int_{\mathbb{T}} \int_{\mathbb{R}} f(t,x,v)dvdx = \rho_0, \quad \forall t \ge 0.
$$

Un grand nombre de méthodes numériques ont été proposées pour approcher les so-
lutions de ce problème. Ici, l'idée est d'écrire le système (6.1) sous la forme d'un système
hyperbolique en utilisant des polynômes de Hermite en la variable vitesse. Des méthodes de
Galerkin avec un petit nombre de fonctions de base sont alors utilisées. Deux approches de
discrétisations ont été formalisées dans [Hol96] : la première, dite symétriquement pondé-
rée, utilise les fonctions de Hermite standards à la fois comme fonctions de base et comme
fonctions tests, alors que la deuxième, dite asymétriquement pondérée, utilise un ensemble
de fonctions tests différent de l'ensemble des fonctions de base. L'approche symétrique ne
permet pas de conserver simultanément la masse et la quantité de mouvement, mais sa stabilité $L^2$ peut être établie. Au contraire, l'approche asymétrique assure la conservation simultanée de la masse, quantité de mouvement et énergie totale, mais sa stabilité $L^2$ n'est pas obtenue.

Dans les travaux [10, 9], obtenus en collaboration avec Francis Filbet, nous étudions la stabilité et la convergence d'une méthode de Hermite asymétrique. L'idée principale pour obtenir la stabilité est d'introduire une fonction de mise à l'échelle $t \mapsto \alpha(t)$ bien adaptée aux variations de la fonction de distribution au cours du temps. Cette idée a déjà été employée pour stabiliser des méthodes spectrales de Hermite dans le contexte d'équations de diffusion linéaires [MST05] et d'équations de convection-diffusion non linéaires [MZ07]. Plus précisément, soit $\mu_t$ la mesure définie par

$$ (6.3) \qquad d\mu_t = \omega(t, v)dxdv $$

avec

$$ (6.4) \qquad \omega(t, v) := \sqrt{2\pi} \exp \left( \frac{\alpha^2(t) |v|^2}{2} \right), $$

et introduisons l'espace $L^2$ à poids suivant :

$$ L^2(d\mu_t) := \left\{ g : \mathbb{T} \times \mathbb{R} \to \mathbb{R} : \iint_{\mathbb{T} \times \mathbb{R}} |g(x, v)|^2 d\mu_t < +\infty \right\}. $$

Remarquons que pour toutes fonctions $\alpha, \tilde{\alpha} : \mathbb{R}_+ \to \mathbb{R}_+$ telles que $\tilde{\alpha}(t) \le \alpha(t)$ pour tout $t \ge 0$, les poids correspondants $\omega, \tilde{\omega}$ définis par (6.4) satisfont $\tilde{\omega}(t, v) \le \omega(t, v)$ pour tout $t, v$, et donc la norme $L^2$ à poids $\omega$ contrôle la norme $L^2$ à poids $\tilde{\omega}$. En particulier pour $\tilde{\alpha} = 0$, la norme $L^2(dx \, dv)$ classique est majorée par la norme $L^2(d\mu_t)$.

L'objectif est de déterminer une fonction positive $\alpha$ permettant d'établir un résultat de stabilité dans la norme associée à cet espace. Si $(f, E)$ est solution du système de Vlasov-Poisson, quelques manipulations impliquant (6.1) ainsi que l'application d'une inégalité de Young conduisent à l'estimation suivante, pour tous $\eta > 0$ :

$$ (6.5) \qquad \frac{1}{2} \frac{d}{dt} \|f\|_{L^2(d\mu_t)}^2 \le \frac{\alpha}{2} \left( \frac{\eta}{2} \alpha^3 \|E\|_{L^\infty}^2 + \alpha' \right) \iint_{\mathbb{T} \times \mathbb{R}} f^2 |v|^2 d\mu_t + \frac{1}{4\eta} \|f\|_{L^2(d\mu_t)}^2. $$

Si $\alpha$ vérifie

$$ (6.6) \qquad \frac{\alpha'(t)}{\alpha(t)^3} \le -\frac{\eta}{2} \|E(t)\|_{L^\infty}^2 \quad \forall t \ge 0, $$

alors pour tous $t \ge 0$,

$$ (6.7) \qquad \|f(t)\|_{L^2(d\mu_t)} \le \|f_0\|_{L^2(d\mu_0)} e^{t/4\eta}. $$

Dans la section 6.1, j'introduis la formulation de l'équation de Vlasov à l'aide d'une base de Hermite en vitesse, et je définis une classe de discrétisation de type Galerkin discontinue en espace, pour lesquelles les conservations sont discutées. La section 6.2 est dédiée aux propriétés de stabilité des discrétisations, en lien avec la définition de la fonction $\alpha$. Nous établissons finalement un résultat de convergence avec estimation d'erreur entre la solution approchée et la solution régulière du système de Vlasov-Poisson dans la section 6.3.
## 6.1 Méthode spectrale de Hermite en vitesse, Galerkin discontinue en espace

Nous nous plaçons dans le même cadre que [FX20] : l'approximation en vitesse s'effectue en développant la fonction de distribution à l'aide d'une base de Hermite, tandis que la discrétisation en espace de l'équation de Vlasov est traitée par une méthode de Galerkin discontinue. La méthode est construite en supposant la fonction strictement positive $t \mapsto \alpha(t)$ donnée, sa définition étant discutée dans la section 6.2. Les propriétés de conservation présentées dans cette section ne sont pas impactées par le choix de $\alpha$.

### 6.1.1 Forme spectrale de Hermite

Nous considérons les fonctions de base de Hermite suivantes, à poids asymétrique dépendant du temps :

$$ (6.8) \qquad \Psi_n(t, v) = \alpha(t) H_n(\alpha(t)v) \omega(t, v)^{-1}, $$

où $H_n$ sont les polynômes de Hermite définis par $H_{-1}(\xi) = 0$, $H_0(\xi) = 1$ et pour $n \ge 1$,

$$ \sqrt{n} H_n(\xi) = \xi H_{n-1}(\xi) - \sqrt{n-1} H_{n-2}(\xi). $$

L'ensemble $(\Psi_n)_n$ est un système orthogonal dans l'espace $L^2(\omega(t)dv)$ :

$$ (6.9) \qquad \langle \Psi_n, \Psi_m \rangle_{L^2(\omega(t)dv)} = \alpha(t) \int_{\mathbb{R}} \Psi_n(v) H_m(\alpha(t)v) dv = \alpha(t) \delta_{n,m}. $$

Notons finalement $V_N$ le sous-espace de $L^2(\omega(t)dv)$ défini par

$$ (6.10) \qquad V_N := \operatorname{Span}\{\Psi_n(t), \quad 0 \le n \le N-1\}. $$

Nous construisons une approximation $f_N$ de la solution $f$ de (6.1) sous la forme

$$ (6.11) \qquad f_N(t, x, v) = \sum_{n=0}^{N-1} C_n(t, x) \Psi_n(t, v), $$

où le système vérifié par les coefficients $(C_n)_{0 \le n \le N-1}$ est obtenu en utilisant la propriété d'orthogonalité (6.9) et en prenant $H_n(\alpha v)$ comme fonction test dans (6.1) : pour $n = 0, \dots, N-1$,

$$ (6.12) \qquad \left\{ \begin{array}{l} \partial_t C_n + \mathcal{T}_n[C] = S_n[C, E_N], \\ \mathcal{T}_n[C] = \frac{1}{\alpha} \left( \sqrt{n} \partial_x C_{n-1} + \sqrt{n+1} \partial_x C_{n+1} \right), \\ S_n[C, E_N] = \frac{\alpha'}{\alpha} \left( n C_n + \sqrt{(n-1)n} C_{n-2} \right) + E_N \alpha \sqrt{n} C_{n-1}, \end{array} \right. $$
avec $C_n = 0$ pour $n < 0$ et $n \ge N$. En remarquant que $\rho_N = C_0$, l'équation de Poisson est
alors approchée par

$$
(6.13) \qquad \frac{\partial E_N}{\partial x} = C_0 - \rho_{0,N}.
$$

Les équations pour $n = 0, 1, 2$ du système (6.12) permettent d'établir la conservation de la masse, de la quantité de mouvement et de l'énergie totale.

**Proposition 6.1.** Soit $f_N$ définie par (6.11), où $((C_n)_n, E_N)$ est solution du système de Vlasov-Poisson écrit sous la forme (6.12)–(6.13). Alors la masse et la quantité de mouvement sont conservées :

$$
\begin{align*}
\frac{d}{dt} \int_{\mathbb{T} \times \mathbb{R}} f_N \, dv \, dx &= \frac{d}{dt} \int_{\mathbb{T}} C_0 \, dx = 0, \\
\frac{d}{dt} \int_{\mathbb{T} \times \mathbb{R}} v \, f_N \, dv \, dx &= \frac{d}{dt} \int_{\mathbb{T}} \frac{C_1}{\alpha} \, dx = 0,
\end{align*}
$$

ainsi que l'énergie totale :

$$
\frac{d}{dt} \left( \frac{1}{2} \int_{\mathbb{T} \times \mathbb{R}} v^2 f_N \, dv \, dx + \int_{\mathbb{T}} |E_N|^2 \, dx \right) = \frac{d}{dt} \left( \frac{1}{2} \int_{\mathbb{T}} \left( \frac{1}{\alpha^2} (\sqrt{2} C_2 + C_0) + |E_N|^2 \right) \, dx \right) = 0.
$$

6.1.2 Discrétisation en espace

Nous considérons une approximation de Galerkin discontinue de l'équation de Vlasov
écrite sous la forme (6.12) grâce à la base de Hermite. Le tore $\mathbb{T}$ est subdivisé en $N_x$
sous-intervalles $I_j = [x_{j-1/2}, x_{j+1/2}]$, $j \in \mathcal{J} = \{0, \dots, N_x - 1\}$, de longueur $h = 1/N_x$.
Notons $\hat{\mathcal{J}} = \{0, \dots, N_x\}$ l'ensemble des indices des interfaces. Introduisons finalement le
paramètre de discrétisation en espace et vitesse $\delta = (h, 1/N)$.

Étant donné un degré $k \in \mathbb{N}$ fixé, définissons l'espace

$$
(6.14) \qquad X_h = \{u \in L^2(\mathbb{T}) : u|_{I_j} \in \mathbb{P}_k(I_j), \ j \in \mathcal{J}\},
$$

où $\mathbb{P}_k(I)$ désigne l'espace local des polynômes de degré au plus $k$ sur l'intervalle $I$. En
notant $u_{j+1/2}^{\pm} = u(x_{j+1/2}^{\pm})$, le saut $[u]_{j+1/2}$ et la moyenne $\{u\}_{j+1/2}$ de $u$ à l'interface
$x_{j+1/2}$ sont définis respectivement par

$$
[u]_{j+1/2} = u(x_{j+1/2}^{+}) - u(x_{j+1/2}^{-}) \quad \text{et} \quad \{u\}_{j+1/2} = \frac{1}{2} (u(x_{j+1/2}^{+}) + u(x_{j+1/2}^{-})), \quad \forall j \in \hat{\mathcal{J}}.
$$

Nous cherchons une approximation $C_{\delta} = (C_{\delta,n})_{0 \le n \le N-1}$ des modes $(C_n)_n$ telle que
$C_{\delta,n}(t, \cdot) \in X_h$ et pour tous $\varphi_n \in X_h$,
(6.15)

$$
\frac{d}{dt} \int_{I_j} C_{\delta,n} \varphi_n \, dx + A_{n,j}(g_n(C_{\delta}), \varphi_n) = \int_{I_j} S_n[C_{\delta}, E_{\delta}] \varphi_n \, dx, \quad j \in \mathcal{J}, \quad 0 \le n \le N-1,
$$
où $A_{n,j}$ est défini par

$$ (6.16) \qquad \left\{ \begin{array}{l} A_{n,j}(g_n, \varphi_n) = - \int_{I_j} g_n \varphi'_n dx + \hat{g}_{n,j+\frac{1}{2}} \varphi^-_{n,j+1/2} - \hat{g}_{n,j-1/2} \varphi^+_{n,j-1/2}, \\ g_n(C_\delta) = \frac{1}{\alpha} \left( \sqrt{n} C_{\delta,n-1} + \sqrt{n+1} C_{\delta,n+1} \right). \end{array} \right. $$

Le flux numérique $\hat{g}_n$ est donné par

$$ (6.17) \qquad \hat{g}_n = \frac{1}{2} \left[ g_n^-(C_\delta) + g_n^+(C_\delta) - \frac{\nu_n}{\alpha} (C_{\delta,n}^+ - C_{\delta,n}^-) \right], $$

avec le coefficient de viscosité numérique $\nu_n$ tel que $\nu_n \in [\underline{\nu}, \bar{\nu}]$ avec $0 < \underline{\nu} \le \bar{\nu} < \infty$.

Nous obtenons donc une solution approchée en espace et vitesse de (6.1) sous la forme

$$ (6.18) \qquad f_{\delta}(t, x, v) = \sum_{n=0}^{N-1} C_{\delta,n}(t, x) \Psi_{n}(t, v), $$

où $(C_{\delta,n})_n$ vérifie le système (6.15) et $(\Psi_n)_n$ sont les fonctions de base définies par (6.8).

Concernant l'approximation $E_\delta$ du champ électrique $E_N$, considérons le potentiel $\Phi_N$ tel que

$$ (6.19) \qquad E_N = -\frac{\partial \Phi_N}{\partial x}, \quad \frac{\partial E_N}{\partial x} = C_0 - \rho_{0,N}, $$

conduisant à l'équation de Poisson unidimensionnelle

$$ -\frac{\partial^2 \Phi_N}{\partial x^2} = C_0 - \rho_{0,N}. $$

Plusieurs choix de discrétisations de ce problème sont possibles.

* **Approximation conforme** : cela correspond à une intégration directe du problème (6.19), qui est évidente en dimension un.

* **Approximation Galerkin discontinue** : déterminer $\Phi_\delta(t, \cdot)$, $E_\delta(t, \cdot) \in X_h$ tels que pour tous $\eta$ et $\zeta \in X_h$, pour tous $i \in \mathcal{J}$,

$$ (6.20) \qquad \left\{ \begin{array}{l} + \int_{I_j} \Phi_\delta \eta' dx - \hat{\Phi}_{\delta,j+\frac{1}{2}} \eta^-_{j+\frac{1}{2}} + \hat{\Phi}_{\delta,j-\frac{1}{2}} \eta^+_{j-\frac{1}{2}} = \int_{I_j} E_\delta \eta dx, \\ - \int_{I_j} E_\delta \zeta' dx + \hat{E}_{\delta,j+\frac{1}{2}} \zeta^-_{j+\frac{1}{2}} - \hat{E}_{\delta,j-\frac{1}{2}} \zeta^+_{j-\frac{1}{2}} = \int_{I_j} (C_{\delta,0} - \rho_{0,\delta}) \zeta dx, \end{array} \right. $$

où les flux numériques $\hat{\Phi}_\delta$ et $\hat{E}_\delta$ sont donnés par

$$ (6.21) \qquad \left\{ \begin{array}{l} \hat{\Phi}_{\delta} = \{\Phi_{\delta}\}, \\ \hat{E}_{\delta} = \{E_{\delta}\} - \beta[\Phi_{\delta}], \end{array} \right. $$

avec $\beta > 0$ une constante ou une constante multipliée par $1/h$.
<figure><img src="image_30.png" /><figcaption>FIGURE 6.1 : Variation de la masse, de la quantité de mouvement et de l'énergie totale au cours du temps, pour l'instabilité *bump-on-tail*.</figcaption></figure>

• Éléments finis mixtes : en notant

$$Y_h := \{u \in C^0(\mathbb{T}) : u|_{I_i} \in \mathbb{P}_{k+1}(I_i), \forall i\},$$

déterminer un couple $(E_\delta, \Phi_\delta) \in Y_h \times X_h$ tel que pour tous $(\eta, \zeta) \in Y_h \times X_h$,

$$ (6.22) \qquad \left\{ \begin{array}{l} \displaystyle \int_{\mathbb{T}} \Phi_{\delta} \eta' \, dx = \int_{\mathbb{T}} E_{\delta} \eta \, dx, \\[1em] \displaystyle \int_{\mathbb{T}} \partial_x E_{\delta} \zeta \, dx = \int_{\mathbb{T}} (C_{\delta,0} - \rho_{0,\delta}) \zeta \, dx. \end{array} \right. $$

Quelle que soit l'approximation choisie pour le problème de Poisson, l'équation $n = 0$ du système (6.15) assure la conservation de la masse :

$$ \frac{d}{dt} \int_{\mathbb{T} \times \mathbb{R}} f_{\delta} \, dv \, dx = \frac{d}{dt} \int_{\mathbb{T}} C_{\delta,0} \, dx = 0. $$

Dans le cas d'une discrétisation Galerkin discontinue de l'équation de Poisson (6.20), avec $\nu_0 = 0$ dans (6.17) (correspondant à un flux centré pour le schéma sur $C_{\delta,0}$), la conservation de l'énergie totale est prouvée rigoureusement, mais pas celle de la quantité de mouvement. Cependant, comme illustré sur la figure 6.1 dans le cas *bump-on-tail* [10], nous constatons numériquement que la variation au cours du temps de la quantité de mouvement reste très faible.

## 6.2 Propriétés de stabilité

La définition (6.4) du poids $\omega(t, v)$, avec une fonction $t \mapsto \alpha(t)$ vérifiant un analogue de la condition (6.6), permet d'établir rigoureusement la stabilité $L^2(d\mu_t)$ des approximations $f_N$ et $f_\delta$, en adaptant les idées du cadre continu au cadre semi-discret.
Plus précisément, indépendamment du choix de discrétisation pour l'équation de Poisson, nous obtenons le résultat de stabilité suivant pour l'approximation $f_\delta$ obtenue par la méthode de Hermite / Galerkin discontinue.

**Proposition 6.2.** Soit $(f_\delta, E_\delta)$ la solution approchée de (6.1) définie par (6.15)-(6.18), avec la fonction $\alpha$ définie par

$$ (6.23) \qquad \alpha(t) := \alpha_0 \left( 1 + \gamma \alpha_0^2 \int_0^t \max(1, \|E_\delta(s)\|_{L^\infty}^2) ds \right)^{-1/2}. $$

Supposons que $\|f_\delta(0)\|_{L^2(d\mu_0)} < +\infty$. Alors, pour tous $t \ge 0$,

$$ (6.24) \qquad \begin{aligned} \frac{d}{dt} \|f_\delta(t)\|_{L^2(d\mu_t)}^2 &:= \frac{d}{dt} \left( \alpha(t) \sum_{n=0}^{N-1} \int_{\mathbb{T}} |C_{\delta,n}|^2 dx \right) \\ &\le - \sum_{n=0}^{N-1} \sum_{j \in \hat{\mathcal{J}}} \nu_n [C_{\delta,n}]_{j-\frac{1}{2}}^2 + \frac{1}{2\gamma} \|f_\delta(t)\|_{L^2(d\mu_t)}^2, \end{aligned} $$

d'où nous déduisons

$$ (6.25) \qquad \|f_\delta(t)\|_{L^2(d\mu_t)} \le \|f_\delta(0)\|_{L^2(d\mu_0)} e^{t/4\gamma}, $$

où le paramètre fixé $\gamma > 0$, pouvant être choisi arbitrairement, est celui apparaissant dans la définition (6.23) de $\alpha$.

Remarquons que l'estimation (6.24) est un analogue discret de (6.5), dans lequel apparaît un terme de dissipation additionnel provenant de la discrétisation Galerkin discontinue en espace.

Pour étudier la convergence de la méthode numérique, nous avons besoin d'une estimation de stabilité pour la solution exacte $f$ du système de Vlasov-Poisson (6.1) dans la même norme à poids $L^2(d\mu_t)$. Ce résultat n'est pas évident puisque le poids $\omega$ dépend de la solution approchée, du fait de la dépendance de la fonction $\alpha$ en $\|E_\delta\|_{L^\infty}^2$. En reprenant la preuve de (6.7) esquissée dans l'introduction de ce chapitre et en ajustant de manière appropriée le paramètre $\eta$ impliqué dans l'inégalité de Young, la proposition suivante peut être établie.

**Proposition 6.3.** Soit $(f, E)$ une solution régulière du système de Vlasov-Poisson (6.1). Supposons que $f_0 \in L^2(d\mu_0)$. Alors il existe $\eta > 0$, dépendant uniquement de la masse initiale de $f_0$, tel que pour tous $t \ge 0$,

$$ \|f(t)\|_{L^2(d\mu_t)} \le \|f_0\|_{L^2(d\mu_0)} e^{t/4\eta}, $$

où $d\mu_t$ est défini par (6.3), la fonction $\alpha$ apparaissant dans le poids $\omega$ étant donnée par (6.23).

Plus précisément, le paramètre $\eta > 0$ est choisi tel que $\eta C/\gamma < 1$, où $C > 0$ est une constante dépendant uniquement de la masse de $f_0$, telle que $\|E(t)\|_{L^\infty}^2 \le C$ pour tous $t \ge 0$, et $\gamma > 0$ est le paramètre fixé apparaissant dans la définition (6.23) de $\alpha$.
Notons que l'espace $L^2(d\mu_t)$ dépend de $\alpha$ défini par (6.23), et donc du paramètre de discrétisation $\delta$ du fait du terme $\|E_\delta\|_{L^\infty}$ impliqué dans cette définition. Ainsi, pour établir un résultat de convergence, il est indispensable de contrôler $\alpha$ uniformément par rapport à $\delta$. Ce contrôle est obtenu en bornant uniformément par rapport à $\delta$ la norme $L^\infty$ du champ électrique discret $E_\delta$.

**Proposition 6.4.** Soit $(f_\delta, E_\delta)$ la solution approchée de (6.1) définie par (6.15)-(6.18) et une approximation conforme de l'équation de Poisson. La fonction $\alpha$ est donnée par (6.23). Supposons $\|f_\delta(0)\|_{L^2(d\mu_0)} < +\infty$. Soit $T > 0$ un temps final fixé.

Alors il existe une constante $C_T > 0$, indépendante du paramètre de discrétisation $\delta$, telle que

$$ (6.26) \qquad \|E_\delta(t)\|_{L^\infty} \le C_T, \quad \forall t \in [0, T]. $$

De plus, la fonction $\alpha$ vérifie

$$ (6.27) \qquad 0 < \underline{\alpha}_T \le \alpha(t) \le \alpha_0, \quad \forall t \in [0, T], $$

où la constante $\underline{\alpha}_T$ est indépendante de $\delta$ et donnée par

$$ \underline{\alpha}_T := \alpha_0 (1 + \gamma \alpha_0^2 (1 + C_T) T)^{-1/2}. $$

Les inégalités de Sobolev et Poincaré-Wirtinger permettent d'établir qu'il existe une constante $c_0 > 0$ telle que

$$ (6.28) \qquad \|E_\delta(t)\|_{L^\infty}^2 \le \frac{c_0}{\alpha(t)} \|f_\delta(t)\|_{L^2(d\mu_t)}^2. $$

En utilisant l'estimation de stabilité (6.25) et la définition (6.23) de $\alpha$, le lemme de Gronwall permet d'obtenir la borne (6.26) sur la norme $L^\infty$ de $E_\delta$. Grâce à cette borne, l'obtention de la borne inférieure $\underline{\alpha}_T$ de $\alpha(t)$ est directe, et la borne supérieure $\alpha_0$ est claire par décroissance de $\alpha$.

Les résultats numériques que nous avons obtenus [10] indiquent que ce cadre mathématique avec un poids dépendant du temps est approprié pour étudier la stabilité. La figure 6.2 présente l'évolution au cours du temps de la norme $L^2$ à poids de $f_\delta$ et de la fonction $\alpha$ pour l'instabilité double faisceau. Les variations de $\alpha$ restent petites même si le champ électrique varie fortement du fait de l'instabilité. Dans le cas où $\alpha$ est choisi constant au cours du temps, nous constatons que la norme de $f_\delta$ croît plus rapidement, alors que le choix de $\alpha(t)$ permet d'amortir cette norme quand les effets non linéaires dominent. La comparaison avec une solution de référence obtenue par le schéma PFC [FSB01] sur un maillage fin permettent de valider nos résultats.

## 6.3 Convergence de la méthode

Dans cette section, je me place dans le cas d'une approximation conforme de l'équation de Poisson (6.19). À partir des résultats de stabilité obtenus, nous sommes en mesure d'établir la convergence de la semi-discrétisation Hermite / Galerkin discontinue proposée, avec une estimation d'erreur.
<figure><img src="image_31.png" /><figcaption>FIGURE 6.2 : Évolution de $||f_δ||_{L^2(dμ_t)}$ (gauche) et de la fonction $α_δ$ (droite), et comparaison avec une solution de référence obtenue avec le schéma PFC, pour l'instabilité double faisceau.</figcaption></figure>

**Théorème 6.1.** Pour tous $t \in [0, T]$, en considérant la fonction $α$ définie par (6.23), soit $f(t, \cdot) \in H^m(dμ_t)$ la solution du système de Vlasov-Poisson (6.1) où $m \ge k + 1$, et $f_δ$ l'approximation définie par (6.15)-(6.18). Alors il existe une constante $C > 0$ indépendante de $δ = (h, 1/N)$ mais dépendant du temps final $T$, telle que

$$ (6.29) \quad \|f(t) - f_\delta(t)\|_{L^2(dxdv)} \le \|f(t) - f_\delta(t)\|_{L^2(d\mu_t)} \le C \left[ \frac{1}{N^{(m-1)/2}} + h^{k+1/2} \right]. $$

La preuve du théorème 6.1 s'appuie classiquement sur l'étude de l'erreur totale décomposée comme la somme de l'erreur de projection et de l'erreur de consistance. Plus précisément, soit $W_δ$ le sous-espace de $L^2(dμ_t)$ défini par $W_δ := V_N \otimes X_h$, où $V_N$ est donné par (6.10) et $X_h$ est défini par (6.14). En notant $P_{W_δ}$ la projection orthogonale sur $W_δ$, l'erreur totale est estimée comme

$$ (6.30) \quad \|f(t) - f_\delta(t)\|_{L^2(d\mu_t)} \le \|f(t) - P_{W_\delta}f(t)\|_{L^2(d\mu_t)} + \|P_{W_\delta}f(t) - f_\delta(t)\|_{L^2(d\mu_t)}. $$

Le premier terme, erreur de projection, s'estime en combinant les propriétés de $P_{V_N}$ la projection sur $V_N$, obtenues en généralisant les résultats obtenus dans [FGT01] au cas où $α$ dépend du temps, et les propriétés classiques de $P_{X_h}$ la projection sur $X_h$. Le second terme, erreur de consistance, s'estime quant à lui en étudiant sa dérivée en temps, et en utilisant les arguments de stabilité sur $f$ et $f_δ$ ainsi que les propriétés d'interpolation.

Le théorème 6.1 met en évidence la précision spectrale de la discrétisation en vitesse basée sur des polynômes de Hermite, tandis que l'ordre de convergence classique de la méthode de Galerkin discontinue est retrouvée pour la discrétisation spatiale. Nous vérifions bien ce résultat en pratique, ce qui est mis en évidence dans la table 6.1, où nous présentons les résultats obtenus pour l'amortissement Landau avec $k = 1$ et $2$ [9].
<table><caption>TABLE 6.1 : Erreur $||f - f_δ||_{L^2(dμ_t)}$ et ordre de convergence pour l'amortissement Landau.</caption><thead><tr><th rowspan="2">N<sub>x</sub> × N</th><th colspan="2">P<sub>1</sub></th><th colspan="2">P<sub>2</sub></th></tr><tr><th>L<sup>2</sup>(dμ<sub>t</sub>) error</th><th>Order</th><th>L<sup>2</sup>(dμ<sub>t</sub>) error</th><th>Order</th></tr></thead><tbody><tr><td>16 × 16</td><td>5.12E-4</td><td>–</td><td>1.44E-5</td><td>–</td></tr><tr><td>32 × 32</td><td>1.05E-4</td><td>2.28</td><td>1.68E-6</td><td>3.09</td></tr><tr><td>64 × 64</td><td>2.31E-5</td><td>2.18</td><td>2.05E-7</td><td>3.04</td></tr><tr><td>128 × 128</td><td>5.42E-6</td><td>2.09</td><td>2.48E-8</td><td>3.04</td></tr></tbody></table>

## 6.4 Perspectives

**Discrétisation en temps.** Dans nos travaux, nous avons considéré uniquement une semi-discrétisation en espace-vitesse. Dans [10], nous proposons pour la mise en œuvre de notre schéma une discrétisation temporelle de type Runge-Kutta d'ordre 2 pour l'approximation (6.15)–(6.17) de l'équation de Vlasov, couplée avec une discrétisation de $\alpha$ comme solution de

$$ \alpha' = -\frac{\gamma}{2} \|E\|_{L^{\infty}}^2 \alpha^3. $$

Cette méthode généralise au cas d'un poids dépendant du temps la discrétisation temporelle définie dans [FX20]. Cependant, elle ne permet pas d'obtenir rigoureusement de propriétés de conservation et de stabilité. Une perspective naturelle est de construire un schéma temporel de type Crank-Nicolson pour $f$ et $\alpha$ adapté de sorte à pouvoir démontrer rigoureusement la stabilité pour une discrétisation complète du problème.

**Généralisation au cas multidimensionnel.** Pour l'instant, notre analyse se limite au cas unidimensionnel, dans lequel le contrôle du champ électrique $E$ est aisé du fait de l'injection de $H^1$ dans $L^\infty$. À l'avenir, nous voudrions adapter notre approche au cadre multidimensionnel en s'appuyant sur les idées développées dans [ACS12] pour le contrôle de la norme $L^\infty$ du champ électrique en dimensions 2 et 3.
[1] C. Berthon, M. Bessemoulin-Chatard, and M. Mathis. Numerical convergence rate for a diffusive limit of hyperbolic systems : *p*-system with damping. *SMAI J. Comput. Math.*, 2 :99–119, 2016. https://hal.archives-ouvertes.fr/hal-01360107/document.

[2] M. Bessemoulin-Chatard. A finite volume Scheme for convection-diffusion equations with nonlinear diffusion derived from the Scharfetter-Gummel scheme. *Numer. Math*, 121 :637–670, 2012.

[3] M. Bessemoulin-Chatard and C. Chainais-Hillairet. Exponential decay of a finite volume scheme to the thermal equilibrium for drift-diffusion systems. *J. Numer. Math.*, 25(3) :147–168, 2017. https://hal.archives-ouvertes.fr/hal-01250709v2/document.

[4] M. Bessemoulin-Chatard and C. Chainais-Hillairet. Uniform-in-time bounds for approximate solutions of the drift-diffusion system. *Numerische Mathematik*, 141(4) :881–916, 2019. https://hal.archives-ouvertes.fr/hal-01659418v2/document.

[5] M. Bessemoulin-Chatard, C. Chainais-Hillairet, and F. Filbet. On discrete functional inequalities for some finite volume schemes. *IMA Journal of Numerical Analysis*, 35(3) :1125–1149, 2015.

[6] M. Bessemoulin-Chatard, C. Chainais-Hillairet, and A. Jüngel. Uniform $L^\infty$ Estimates for Approximate Solutions of the Bipolar Drift-Diffusion System. In *International Conference on Finite Volumes for Complex Applications*, pages 381–389. Springer, 2017. https://hal.archives-ouvertes.fr/hal-01472643/document.

[7] M. Bessemoulin-Chatard, C. Chainais-Hillairet, and H. Mathis. Numerical schemes for semiconductors energy-transport models. In *Finite volumes for complex applications IX—methods, theoretical aspects, examples—FVCA 9, Bergen, Norway, June*
2020, volume 323 of Springer Proc. Math. Stat., pages 75–90. Springer, Cham, [2020]
©2020.

[8] M. Bessemoulin-Chatard, C. Chainais-Hillairet, and H. Mathis. Analysis of numerical schemes for semiconductor energy-transport models. Appl. Numer. Math., 168:143–169, 2021.

[9] M. Bessemoulin-Chatard and F. Filbet. On the convergence of discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system. 2022, https://hal.archives-ouvertes.fr/hal-03760711/document, soumis.

[10] M. Bessemoulin-Chatard and F. Filbet. On the stability of conservative discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system. J. Comput. Phys., 451:Paper No. 110881, 28, 2022.

[11] M. Bessemoulin-Chatard, M. Herda, and T. Rey. Hypocoercivity and diffusion limit of a finite volume scheme for linear kinetic equations. Math. Comp., 89(323):1093–1133, 2020.

[12] M. Bessemoulin-Chatard, G. Lissoni, and H. Mathis. Numerical analysis of DDFV schemes for semiconductors energy-transport models. Comput. Appl. Math., 41(1):Paper No. 15, 29, 2022.

[13] S. Bulteau, M. Badsi, C. Berthon, and M. Bessemoulin-Chatard. A fully well-balanced and asymptotic preserving scheme for the shallow-water equations with a generalized Manning friction source term. Calcolo, 58(4):Paper No. 41, 37, 2021.

[14] S. Bulteau, C. Berthon, and M. Bessemoulin-Chatard. An Asymptotic Preserving scheme for the shallow-water equations with Manning friction using viscous correction of the HLL scheme. 2019, https://hal.archives-ouvertes.fr/hal-02479081v1/document, soumis.

[15] S. Bulteau, C. Berthon, and M. Bessemoulin-Chatard. Convergence rate of an Asymptotic Preserving scheme for the diffusive limit of the $p$-system with damping. Communications in Mathematical Sciences, 17(6):1459–1486, 2019. https://hal.archives-ouvertes.fr/hal-01895417/document.

[16] M. Chatard. Asymptotic Behavior of the Scharfetter–Gummel Scheme for the Drift-Diffusion Model. In Jaroslav Fořt, Jiří Fürst, Jan Halama, Raphaèle Herbin, and Florence Hubert, editors, Finite Volumes for Complex Applications VI Problems & Perspectives, volume 4 of Springer Proceedings in Mathematics, pages 235–243. Springer Berlin Heidelberg, 2011.
[ABH07] B. Andreianov, F. Boyer, and F. Hubert. Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes. *Numer. Methods Partial Differential Equations*, 23(1):145-195, 2007.

[ACD$^{+}$04] A. Arnold, J. A. Carrillo, L. Desvillettes, J. Dolbeault, A. Jüngel, C. Lederman, P. A. Markowich, G. Toscani, and C. Villani. Entropies and equilibria of many-particle systems : an essay on recent research. *Monatsh. Math.*, 142(1-2):35-43, 2004.

[ACHFH] D. Abdel, C. Chainais-Hillairet, P. Farrell, and M. Herda. Numerical analysis of a finite volume scheme for charge transport in perovskite solar cells. 2022, https://hal.science/hal-03779230v2.

[ACM17] B. Andreianov, C. Cancès, and A. Moussa. A nonlinear time compactness result and applications to discretization of degenerate parabolic-elliptic PDEs. *J. Funct. Anal.*, 273(12):3633-3670, 2017.

[ACS12] B. Ayuso, J. A. Carrillo, and C.-W. Shu. Discontinuous Galerkin methods for the multi-dimensional Vlasov-Poisson problem. *Mathematical Models and Methods in Applied Sciences*, 22(12):1250042, 2012.

[ADLT21] L. Addala, J. Dolbeault, X. Li, and M. L. Tayeb. $L^2$-hypocoercivity and large time asymptotics of the linearized Vlasov-Poisson-Fokker-Planck system. *J. Stat. Phys.*, 184(1):Paper No. 4, 34, 2021.

[AMTU01] A. Arnold, P. Markowich, G. Toscani, and A. Unterreiter. On convex Sobolev Inequalities and the rate of convergence to equilibrium for Fokker-Planck type equations. *Comm. Partial Differential Equations*, 26:43-100, 2001.

[BC16] C. Berthon and C. Chalons. A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations. *Math. Comp.*, 85(299):1281-1307, 2016.
[BF] A. Blaustein and F. Filbet. On a discrete framework of hypocoercivity for kinetic equations. 2022, https://hal.science/hal-03792511.

[BGK54] P.L. Bhatnagar, E.P. Gross, and M. Krook. A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-component Systems. *Phys. Rev.*, 94(3):511–525, 1954.

[BHN07] S. Bianchini, B. Hanouzet, and R. Natalini. Asymptotic behavior of smooth solutions for partially dissipative hyperbolic systems with a convex entropy. *Comm. Pure Appl. Math.*, 60(11):1559–1622, 2007.

[Bia19] R. Bianchini. Strong convergence of a vector-BGK model to the incompressible Navier-Stokes equations via the relative entropy method. *J. Math. Pures Appl. (9)*, 132:280–307, 2019.

[Bia21] R. Bianchini. Relative entropy in diffusive relaxation for a class of discrete velocities BGK models. *Commun. Math. Sci.*, 19(1):39–54, 2021.

[BLP79] A. Bensoussan, J. L. Lions, and G. C. Papanicolaou. Boundary layers and homogenization of transport processes. *Publications of the Research Institute for Mathematical Sciences*, 15(1):53–157, 1979.

[BLT13] C. Berthon, P. LeFloch, and R. Turpault. Late-time/stiff-relaxation asymptotic-preserving approximations of hyperbolic equations. *Mathematics of Computation*, 82(282):831–860, 2013.

[BSS84] C. Bardos, R. Santos, and R. Sentis. Diffusion approximation and computation of the critical size. *Trans. Amer. Math. Soc.*, 284(2):617–649, 1984.

[BT11] C. Berthon and R. Turpault. Asymptotic preserving HLL schemes. *Numer. Methods Partial Differential Equations*, 27:1396–1422, 2011.

[CCC19] C. Calgaro, C. Colin, and E. Creusé. A combined finite volume - finite element scheme for a low-Mach system involving a Joule term. *AIMS Mathematics*, 5(1):311–331, 2019.

[CCHFG21] C. Cancès, C. Chainais-Hillairet, J. Fuhrmann, and B. Gaudeul. A numerical-analysis-focused comparison of several finite volume schemes for a unipolar degenerate drift-diffusion model. *IMA J. Numer. Anal.*, 41(1):271–314, 2021.

[CH09] C. Chainais-Hillairet. Discrete duality finite volume schemes for two-dimensional drift-diffusion and energy-transport models. *Internat. J. Numer. Methods Fluids*, 59(3):239–257, 2009.

[CH16] J.-S. Camier and F. Hermeline. A monotone nonlinear finite volume method for approximating diffusion operators on general meshes. *Internat. J. Numer. Methods Enggr.*, 107(6):496–519, 2016.
[CHH20] C. Chainais-Hillairet and M. Herda. Large-time behaviour of a family of finite volume schemes for boundary-driven convection-diffusion equations. *IMA J. Numer. Anal.*, 40(4):2473-2504, 2020.

[CJM+01] J. A. Carrillo, A. Jüngel, P. Markowich, G. Toscani, and A. Unterreiter. Entropy dissipation methods for degenerate parabolic problems and generalized Sobolev inequalities. *Monatsh. Math.*, 133:1-82, 2001.

[DGJ97] P. Degond, S. Génieys, and A. Jüngel. A system of parabolic equations in nonequilibrium thermodynamics including thermal and electrical effects. *J. Math. Pures Appl. (9)*, 76(10):991-1015, 1997.

[DGP00] P. Degond, T. Goudon, and F. Poupaud. Diffusion limit for nonhomogeneous and non-micro-reversible processes. *Indiana Univ. Math. J.*, 49(3):1175-1198, 2000.

[DJP00] P. Degond, A. Jüngel, and P. Pietra. Numerical discretization of energy-transport models for semiconductors with nonparabolic band structure. *SIAM J. Sci. Comput.*, 22(3):986-1007, 2000.

[DMS15] J. Dolbeault, C. Mouhot, and C. Schmeiser. Hypocoercivity for linear kinetic equations conserving mass. *Trans. Amer. Math. Soc.*, 367(6):3807-3828, 2015.

[DMTB15] A. Duran, F. Marche, R. Turpault, and C. Berthon. Asymptotic preserving scheme for the shallow water equations with source terms on unstructured meshes. *J. Comput. Phys.*, 287:184-206, 2015.

[EFG06] R. Eymard, J. Fuhrmann, and K. Gärtner. A finite volume scheme for nonlinear parabolic equations derived from one-dimensional local Dirichlet problems. *Numer. Math.*, 102(3):463-495, 2006.

[EGH00] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In *Handbook of numerical analysis*, volume VII, pages 713-1020. North-Holland, Amsterdam, 2000. 7

[FGT01] J. C. M. Fok, B. Guo, and T. Tang. Combined Hermite spectral-finite difference method for the Fokker-Planck equation. *Mathematics of computation*, 71(240):1497-1528, 2001.

[Fou02] M. Fournié. Numerical discretization of energy-transport model for semiconductors using high-order compact schemes. *Appl. Math. Lett.*, 15(6):721-726, 2002.

[FSB01] F. Filbet, E. Sonnendrücker, and P. Bertrand. Conservative numerical schemes for the Vlasov equation. *Journal of Computational Physics*, 172(1):166-187, 2001.
[FX20] F. Filbet and T. Xiong. Conservative Discontinuous Galerkin/Hermite Spectral Method for the Vlasov-Poisson System. Commun. Appl. Math. Comput., 2020.

[GG86] H. Gajewski and K. Gröger. On the basic equations for carrier transport in semiconductors. J. Math. Anal. Appl., 113:12-35, 1986.

[GG89] H. Gajewski and K. Gröger. Semiconductor equations for variable mobilities based on Boltzmann statistics or Fermi-Dirac statistics. Math. Nachr., 140:7-36, 1989.

[Gou05] T. Goudon. Hydrodynamic limit for the Vlasov-Poisson-Fokker-Planck system: Analysis of the two-dimensional case. Mathematical Models and Methods in Applied Sciences, 15(05):737-752, 2005.

[HJP03] S. Holst, A. Jüngel, and P. Pietra. A mixed finite-element discretization of the energy-transport model for semiconductors. SIAM J. Sci. Comput., 24(6):2058-2075, 2003.

[HJP04] S. Holst, A. Jüngel, and P. Pietra. An adaptive mixed scheme for energy-transport simulations of field-effect transistors. SIAM J. Sci. Comput., 25(5):1698-1716, 2004.

[HL92] L. Hsiao and T.-P. Liu. Convergence to nonlinear diffusion waves for solutions of a system of hyperbolic conservation laws with damping. Comm. Math. Phys., 143(3):599-605, 1992.

[HLvL83] A. Harten, P. D. Lax, and B. van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev., 25(1):35-61, 1983.

[HN04] F. Hérau and F. Nier. Isotropic hypoellipticity and trend to equilibrium for the Fokker-Planck equation with a high-degree potential. Arch. Ration. Mech. Anal., 171(2):151-218, 2004.

[Hol96] J. P. Holloway. Spectral velocity discretizations for the Vlasov-Maxwell equations. Transport theory and statistical physics, 25(1):1-32, 1996.

[Il'69] A.M. Il'in. A difference scheme for a differential equation with a small parameter multiplying the highest derivative. Math. Zametki, 6:237-248, 1969.

[JÖ1] A. Jüngel. Quasi-hydrodynamic semiconductor equations, volume 41 of Progress in Nonlinear Differential Equations and their Applications. Birkhäuser Verlag, Basel, 2001.

[Jin99] S. Jin. Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. SIAM J. Sci. Comput., 21(2):441-454 (electronic), 1999.
[JPT98] S. Jin, L. Pareschi, and G. Toscani. Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations. *SIAM J. Numer. Anal.*, 35(6):2405-2439, 1998.

[Jün95] A. Jüngel. Numerical approximation of a drift-diffusion model for semiconductors with nonlinear diffusion. *Z. Angew. Math. Mech.*, 75(10):783-799, 1995.

[Jün11] A. Jüngel. *Quasi-hydrodynamic semiconductor equations*, volume 41. Birhäuser, 2011.

[Jü94] A. Jüngel. On the existence and uniqueness of transient solutions of a degenerate nonlinear drift-diffusion model for semiconductors. *Math. Models Methods Appl. Sci.*, 4(5):677-703, 1994.

[Jü95] A. Jüngel. Qualitative behavior of solutions of a degenerate nonlinear drift-diffusion model for semiconductors. *Math. Models Methods Appl. Sci.*, 5(5):497-518, 1995.

[Kan20] M. Kantner. Generalized Scharfetter-Gummel schemes for electro-thermal transport in degenerate semiconductors using the Kelvin formula for the Seebeck coefficient. *J. Comput. Phys.*, 402:109091, 24, 2020.

[KRF+15] T. Koprucki, N. Rotundo, P. Farrell, D. Doan, and J. Fuhrmann. On thermodynamic consistency of a Scharfetter-Gummel scheme based on a modified thermal voltage for drift-diffusion equations with diffusion enhancement. *Optical and Quantum Electronics*, 47(6):1327-1332, 2015.

[LT13] C. Lattanzio and A.E. Tzavaras. Relative entropy in diffusive relaxation. *SIAM J. Math. Anal.*, 45(3):1563-1584, 2013.

[Man90] R. Manning. On the flow of water in open channels and pipes. *Transactions of the Institution of Civil Engineers of Ireland*, 20:161-207, 1890.

[MDBCF17] V. Michel-Dansac, C. Berthon, S. Clain, and F. Foucher. A well-balanced scheme for the shallow-water equations with topography or Manning friction. *J. Comput. Phys.*, 335:115-154, 2017.

[Mei10] M. Mei. Best asymptotic profile for hyperbolic $p$-system with damping. *SIAM J. Math. Anal.*, 42(1):1-23, 2010.

[MMS88] P. Marcati, A.J. Milani, and P. Secchi. Singular convergence of weak solutions for a quasilinear nonhomogeneous hyperbolic system. *Manuscripta Math.*, 60(1):49-69, 1988.

[MN06] C. Mouhot and L. Neumann. Quantitative perturbative study of convergence to equilibrium for collisional kinetic models in the torus. *Nonlinearity*, 19(4):969-998, 2006.