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

HAL
open science

# On a discontinuity capturing Finite Volume method that can solve fully incompressible flows

Esteban Coiffier, Michael Ndjinga

► To cite this version:

Esteban Coiffier, Michael Ndjinga. On a discontinuity capturing Finite Volume method that can solve fully incompressible flows. SMAI 2023 - Congrès de la Société des Mathématiques Appliquées et Industrielles, May 2023, Le Gosier, Guadeloupe, France. hal-04397820

HAL Id: hal-04397820

https://inria.hal.science/hal-04397820v1

Submitted on 16 Jan 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 1. Context and Aim

**Nuclear thermal hydraulics** : simulation for safety in nuclear reactors is based on multiphase flows models (liquid, vapor), **Low-Mach flows**, quasi-incompressible or totally **incompressible** flows (historically approximated by staggered schemes [2]). However, the presence of two phases can generate **discontinuities** (treated classicaly with Riemann Solvers [5] [1]).

→ In need of a scheme that can solve incompressible flows and capturing discontinuities.

Multiphase models are very complicated, we started from Euler barotropic equations :

$$ \left\{ \begin{array}{l} \partial_t \rho + \nabla \cdot \mathbf{q} = 0 \\ \partial_t \mathbf{q} + \nabla \cdot \left( \frac{\mathbf{q} \otimes \mathbf{q}}{\rho} \right) + \nabla p = \mathbf{f} \\ p = g(\rho) \end{array} \right. \qquad (1) $$

with $\rho$ density of the fluid, $q$ its momentum, $p$ the pressure and $f$ an external force

# 3. Numerical results

## Riemann Problem : computation of discontinuities

Numerical tests description : 1D problem, initial condition is composed of two constant state sepa-

PStaggered scheme for isothermal Euler equations

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

rated by a discontinuity, fluid is water.

## Lid Driven Cavity : Low Mach 2D

Numerical tests description : $u = 1m/s$ on the top, and equal to 0 on the other sides on a 50 x 50 cells cartesian mesh. Fluid is water again.

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

We converged to a stationary state with **viscosity** that was also implemented in the manner of the previous paragraph, we can observe the formation of vortices, which is not guaranteed by classical schemes at Low Mach numbers

## Stationary with gravity : Incompressible 1D

Numerical tests description : incompressible flow, and stationary test case, with the external force $f = \rho g$ being the gravity.

Implicit staggered scheme for isothermal Euler equations
1D channel under gravity

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

# 2. Process

The two main ingredients are :

## Numerical Dissipation of staggered schemes

We recall where the expression "numerical dissipation" comes from. For a two point semi-discretization

$$ \partial_t U + \partial_x \mathcal{F}(U) = 0 \quad \text{becomes} \quad U_i'(t) + \frac{F_{i,i+1} - F_{i-1,i}}{\Delta x} = 0 \qquad (2) $$

where the discrete flux is defined as

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

so that, when $\Delta x \to 0$ : $\partial_t U + \partial_x \mathcal{F}(U) = \Delta x \mathcal{D} \partial_{xx} U + \mathcal{O}((\Delta x)^2)$

$\mathcal{D}$ is the **numerical dissipation operator**.

## Example on wave system : Godunov scheme

$$ \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} \rho_i + \frac{(q_{i+1} - q_{i-1})}{2\Delta x} = \frac{c(\rho_{i+1} - 2\rho_i + \rho_{i-1})}{2\Delta x} \\ \\ \displaystyle \frac{d}{dt} q_i + \frac{c^2(\rho_{i+1} - \rho_{i-1})}{2\Delta x} = \frac{c(q_{i+1} - 2q_i + q_{i-1})}{2\Delta x}, \quad U_i := (\rho_i, q_i)^t \text{ and } \mathcal{D} := \begin{pmatrix} c & 0 \\ 0 & c \end{pmatrix} \end{array} \right. \qquad (3) $$

## Example on wave equation : Staggered scheme

$$ \left\{ \begin{array}{l} \Delta x \frac{d}{dt} \rho_i + (q_{i+1/2} - q_{i-1/2}) = 0 \\ \Delta x \frac{d}{dt} q_i + c^2(\rho_{i+1} - \rho_i) = 0 \end{array} \right. \quad U_i := (\rho_i, q_{i+1/2})^t \text{ and } \mathcal{D} := \begin{pmatrix} 0 & 1 \\ -c^2 & 0 \end{pmatrix} $$

## Primitive variables

In the incompressible regime the density is constant by definition. **It cannot be an unknown of the problem.** We need to replace the density by the pressure. Instead of solving (1) we solve the following system :

$$ \left\{ \begin{array}{l} \partial_t (\rho \mathbf{u}) + \nabla \cdot (\mathbf{u} \rho) = 0 \\ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\mathbf{u} \rho \otimes \mathbf{u}) + \nabla p = \mathbf{f} \end{array} \right. \qquad (4) $$

## Putting it together

Putting (4) into the form of (2) we get the following discretization,

$$ U(V_i)'(t) + \frac{G_{i,i+1} - G_{i-1,i}}{\Delta x} = 0, \quad U(V_i) := (\rho(p_i), q_i)^t, \quad V_i = (p_i, u_i)^t $$

where the discrete flux is defined as

$$ G_{i,j} = \frac{\mathcal{F}(U(V_i)) + \mathcal{F}(U(V_j))}{2} - \mathcal{D} \frac{V_i - V_j}{2}, \quad \mathcal{F}(U(V_i)) := \begin{pmatrix} \rho(p_i)\mathbf{u}_i \\ \rho(p_i)\mathbf{u}_i \otimes \mathbf{u}_i + p_i \mathbb{I}_d \end{pmatrix}, $$

with the numerical dissipation extracted from classical staggered schemes [4] (in multi-D) (we differ from [3] for example)

$$ \mathcal{D} = \text{sign}(\mathbf{u}_{ij} \cdot \mathbf{n}) \begin{pmatrix} \frac{\mathbf{u}_{ij} \cdot \mathbf{n}}{c_{ij}^2} & \rho(p)^t \mathbf{n} \\ -\mathbf{n} - \frac{(\mathbf{u}_{ij} \cdot \mathbf{n})}{c_{ij}^2} \mathbf{u}_{ij} & \rho(p)(\mathbf{u}_{ij} \otimes \mathbf{n} + (\mathbf{u}_{i,j} \cdot \mathbf{n})\mathbb{I}_d) \end{pmatrix} $$

with $\mathbf{u}_{ij} := \frac{\mathbf{u}_i \sqrt{\rho(p_i)} + \mathbf{u}_j \sqrt{\rho(p_j)}}{\sqrt{\rho(p_i)} + \sqrt{\rho(p_j)}}$ and $c_{ij} := \frac{p_i - p_j}{\rho(p_i) - \rho(p_j)}$ if $\rho(p_i) \neq \rho(p_j)$ and $p'(\rho(p_i))$ else.

Finally, we integrate the time implicitly we solve the non-linear system with a Newton scheme

# 4. Conclusions

* Computations of discontinuities are very good

* Low-mach regime is handled effectively

* Incompressible regime can be solved with the same scheme but the linear solver has trouble when the system is singular (Newmann and Wall only boundary conditions)

# 5. References

[1] Sergei K Godunov and I Bohachevsky. "Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics". In: *Matematičeskij sbornik* 47.3 (1959), pp. 271–306.

[2] Francis H Harlow and J Eddie Welch. "Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface". In: *The physics of fluids* 8.12 (1965), pp. 2182–2189.

[3] Raphaele Herbin, Walid Kheriji, and J-C Latché. "On some implicit and semi-implicit staggered schemes for the shallow water and Euler equations". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 48.6 (2014), pp. 1807–1857.

[4] Andrea Prosperetti and Grétar Tryggvason. *Computational methods for multiphase flow*. Cambridge university press, 2009.

[5] Philip L Roe. "Approximate Riemann solvers, parameter vectors, and difference schemes". In: *Journal of computational physics* 43.2 (1981), pp. 357–372.