# Une classe de schémas conservatifs sur grilles décalées pour les systèmes hyperboliques

Michael Ndjinga, Katia Ait Ameur, Coraline Mounier

► To cite this version:

Michael Ndjinga, Katia Ait Ameur, Coraline Mounier. Une classe de schémas conservatifs sur grilles décalées pour les systèmes hyperboliques. SMAI 2021 - 10ème Biennale Française des Mathématiques Appliquées et Industrielles, Jun 2021, La Grande Motte, France. cea-04447208

HAL Id: cea-04447208

https://cea.hal.science/cea-04447208

Submitted on 8 Feb 2024

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_1.png" /></figure>

DE LA RECHERCHE À L'INDUSTRIE

# UNE CLASSE DE SCHÉMAS CONSERVATIFS SUR GRILLES DÉCALÉES POUR LES SYSTÈMES HYPERBOLIQUES

22 Juin 2021

M. Ndjinga, K. Ait Ameur, C. Mounier, CEA Saclay et X-CMAP

DES/ISAS/DM2S/STMF/LMEC
<figure><img src="image_1.png" /></figure>

DE LA RECHERCHE À L'INDUSTRIE

# UNE CLASSE DE SCHÉMAS CONSERVATIFS SUR GRILLES DÉCALÉES POUR LES SYSTÈMES HYPERBOLIQUES

22 Juin 2021

M. Ndjinga, K. Ait Ameur, C. Mounier, CEA Saclay et X-CMAP

DES/ISAS/DM2S/STMF/LMEC
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

Une nouvelle classe de schémas décalés

Résultats numériques

Conclusion et perspectives
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

Une nouvelle classe de schémas décalés

Résultats numériques

Conclusion et perspectives
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

Une nouvelle classe de schémas décalés

Résultats numériques

Conclusion et perspectives
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

Une nouvelle classe de schémas décalés

Résultats numériques

Conclusion et perspectives
## Il y a 10 ans ... à Fukushima

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

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

<figure><img src="image_4.png" /><figcaption>Lundi 21 mars 2011 - 12h00 UTC</figcaption></figure>

<figure><img src="image_5.png" /></figure>
# Il y a 10 ans ... à Fukushima

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

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

<figure><img src="image_4.png" /><figcaption>Lundi 21 mars 2011 - 12h00 UTC</figcaption></figure>

<figure><img src="image_5.png" /></figure>
* Multicomponent flow : liquid, vapor, solid obstacles

* Stiff source terms (heat source, boiling, singular pressure loss)

* Flow in subsonic or supersonic regime

- Note : define the Mach number for two phase flows

* Specifics of two phase flows

- non conservative models (e.g. 6 and 7 equation model for two-phase flows)

- characteristic fields neither GNL neither LD (e.g. 3 equation model with drift)

- not always admit entropies (e.g. 3 equation model with drift)

- sonic points in counter-current flows

**Conséquence : la structure mathématique et le comportement numérique sont TRES DIFFERENTS du cas des écoulements monophasiques**
* Multicomponent flow : liquid, vapor, solid obstacles

* Stiff source terms (heat source, boiling, singular pressure loss)

* Flow in subsonic or supersonic regime

- Note : define the Mach number for two phase flows

* Specifics of two phase flows

- non conservative models (e.g. 6 and 7 equation model for two-phase flows)

- characteristic fields neither GNL neither LD (e.g. 3 equation model with drift)

- not always admit entropies (e.g. 3 equation model with drift)

- sonic points in counter-current flows

**Conséquence : la structure mathématique et le comportement numérique sont TRES DIFFERENTS du cas des écoulements monophasiques**
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

Une nouvelle classe de schémas décalés

Résultats numériques

Conclusion et perspectives
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

Une nouvelle classe de schémas décalés

Résultats numériques

Conclusion et perspectives
<figure><img src="image_6.png" /></figure>

$$ \left\{ \begin{aligned} \frac{\partial}{\partial t}(\rho) &+ \nabla \cdot (\rho \mathbf{u}) & = 0, \\ \frac{\partial}{\partial t}(\rho \mathbf{u}) &+ \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla p & = \mu \Delta \mathbf{u} \\ \frac{\partial}{\partial t}(\rho E) &+ \nabla \cdot (\rho H \mathbf{u}) & = \lambda \Delta T \end{aligned} \right. \quad (1) $$

* $\rho$, $\mathbf{u}$, $p$ may have very different magnitudes

* Use reference length $L_{ref}$, velocity $u_{ref}$, pressure $p_{ref}$ to scale the equations

* Dimensionless numbers characterise the flow regimes

* Recover the low Mach asymptotic
<figure><img src="image_6.png" /></figure>

$$ \left\{ \begin{aligned} \frac{\partial}{\partial t}(\rho) &+ \nabla \cdot (\rho \mathbf{u}) & = 0, \\ \frac{\partial}{\partial t}(\rho \mathbf{u}) &+ \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla p & = \mu \Delta \mathbf{u} \\ \frac{\partial}{\partial t}(\rho E) &+ \nabla \cdot (\rho H \mathbf{u}) & = \lambda \Delta T \end{aligned} \right. \quad (1) $$

* $\rho$, $\mathbf{u}$, $p$ may have very different magnitudes

* Use reference length $L_{ref}$, velocity $u_{ref}$, pressure $p_{ref}$ to scale the equations

* Dimensionless numbers characterise the flow regimes

* Recover the low Mach asymptotic
New variables $\rho := \frac{\rho}{\rho_{ref}}$, $p := \frac{\rho}{\rho_{ref}}$, $\mathbf{u} := \frac{\mathbf{u}}{u_{ref}}$

$$ \left\{ \begin{aligned} \frac{\partial}{\partial t}(\rho) &+ \nabla \cdot (\rho \mathbf{u}) & = 0, \\ \frac{\partial}{\partial t}(\rho \mathbf{u}) &+ \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \frac{1}{M^2} \nabla p & = \frac{1}{\text{Re}} \Delta \mathbf{u} \\ \frac{\partial}{\partial t}(\rho E) &+ \nabla \cdot (\rho H \mathbf{u}) & = \frac{1}{\text{Pe}} \Delta T \end{aligned} \right. \quad (2) $$

$$ M = \frac{u_{ref}}{c_{ref}}, \quad \text{Re} = \frac{L_{ref} u_{ref}}{\frac{\mu}{\rho_{ref}}} \quad \text{Pe} = \frac{L_{ref} u_{ref}}{\frac{\lambda}{\rho_{ref} c_{p}}} $$

* In a nuclear reactor, $\text{Re} \gg 1$, $\text{Pe} \gg 1$, $M \ll 1$ ($\|\mathbf{u}\| \ll c$)

* We need more precision on $p$ than on other variables

Note : the Mach number is not clearly defined in two-phase flow
New variables $\rho := \frac{\rho}{\rho_{ref}}$, $\boldsymbol{p} := \frac{\rho}{\rho_{ref}}$, $\boldsymbol{u} := \frac{\boldsymbol{u}}{u_{ref}}$

$$ \left\{ \begin{aligned} \frac{\partial}{\partial t}(\rho) &+ \nabla \cdot (\rho \mathbf{u}) & = 0, \\ \frac{\partial}{\partial t}(\rho \mathbf{u}) &+ \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \frac{1}{M^2} \nabla p & = \frac{1}{\mathrm{Re}} \Delta \mathbf{u} \\ \frac{\partial}{\partial t}(\rho E) &+ \nabla \cdot (\rho H \mathbf{u}) & = \frac{1}{\mathrm{Pe}} \Delta T \end{aligned} \right. \quad (2) $$

$$ M = \frac{u_{ref}}{c_{ref}}, \quad \mathrm{Re} = \frac{L_{ref} u_{ref}}{\frac{\mu}{\rho_{ref}}} \quad \mathrm{Pe} = \frac{L_{ref} u_{ref}}{\frac{\lambda}{\rho_{ref} c_{p}}} $$

* In a nuclear reactor, $\mathrm{Re} \gg 1$, $\mathrm{Pe} \gg 1$, $M \ll 1$ ($||\mathbf{u}|| \ll c$)

* We need more precision on $\rho$ than on other variables

Note : the Mach number is not clearly defined in two-phase flow
For well prepared initial data ($M \ll 1$ initially)

$$ egin{cases} 
\rho(\vec{x}, t=0) = \rho_* + O(M)(\vec{x}) \\ 
\rho(\vec{x}, t=0) = p_* + O(M^2)(\vec{x}) \\ 
u(\vec{x}, t=0) = u_*(\vec{x}) + O(M)(\vec{x}) 
\end{cases} $$

the solution $(\rho, p, u)(t, \vec{x})$ to the compressible Euler equations takes the form

$$ egin{cases} 
\rho(\vec{x}, t) = \rho_* + O(M)(\vec{x}, t) \\ 
p(\vec{x}, t) = p_* + O(M^2)(\vec{x}, t) \\ 
u(\vec{x}, t) = \tilde{u}(\vec{x}, t) + O(M)(\vec{x}, t) 
\end{cases} $$

where $\tilde{u}$ is the solution of the incompressible Euler equations with initial data $\tilde{u}(\vec{x}, 0) = u_*(\vec{x})$

► Numerical schemes should preserve this asymptotics
For well prepared initial data ($M \ll 1$ initially)

$$ egin{cases} 
\rho(\vec{x}, t=0) = \rho_* + O(M)(\vec{x}) \\ 
p(\vec{x}, t=0) = p_* + O(M^2)(\vec{x}) \\ 
u(\vec{x}, t=0) = u_*(\vec{x}) + O(M)(\vec{x}) 
\end{cases} $$

the solution $(\rho, p, u)(t, \vec{x})$ to the compressible Euler equations takes the form

$$ egin{cases} 
\rho(\vec{x}, t) = \rho_* + O(M)(\vec{x}, t) \\ 
p(\vec{x}, t) = p_* + O(M^2)(\vec{x}, t) \\ 
u(\vec{x}, t) = \tilde{u}(\vec{x}, t) + O(M)(\vec{x}, t) 
\end{cases} $$

where $\tilde{u}$ is the solution of the incompressible Euler equations with initial data $\tilde{u}(\vec{x}, 0) = u_*(\vec{x})$

► Numerical schemes should preserve this asymptotics
► Upwind based FV colocated methods (Godunov, Roe, ...)

- do not satisfy the low Mach asymptotics X

► Centered FV colocated methods

- satisfy the low Mach asymptotics... ✓

<figure><img src="image_7.png" /><figcaption>Figure: Implicit **upwind** scheme (Lid driven cavity)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure: Implicit **centered** scheme (Lid driven cavity)</figcaption></figure>

- but generate spurious **checkerboard** oscillations X
* Upwind based FV colocated methods (Godunov, Roe, ...)

- do not satisfy the low Mach asymptotics X

* Centered FV colocated methods

- satisfy the low Mach asymptotics... ✓

<figure><img src="image_7.png" /><figcaption>Figure: Implicit upwind scheme (Lid driven cavity)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure: Implicit centered scheme (Lid driven cavity)</figcaption></figure>

- but generate spurious checkerboard oscillations X
* Modèle diphasique à 4 équations

* Canal 1D Chauffé

* Termes sources raides : puissance de chauffage, changement de phase

* Comparaison des schémas VF colocalisés avec et sans traitement spécifique du terme source

$$ 
\left\{
\begin{aligned}
& \frac{\partial(\phi\rho)}{\partial t} + \nabla \cdot (\phi\rho \vec{V}) = 0, \\
& \frac{\partial(\phi\rho C)}{\partial t} + \nabla \cdot (\phi\rho C\vec{V} + \phi\rho C(1-C)\vec{V}_r) = \nabla \cdot (\phi K_{cv} \nabla C) + \phi \Gamma_v, \\
& \frac{\partial(\phi\rho \vec{V})}{\partial t} + \nabla \cdot (\phi\rho \vec{V} \otimes \vec{V} + \phi\rho C(1-C)\vec{V}_r \otimes \vec{V}_r + \phi P I_d) = P \nabla\phi + \nabla \cdot (\phi \bar{\sigma}) + \phi \vec{\tau} + \phi\rho \vec{g} \\
& \frac{\partial(\phi\rho E)}{\partial t} + \nabla \cdot (\phi\rho H \vec{V} + \phi\rho C(1-C)(H_v - H_l)\vec{V}_r) = \phi\rho \vec{g} \cdot \vec{V} + \phi Q + \nabla \cdot (\phi \vec{d})
\end{aligned}
\right.
\quad (3)
$$
▶ Modèle diphasique à 4 équations

▶ Canal 1D Chauffé

▶ Termes sources raides : puissance de chauffage, changement de phase

▶ Comparaison des schémas VF colocalisés avec et sans traitement spécifique du terme source

$$
\left\{
\begin{array}{l}
\displaystyle \frac{\partial(\phi \rho)}{\partial t} + \nabla \cdot (\phi \rho \vec{V}) = 0, \\[1.5em]
\displaystyle \frac{\partial(\phi \rho C)}{\partial t} + \nabla \cdot (\phi \rho C \vec{V} + \phi \rho C(1-C) \vec{V}_r) = \nabla \cdot (\phi K_{cv} \nabla C) + \phi \Gamma_v, \\[1.5em]
\displaystyle \frac{\partial(\phi \rho \vec{V})}{\partial t} + \nabla \cdot (\phi \rho \vec{V} \otimes \vec{V} + \phi \rho C(1-C) \vec{V}_r \otimes \vec{V}_r + \phi P \vec{I}_d) = P \nabla \phi + \nabla \cdot (\phi \bar{\sigma}) + \phi \vec{\tau} + \phi \rho \vec{g} \\[1.5em]
\displaystyle \frac{\partial(\phi \rho E)}{\partial t} + \nabla \cdot (\phi \rho H \vec{V} + \phi \rho C(1-C) (H_v - H_l) \vec{V}_r) = \phi \rho \vec{g} \cdot \vec{V} + \phi Q + \nabla \cdot (\phi \vec{d})
\end{array}
\right.
\qquad (3)
$$
<figure><img src="image_9.png" /></figure>

► La pression doit être décroissante

► Schéma upwind → Oscillations parasites X

► Schéma centré → Oscillations parasites X

► Schéma upwind + traitement spécifique du terme source ✓
<figure><img src="image_9.png" /></figure>

► La pression doit être décroissante

► Schéma upwind → Oscillations parasites X

► Schéma centré → Oscillations parasites X

► Schéma upwind + traitement spécifique du terme source ✓
Observation empirique :

* **Solveurs VF colocalisés** → oscillations parasites et/ou pb de précision à faible nombre de Mach

- Solveurs de Riemann basés sur la théorie des sytèmes hyperboliques de lois de conservation

- Solveurs de pression colocalisés inspirés du cas incompressible (méthodes de projection)

- Oscillations d'autant plus importantes qu'un terme source est présent

* **Solveurs VF décalés** → plus robustes aux termes sources mais

- Souvent basé sur des formulations non conservatives

- Extension aux maillages généraux délicate

Nécessité d'un cadre théorique permettant la construction et l'analyse de schémas décalés stables.
Observation empirique :

► **Solveurs VF colocalisés** → oscillations parasites et/ou pb de précision à faible nombre de Mach

- Solveurs de Riemann basés sur la théorie des sytèmes hyperboliques de lois de conservation

- Solveurs de pression colocalisés inspirés du cas incompressible (méthodes de projection)

- Oscillations d'autant plus importantes qu'un terme source est présent

► **Solveurs VF décalés** → plus robustes aux termes sources mais

- Souvent basé sur des formulations non conservatives

- Extension aux maillages généraux délicate

Nécessité d'un cadre théorique permettant la construction et l'analyse de schémas décalés stables.
Contexte

Problématique bas Mach

**Diffusion numérique de schémas classiques**

Une nouvelle classe de schémas décalés

Résultats numériques

Conclusion et perspectives
Contexte

Problématique bas Mach

**Diffusion numérique de schémas classiques**

Une nouvelle classe de schémas décalés

Résultats numériques

Conclusion et perspectives
Soit une EDP non linéaire

$$ \partial_t U + A(U)\partial_x U = 0. \quad (4) $$

Pour une semi-discrétisation *disc* à 2 points donnée

$$ \partial_t U_i + A_{disc}^{+} \frac{U_i - U_{i-1}}{\Delta x} + A_{disc}^{-} \frac{U_{i+1} - U_i}{\Delta x} = 0. \quad (5) $$

on trouve lorsque $\Delta x \to 0$ :

$$ \partial_t U + A \partial_x U = D(U, \Delta x) + o(\Delta x) \quad (6) $$

$$ = -\Delta x D_{disc} \partial_{xx} U + o(\Delta x). \quad (7) $$

où

* $A = A_{disc}^{+} + A_{disc}^{-}$ (consistance de (5) avec (4))

* $D = A_{disc}^{+} - A_{disc}^{-}$ représente la diffusion/dissipation numérique → stabilité et précision du schéma.
Soit une EDP non linéaire

$$ \partial_t U + A(U)\partial_x U = 0. \quad (4) $$

Pour une semi-discrétisation *disc* à 2 points donnée

$$ \partial_t U_i + A_{disc}^{+} \frac{U_i - U_{i-1}}{\Delta x} + A_{disc}^{-} \frac{U_{i+1} - U_i}{\Delta x} = 0. \quad (5) $$

on trouve lorsque $\Delta x \to 0$ :

$$ \partial_t U + A \partial_x U = D(U, \Delta x) + o(\Delta x) \quad (6) $$

$$ = -\Delta x D_{disc} \partial_{xx} U + o(\Delta x). \quad (7) $$

où

* $A = A_{disc}^{+} + A_{disc}^{-}$ (consistance de (5) avec (4))

* $D = A_{disc}^{+} - A_{disc}^{-}$ représente la diffusion/dissipation numérique → stabilité et précision du schéma.
On néglige le terme $q^2/\rho = \rho u^2$ dans les équations d'Euler:

$$ egin{cases} \partial_t \rho + \partial_x q = 0 \\ \partial_t q + c^2 \partial_x \rho = 0 \end{cases} $$

C'est un système hyperbolique linéaire de valeurs propres $\pm c$.

Discrétisation centrée :

$$ egin{cases} \partial_t \rho_i + \frac{q_{i+1} - q_{i-1}}{2\Delta x} = 0 \\ \partial_t q_i + c^2 \frac{\rho_{i+1} - \rho_{i-1}}{2\Delta x} = 0 \end{cases} $$

* $D_{centre} = 0$ (schéma d'ordre 2)

* $D_{centre}$ est diagonalisable

* Découplage pair/impair des cellules
On néglige le terme $q^2/\rho = \rho u^2$ dans les équations d'Euler:

$$ \begin{cases} \partial_t \rho + \partial_x q = 0 \\ \partial_t q + c^2 \partial_x \rho = 0 \end{cases} $$

C'est un système hyperbolique linéaire de valeurs propres $\pm c$.

Discrétisation centrée :

$$ \begin{cases} \partial_t \rho_i + \frac{q_{i+1} - q_{i-1}}{2\Delta x} = 0 \\ \partial_t q_i + c^2 \frac{\rho_{i+1} - \rho_{i-1}}{2\Delta x} = 0 \end{cases} $$

* $D_{centre} = 0$ (schéma d'ordre 2)

* $D_{centre}$ est diagonalisable

* Découplage pair/impair des cellules
## Discrétisation upwind

$$ egin{cases} 
\partial_t \rho_i + \frac{q_{i+1} - q_{i-1}}{2\Delta x} = c \frac{\rho_{i+1} - 2\rho_i + \rho_{i-1}}{\Delta x} \\ 
\partial_t q_i + c^2 \frac{\rho_{i+1} - \rho_{i-1}}{2\Delta x} = c \frac{q_{i+1} - 2q_i + q_{i-1}}{\Delta x} 
\end{cases} $$

Matrice de diffusion numérique associée :

$$ D_{\text{upwind}} = |A_{\text{wave}}| = \begin{pmatrix} c & 0 \\ 0 & c \end{pmatrix}, \quad A_{\text{wave}} = \begin{pmatrix} 0 & 1 \\ c^2 & 0 \end{pmatrix}. $$

* On applique la même dissipation numérique aux deux équations.

* la matrice de diffusion numérique $D_{\text{upwind}}$ est diagonalisable sur $\mathbb{R}$
## Discrétisation upwind

$$ egin{cases} \partial_t \rho_i + \frac{q_{i+1} - q_{i-1}}{2\Delta x} = c \frac{\rho_{i+1} - 2\rho_i + \rho_{i-1}}{\Delta x} \\ \partial_t q_i + c^2 \frac{\rho_{i+1} - \rho_{i-1}}{2\Delta x} = c \frac{q_{i+1} - 2q_i + q_{i-1}}{\Delta x} \end{cases} $$

Matrice de diffusion numérique associée :

$$ D_{upwind} = |A_{wave}| = \begin{pmatrix} c & 0 \\ 0 & c \end{pmatrix}, \quad A_{wave} = \begin{pmatrix} 0 & 1 \\ c^2 & 0 \end{pmatrix}. $$

* On applique la même dissipation numérique aux deux équations.

* la matrice de diffusion numérique $D_{upwind}$ est diagonalisable sur $\mathbb{R}$
### Discrétisation décalée “simple”

$$ egin{cases} \partial_t \rho_i + \frac{q_{i+1/2} - q_{i-1/2}}{\Delta x} = 0 \\ \partial_t q_{i+1/2} + c^2 \frac{\rho_{i+1} - \rho_i}{\Delta x} = 0 \end{cases}, \quad U_i = \begin{pmatrix} \rho_i \\ q_{i+1/2} \end{pmatrix}. $$

Matrice de diffusion numérique associée :

$$ D_{stag} = \begin{pmatrix} 0 & 1 \\ -c^2 & 0 \end{pmatrix}, \quad A_{wave} = \begin{pmatrix} 0 & 1 \\ c^2 & 0 \end{pmatrix}. $$

* L'essentiel de la diffusion numérique porte sur la deuxième équation.

* la matrice de diffusion numérique $D_{stag}$ n'est pas diagonalisable sur $\mathbb{R}$
### Discrétisation décalée “simple”

$$ \left\{ \begin{array}{l} \partial_t \rho_i + \frac{q_{i+1/2} - q_{i-1/2}}{\Delta x} = 0 \\ \partial_t q_{i+1/2} + c^2 \frac{\rho_{i+1} - \rho_i}{\Delta x} = 0 \end{array} \right. , \quad U_i = \begin{pmatrix} \rho_i \\ q_{i+1/2} \end{pmatrix}. $$

Matrice de diffusion numérique associée :

$$ D_{stag} = \begin{pmatrix} 0 & 1 \\ -c^2 & 0 \end{pmatrix}, \quad A_{wave} = \begin{pmatrix} 0 & 1 \\ c^2 & 0 \end{pmatrix}. $$

* L'essentiel de la diffusion numérique porte sur la deuxième équation.

* la matrice de diffusion numérique $D_{stag}$ n'est pas diagonalisable sur $\mathbb{R}$
## Discrétisation décalée de la littérature (plus complexe)

$$ egin{cases} \partial_t \rho_i + \frac{1}{\Delta x} (\rho_{i+\frac{1}{2}}^{up} u_{i+\frac{1}{2}} - \rho_{i-\frac{1}{2}}^{up} u_{i-\frac{1}{2}}) = 0 \\ \partial_t (\bar{\rho}_{i+\frac{1}{2}} u_{i+\frac{1}{2}}) + c^2 \frac{1}{\Delta x} (\rho_{i+1} - \rho_i) = 0, \end{cases} \quad (8) $$

avec

$$ egin{cases} \rho_{i+\frac{1}{2}}^{up} = \frac{\rho_i + \rho_{i+1}}{2} + \text{sign}(u_{i+\frac{1}{2}}) \frac{\rho_i - \rho_{i+1}}{2}, \\ \bar{\rho}_{i+\frac{1}{2}} = \frac{1}{2}(\rho_i + \rho_{i+1}). \end{cases} \quad (9) $$

$$ \mathcal{D}(U, \Delta x) = \Delta x \begin{pmatrix} \operatorname{sign}(u) & 0 \\ 0 & 1 \end{pmatrix} \partial_x \left[ \begin{pmatrix} u & 0 \\ c^2 & 0 \end{pmatrix} \partial_x \begin{pmatrix} \rho \\ q \end{pmatrix} \right]. \quad (10) $$

* Si $\operatorname{sign}(u)$ est constant le schéma est linéarisable et stable

* Si $\operatorname{sign}(u)$ varie le schéma n'est pas linéarisable
## Discrétisation décalée de la littérature (plus complexe)

$$ egin{cases} \partial_t \rho_i + \frac{1}{\Delta x} (\rho_{i+\frac{1}{2}}^{up} u_{i+\frac{1}{2}} - \rho_{i-\frac{1}{2}}^{up} u_{i-\frac{1}{2}}) = 0 \\ \partial_t (\bar{\rho}_{i+\frac{1}{2}} u_{i+\frac{1}{2}}) + c^2 \frac{1}{\Delta x} (\rho_{i+1} - \rho_i) = 0, \end{cases} \quad (8) $$

avec

$$ egin{cases} \rho_{i+\frac{1}{2}}^{up} = \frac{\rho_i + \rho_{i+1}}{2} + \text{sign}(u_{i+\frac{1}{2}}) \frac{\rho_i - \rho_{i+1}}{2}, \\ \bar{\rho}_{i+\frac{1}{2}} = \frac{1}{2}(\rho_i + \rho_{i+1}). \end{cases} \quad (9) $$

$$ \mathcal{D}(U, \Delta x) = \Delta x \begin{pmatrix} \operatorname{sign}(u) & 0 \\ 0 & 1 \end{pmatrix} \partial_x \left[ \begin{pmatrix} u & 0 \\ c^2 & 0 \end{pmatrix} \partial_x \begin{pmatrix} \rho \\ q \end{pmatrix} \right]. \quad (10) $$

* Si $\operatorname{sign}(u)$ est constant le schéma est linéarisable et stable

* Si $\operatorname{sign}(u)$ varie le schéma n'est pas linéarisable
Pour un **système conservatif** d'EDPs

$$ \partial_t U + \partial_x F(U) = 0 \quad \text{forme conservative} $$

un **schéma conservatif** prend la forme

$$ \partial_t U_i + \frac{F_{i,j+1} - F_{i-1,j}}{\Delta x} = 0 \quad (11) $$

On peut écrire les **flux numériques** $F_{ij}$ comme la somme d'une partie centrée et d'un partie décentrée

$$ F_{i,j} = \frac{F(U_i) + F(U_j)}{2} + D \frac{U_i - U_j}{2}, \quad (12) $$

Lorsque $\Delta x \to 0$, (11) donne

$$ \partial_t U + \partial_x F(U) = \Delta x D \partial_{xx} U + \mathcal{O}((\Delta x)^2). $$

► Pour un opérateur $D$ donné, il existe un schéma conservatif ayant pour diffusion numérique $D$
Exemple 5 : cas d'un schéma conservatif

Pour un **système conservatif** d'EDPs

$$ \partial_t U + \partial_x F(U) = 0 \quad \text{forme conservative} $$

un **schéma conservatif** prend la forme

$$ \partial_t U_i + \frac{F_{i,j+1} - F_{i-1,j}}{\Delta x} = 0 \quad (11) $$

On peut écrire les **flux numériques** $F_{ij}$ comme la somme d'une partie centrée et d'un partie décentrée

$$ F_{i,j} = \frac{F(U_i) + F(U_j)}{2} + D \frac{U_i - U_j}{2}, \quad (12) $$

Lorsque $\Delta x \to 0$, (11) donne

$$ \partial_t U + \partial_x F(U) = \Delta x D \partial_{xx} U + \mathcal{O}((\Delta x)^2). $$

► Pour un opérateur $D$ donné, il existe un schéma conservatif ayant pour diffusion numérique $D$
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

**Une nouvelle classe de schémas décalés**

Résultats numériques

Conclusion et perspectives
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

**Une nouvelle classe de schémas décalés**

Résultats numériques

Conclusion et perspectives
Après linéarisation, on s'intéresse au système hyperbolique

$$ \partial_t U + A\partial_x U = 0, \quad F(U) = AU $$

A doit être diagonalisable ⇒ symétrisation possible

$$ S = P^{-1}AP, \quad {}^tS = S $$

* changement de variables $\bar{U} = PU$ :

$$ \partial_t \bar{U} + S\partial_x \bar{U} = 0 $$

* En variables “symétrisantes”, l'énergie reste constante

$$ \partial_t \int_{\mathbb{R}} ||\bar{U}||^2 dx = 0. $$

* D'où la stabilité en continu pour un système hyperbolique:

$$ ||\bar{U}(\cdot, t)||_{L^2} = ||\bar{U}_0||_{L^2}. $$
Après linéarisation, on s'intéresse au système hyperbolique

$$ \partial_t U + A\partial_x U = 0, \quad F(U) = AU $$

A doit être diagonalisable ⇒ symétrisation possible

$$ S = P^{-1}AP, \quad {}^tS = S $$

► changement de variables $\bar{U} = PU$ :

$$ \partial_t \bar{U} + S\partial_x \bar{U} = 0 $$

► En variables “symétrisantes”, l'énergie reste constante

$$ \partial_t \int_{\mathbb{R}} ||\bar{U}||^2 dx = 0. $$

► D'où la stabilité en continu pour un système hyperbolique:

$$ ||\bar{U}(\cdot, t)||_{L^2} = ||\bar{U}_0||_{L^2}. $$
1 On considère un schéma de diffusion numérique D

2 On calcule l'évolution de la norme dans la base symétrisante

$$ \partial_t \int_{\mathbb{R}} ||\bar{U}_h||^2 dx = -h \int t (\partial_x \bar{U}_h)(\bar{D} + t \bar{D})(\partial_x \bar{U}_h) + o(h). $$

3 On conclut que

- le schéma sera stable si

$$ \bar{D} + t \bar{D} \geq 0 $$

$$ \Rightarrow \| \bar{U}_h(\cdot, t) \|_{L^2} \le \| \bar{U}_{h0} \|_{L^2}. \quad (13) $$

- le taux de dissipation est lié au spectre de $\bar{D} + t \bar{D}$
1 On considère un schéma de diffusion numérique D

2 On calcule l'évolution de la norme dans la base symétrisante

$$ \partial_t \int_{\mathbb{R}} ||\bar{U}_h||^2 dx = -h \int t (\partial_x \bar{U}_h)(\bar{D} + t \bar{D})(\partial_x \bar{U}_h) + o(h). $$

3 On conclut que

- le schéma sera stable si

$$ \bar{D} + t \bar{D} \geq 0 $$

$$ \Rightarrow \| \bar{U}_h(\cdot, t) \|_{L^2} \le \| \bar{U}_{h0} \|_{L^2}. \quad (13) $$

- le taux de dissipation est lié au spectre de $\bar{D} + t \bar{D}$
► Le schéma upwind :

$$
D_{\text{upwind}} = c \mathbb{I}_d, \quad \bar{D}_{\text{upwind}} = c \mathbb{I}_d, \quad \bar{D}_{\text{upwind}} + {}^t \bar{D}_{\text{upwind}} = c \mathbb{I}_d.
$$

- Fort couplage entre cellules ✓

- Diffusion numérique trop forte si $c \gg \|\mathbf{u}\|_X$

► Le schéma centré :

$$
D_{\text{centered}} = 0, \quad \bar{D}_{\text{centered}} = 0, \quad \bar{D}_{\text{centered}} + {}^{t}\bar{D}_{\text{centered}} = 0.
$$

- Diffusion numérique nulle ✓

- Découplage des cellules paires et impaires x

► Le schéma décalé :

$$
D_{\text{stag}} = \begin{pmatrix} 0 & 1 \\ -c^2 & 0 \end{pmatrix}, \bar{D}_{\text{stag}} = \begin{pmatrix} 0 & c \\ -c & 0 \end{pmatrix}, \bar{D}_{\text{stag}} + {}^t \bar{D}_{\text{stag}} = 0.
$$

- Diffusion numérique nulle ✓

- Fort couplage entre cellules ✓
► Le schéma upwind :

$$
D_{\text{upwind}} = c \mathbb{I}_d, \quad \bar{D}_{\text{upwind}} = c \mathbb{I}_d, \quad \bar{D}_{\text{upwind}} + t \bar{D}_{\text{upwind}} = c \mathbb{I}_d.
$$

- Fort couplage entre cellules ✓

- Diffusion numérique trop forte si $c \gg \|\mathbf{u}\|_X$

► Le schéma centré :

$$
D_{\text{centered}} = 0, \quad \bar{D}_{\text{centered}} = 0, \quad \bar{D}_{\text{centered}} + t \bar{D}_{\text{centered}} = 0.
$$

- Diffusion numérique nulle ✓

- Découplage des cellules paires et impaires x

► Le schéma décalé :

$$
D_{\text{stag}} = \begin{pmatrix} 0 & 1 \\ -c^2 & 0 \end{pmatrix}, \bar{D}_{\text{stag}} = \begin{pmatrix} 0 & c \\ -c & 0 \end{pmatrix}, \bar{D}_{\text{stag}} + t \bar{D}_{\text{stag}} = 0.
$$

- Diffusion numérique nulle ✓

- Fort couplage entre cellules ✓
1 Déterminer une matrice $D$ vérifiant les critères de stabilité ET de précision

- $\bar{D} + t\bar{D} \geq 0$ (stabilité)

- $\bar{D} + t\bar{D}$ ne contient pas la vitesse du son $c$ (precision)

2 calculer les flux entre cellules à partir de l'expression :

$$F_{i,j} = \frac{F(U_i) + F(U_j)}{2} + D \frac{U_i - U_j}{2}$$

Remarques :

- Le schéma peut être colocalisé comme décalé, une step ou multistep

- Le schéma est en général implicite si $Sp(\bar{D}) \subset i\mathbb{R}$

- On peut avoir une dissipation nulle $||\bar{U}_h(\cdot, t)||_{L^2} = ||\bar{U}_{h0}||_{L^2}$ et un décentrement $\bar{D}$ non nul
1 Déterminer une matrice $D$ vérifiant les critères de stabilité ET de précision

- $\bar{D} + t\bar{D} \geq 0$ (stabilité)

- $\bar{D} + t\bar{D}$ ne contient pas la vitesse du son $c$ (precision)

2 calculer les flux entre cellules à partir de l'expression :

$$F_{i,j} = \frac{F(U_i) + F(U_j)}{2} + D \frac{U_i - U_j}{2}$$

Remarques :

- Le schéma peut être colocalisé comme décalé, une step ou multistep

- Le schéma est en général implicite si $Sp(\bar{D}) \subset i\mathbb{R}$

- On peut avoir une dissipation nulle $||\bar{U}_h(\cdot, t)||_{L^2} = ||\bar{U}_{h0}||_{L^2}$ et un décentrement $\bar{D}$ non nul
We consider the following Euler equations :

$$ \partial_t \rho + \nabla \cdot (\rho \vec{u}) = 0 \quad (14) $$

$$ \partial_t (\rho \vec{u}) + \nabla \cdot (\rho \vec{u} \otimes \vec{u}) + \vec{\nabla} p = 0 \quad (15) $$

$$ \partial_t (\rho(e + \frac{1}{2}||\vec{u}||^2)) + \nabla \cdot ((\rho e + \frac{1}{2}\rho||\vec{u}||^2 + p)\vec{u}) = 0 \quad (16) $$

Symmetrisation with $\frac{\partial p}{\partial \rho}|_e$, $\frac{\partial p}{\partial e}|_\rho$ not appearing on the diagonal :

$$ S = \begin{pmatrix} u & \sqrt{\left.\frac{\partial p}{\partial \rho}\right|_e} & 0 \\ \sqrt{\left.\frac{\partial p}{\partial \rho}\right|_e} & u & \sqrt{\left.\frac{\partial p}{\partial e}\right|_\rho \frac{p}{\rho^2}} \\ 0 & \sqrt{\left.\frac{\partial p}{\partial e}\right|_\rho \frac{p}{\rho^2}} & u \end{pmatrix} . $$
We consider the following Euler equations :

$$ \partial_t \rho + \nabla \cdot (\rho \vec{u}) = 0 \quad (14) $$

$$ \partial_t (\rho \vec{u}) + \nabla \cdot (\rho \vec{u} \otimes \vec{u}) + \vec{\nabla} p = 0 \quad (15) $$

$$ \partial_t (\rho(e + \frac{1}{2}||\vec{u}||^2)) + \nabla \cdot ((\rho e + \frac{1}{2}\rho||\vec{u}||^2 + p)\vec{u}) = 0 \quad (16) $$

Symmetrisation with $\frac{\partial p}{\partial \rho}|_e$, $\frac{\partial p}{\partial e}|_\rho$ not appearing on the diagonal :

$$ S = \begin{pmatrix} u & \sqrt{\left.\frac{\partial p}{\partial \rho}\right|_e} & 0 \\ \sqrt{\left.\frac{\partial p}{\partial \rho}\right|_e} & u & \sqrt{\left.\frac{\partial p}{\partial e}\right|_\rho \frac{p}{\rho^2}} \\ 0 & \sqrt{\left.\frac{\partial p}{\partial e}\right|_\rho \frac{p}{\rho^2}} & u \end{pmatrix} . $$
Dans la base symétrisante :

$$
\bar{D}_{\text{Stag}} = \begin{pmatrix} |u| & \sqrt{\left.\frac{\partial p}{\partial \rho}\right|_e} & 0 \\ -\sqrt{\left.\frac{\partial p}{\partial \rho}\right|_e} & |u| & -\sqrt{\left.\frac{\partial p}{\partial e}\right|_{\rho} \frac{p}{\rho^2}} \\ 0 & \sqrt{\left.\frac{\partial p}{\partial e}\right|_{\rho} \frac{p}{\rho^2}} & |u| \end{pmatrix}
$$

Dans la base physique :

$$
D_{\text{Stag}} = \begin{pmatrix} |u| - u & 1 & 0 \\ -\left.\frac{\partial p}{\partial \rho}\right|_e - \nu^2 + \frac{1}{\rho} \left.\frac{\partial p}{\partial e}\right|_{\rho} \left(\frac{\nu^2}{2} - e\right) & u + |u| - \frac{u}{\rho} \left.\frac{\partial p}{\partial e}\right|_{\rho} & \frac{1}{\rho} \left.\frac{\partial p}{\partial e}\right|_{\rho} \\ -u \left(c^2 + E - \frac{p}{\rho} + \frac{1}{\rho} \left(\frac{\nu^2}{2} - e\right) \frac{\partial p}{\partial e}\right)_{\rho} & E - \frac{p}{\rho} - \frac{\nu^2}{\rho} \left.\frac{\partial p}{\partial e}\right|_{\rho} & \frac{1}{\rho} \left(u \left.\frac{\partial p}{\partial e}\right|_{\rho} + |u| \rho\right) \end{pmatrix}
$$
Dans la base symétrisante :

$$ \bar{D}_{\text{Stag}} = \begin{pmatrix} |u| & \sqrt{\left.\frac{\partial p}{\partial \rho}\right|_e} & 0 \\ -\sqrt{\left.\frac{\partial p}{\partial \rho}\right|_e} & |u| & -\sqrt{\left.\frac{\partial p}{\partial e}\right|_{\rho} \frac{p}{\rho^2}} \\ 0 & \sqrt{\left.\frac{\partial p}{\partial e}\right|_{\rho} \frac{p}{\rho^2}} & |u| \end{pmatrix} $$

Dans la base physique :

$$ D_{\text{Stag}} = \begin{pmatrix} |u| - u & 1 & 0 \\ -\left.\frac{\partial p}{\partial \rho}\right|_e - \nu^2 + \frac{1}{\rho} \left.\frac{\partial p}{\partial e}\right|_{\rho} \left(\frac{\nu^2}{2} - e\right) & u + |u| - \frac{u}{\rho} \left.\frac{\partial p}{\partial e}\right|_{\rho} & \frac{1}{\rho} \left.\frac{\partial p}{\partial e}\right|_{\rho} \\ -u \left(c^2 + E - \frac{p}{\rho} + \frac{1}{\rho} \left(\frac{\nu^2}{2} - e\right) \left.\frac{\partial p}{\partial e}\right|_{\rho}\right) & E - \frac{p}{\rho} - \frac{\nu^2}{\rho} \left.\frac{\partial p}{\partial e}\right|_{\rho} & \frac{1}{\rho} \left(u \left.\frac{\partial p}{\partial e}\right|_{\rho} + |u| \rho\right) \end{pmatrix} $$
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

Une nouvelle classe de schémas décalés

**Résultats numériques**

Conclusion et perspectives
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

Une nouvelle classe de schémas décalés

Résultats numériques

Conclusion et perspectives
## Problèmes de Riemann difficiles

* Evaluation de la Tests robustesse de méthodes numériques

* Solutions exactes disponibles → évaluation de la précision

* Evaluation des méthodes numériques en présence d'ondes de choc et de contact

* La loi d'état est celle d'un gaz parfait

<table><thead><tr><th>Test</th><th>ρ<sub>L</sub></th><th>u<sub>L</sub></th><th>ρ<sub>L</sub></th><th>ρ<sub>R</sub></th><th>u<sub>R</sub></th><th>ρ<sub>R</sub></th></tr></thead><tbody><tr><td>Toro1</td><td>1</td><td>0.75</td><td>1</td><td>0.125</td><td>0</td><td>0.1</td></tr><tr><td>Toro2</td><td>1</td><td>-2</td><td>0.4</td><td>1</td><td>2</td><td>0.4</td></tr><tr><td>Toro3</td><td>1</td><td>0</td><td>1000</td><td>1</td><td>0</td><td>0.01</td></tr><tr><td>Toro4</td><td>5.99924</td><td>19.5975</td><td>460.894</td><td>5.99242</td><td>-6.19633</td><td>46.0950</td></tr><tr><td>Toro5</td><td>1</td><td>-19.5975</td><td>1000</td><td>1</td><td>-19.59745</td><td>0.01</td></tr><tr><td>Toro6</td><td>1.4</td><td>0</td><td>1</td><td>1.4</td><td>0</td><td>1</td></tr><tr><td>Toro7</td><td>1.4</td><td>0.1</td><td>1</td><td>1</td><td>0.1</td><td>1</td></tr></tbody></table>

Table: Initial data for the seven tests Toro1, ..., Toro7
## Problèmes de Riemann difficiles

* Evaluation de la Tests robustesse de méthodes numériques

* Solutions exactes disponibles → évaluation de la précision

* Evaluation des méthodes numériques en présence d'ondes de choc et de contact

* La loi d'état est celle d'un gaz parfait

<table><thead><tr><th>Test</th><th>ρ<sub>L</sub></th><th>u<sub>L</sub></th><th>ρ<sub>L</sub></th><th>ρ<sub>R</sub></th><th>u<sub>R</sub></th><th>ρ<sub>R</sub></th></tr></thead><tbody><tr><td>Toro1</td><td>1</td><td>0.75</td><td>1</td><td>0.125</td><td>0</td><td>0.1</td></tr><tr><td>Toro2</td><td>1</td><td>-2</td><td>0.4</td><td>1</td><td>2</td><td>0.4</td></tr><tr><td>Toro3</td><td>1</td><td>0</td><td>1000</td><td>1</td><td>0</td><td>0.01</td></tr><tr><td>Toro4</td><td>5.99924</td><td>19.5975</td><td>460.894</td><td>5.99242</td><td>-6.19633</td><td>46.0950</td></tr><tr><td>Toro5</td><td>1</td><td>-19.5975</td><td>1000</td><td>1</td><td>-19.59745</td><td>0.01</td></tr><tr><td>Toro6</td><td>1.4</td><td>0</td><td>1</td><td>1.4</td><td>0</td><td>1</td></tr><tr><td>Toro7</td><td>1.4</td><td>0.1</td><td>1</td><td>1</td><td>0.1</td><td>1</td></tr></tbody></table>

Table: Initial data for the seven tests Toro1, ..., Toro7
* Convergence de 6 tests sur 7

* Echec du test Toro2 (densité négative car régime proche du vide)

<figure><img src="image_10.png" /><figcaption>Figure: Densité du test Toro4</figcaption></figure>

Nécessité d'un benchmark similaire diphasique
Benchmark 2D et/ou avec termes source pour vérification des
codes de calcul
* Convergence de 6 tests sur 7

* Echec du test Toro2 (densité négative car régime proche du vide)

<figure><img src="image_10.png" /><figcaption>Figure: Densité du test Toro4</figcaption></figure>

Nécessité d'un benchmark similaire diphasique
Benchmark 2D et/ou avec termes source pour vérification des
codes de calcul
* Canal 1D

* chauffage

* changement de phase

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

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

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

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

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

<figure><img src="image_16.png" /></figure>
* Canal 1D

* chauffage

* changement de phase

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

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

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

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

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

<figure><img src="image_16.png" /></figure>
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

Une nouvelle classe de schémas décalés

Résultats numériques

Conclusion et perspectives
Contexte

Problématique bas Mach

Diffusion numérique de schémas classiques

Une nouvelle classe de schémas décalés

Résultats numériques

**Conclusion et perspectives**
Avec le nouveau schéma décalé (grille cartésienne)

* Simulation d'ondes de choc

* Simulation du chauffage

* Simulation du changement de phase

* Simulation d'une perte de charge singulière

Convergence vers la solution analytique (sauf proximité du vide).

Nous pouvons à présent étudier

* Des modèles diphasiques plus complexes (e.g. avec dérive)

* Analyse précise du bilan d'entropie

* Analyse précise du bilan du caractère TVD

* Extension aux maillages quelconques
Avec le nouveau schéma décalé (grille cartésienne)

* Simulation d'ondes de choc

* Simulation du chauffage

* Simulation du changement de phase

* Simulation d'une perte de charge singulière

Convergence vers la solution analytique (sauf proximité du vide).

Nous pouvons à présent étudier

* Des modèles diphasiques plus complexes (e.g. avec dérive)

* Analyse précise du bilan d'entropie

* Analyse précise du bilan du caractère TVD

* Extension aux maillages quelconques
Thanks! Any questions?
<figure><img src="image_17.png" /></figure>

Thanks! Any questions?