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

HAL
open science

# A rough scheme to couple free and porous media

Jean-Marc Hérard

► To cite this version:

Jean-Marc Hérard. A rough scheme to couple free and porous media. International Journal on Finite Volumes, 2006, 3 (2), http://www.latp.univ-mrs.fr/IJFV/spip.php?article17. hal-01114212

HAL Id: hal-01114212

https://hal.science/hal-01114212

Submitted on 10 Feb 2015

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.

Public Domain
# A rough scheme to couple free and porous media

Jean-Marc Hérard†

† Électricité de France, Recherche et Développement,
Département Mécanique des Fluides, Energies et Environnement,
6 quai Watier, 78401 Chatou cedex, FRANCE

jean-marc.herard@edf.fr

## Abstract

This paper is devoted to the computation of flows between free and porous media separated by a thin interface. The basic strategy relies on some ideas developed earlier by J.M. Greenberg and A.Y. Leroux on their work on well balanced schemes. This approach requires introducing a set of partial differential equations at the interface, in order to account for the sudden change of medium. The main features of the interface PDE are investigated. We afterwards propose to compute approximations of solutions with help of an approximate Godunov scheme. A linear interface Riemann solver is introduced, which aims at enforcing the continuity of the two (steady wave-) Riemann invariants. Numerical computations involving shock waves or rarefaction waves are examined and the agreement with the entropy inequality is tracked. Effects of the mesh refinement and the impact of the smoothing of the thin interface are also addressed in the paper.

**Key words** : Porous media / Interfacial coupling / Godunov scheme

## 1 Introduction

Quite recently, considerable attention has been paid to the analysis and the simulation of conservation laws involving discontinuous functions (see [3, 4, 5, 25, 26, 27, 29, 30] among others for instance). The problems involved may for instance occur when focusing on flows in porous media, which contain their intrinsic discontinuities. Actually, the structure of the conservation laws allows a deep investigation, and the connection through the interface between media is rather clear.

Another kind of problems arises in an industrial framework when one aims at performing complex applications involving several codes. When coupling existing
codes through some fictitious interface, in the framework of flows arising in the primary coolant circuit of nuclear power plants, at least three distinct phenomena may appear at the interface. Of course, models on each side may be different. One may for instance use a two-fluid approach on the one hand, and an homogeneous approach on the other hand. Alternatively, one may use a model M1 on the left hand side of the interface, and its counterpart on the right hand side of the interface, but assuming instantaneous relaxation with respect to some variables (typically, the velocity or the pressure). The third case, which seems to be the simplest one, occurs when the governing equations in each code rely on the same modelling assumptions. However the flow may :

(i) come from a one dimensional pipe into a three dimensional component [22],

(ii) come from a free domain into a domain including obstacles.

In that latter case, a standard approach consists in computing the ratio of the volume of obstacles over the total volume, which becomes the 'local' porosity. One advantage is that this enables to get rid of a precise description of the geometry of obstacles (wall boundaries) which would obviously occur when using a 3D local approach. A drawback nonetheless arises when coupling this code with another one devoted to the computation of flows in a free medium.

The present contribution provides some simple investigation on the latter problem (ii). It is part of a whole which aims at providing methods and tools to cope with the interfacial coupling of systems occurring in two-phase flow modelling (see also [1, 2, 12, 14] in that framework).

We will assume that there are no local strong discontinuities of porosity inside the porous medium. We will assume that the basic set of equations which will govern the motion of the flow through both the free and porous medium is :

$$ \frac{\partial W}{\partial t} + \frac{\partial f(W)}{\partial x} + \frac{\partial G(W)}{\partial x} = 0 \qquad (1) $$

when $x < 0, t > 0$, and :

$$ \frac{\partial \epsilon W}{\partial t} + \frac{\partial \epsilon f(W)}{\partial x} + \epsilon \frac{\partial G(W)}{\partial x} = 0 \qquad (2) $$

on the right side of the interface ($x > 0, t > 0$). The given function $\epsilon(x)$ is assumed to be regular. The state variable is $W = (\rho, \rho U)$, and we introduce the fluxes $f(W) = UW$, and $G(W) = (0, P(\rho))$. The set of equations (2) is of course consistent in the following sense: if one assumes that $U(x, t = 0) = 0$ and $\rho(x, t = 0) = \rho_0$, then the flow remains steady, even if the porosity $\epsilon$ is not uniform. Obviously system (2) retrieves a conservative form if $\epsilon$ is uniform.

The latter interface system (2) will be used to model the whole problem (with a discontinuous porosity at $x = 0$), and some of its properties will be examined. We will afterwards introduce the basic cell scheme which is based on Greenberg, Leroux and Gosse ideas ([16, 17, 18, 19, 20, 24]). The porosity will thus be discretised piecewise constant on each cell; moreover, the Riemann problem at each cell interface
will account for $\epsilon$ discontinuity. An approximate interface Riemann solver will be introduced, which ensures the continuity of the two underlying steady-wave Riemann invariants. The interface solver is based on the classical VFRoe-ncv approach ([6]), which is a sequel of the rough Godunov scheme [10]. Some basic experiments involving shock waves or rarefaction waves travelling through the porous interface will be simulated with these schemes. Effects of the mesh refinement will be studied. Eventually, the effect of smoothing of $\epsilon$ around the interface will be discussed, together with its impact on the converged approximations. This paper thus aims at highlighting some of the problems which may occur when computing this kind of flows.

## 2 A global model to describe the whole flow

Throughout the paper, $W = t(\epsilon, \epsilon\rho, \epsilon\rho U)$ will design the so-called conservative variable. Recall that the porosity $\epsilon$ only varies wrt to $x$ and is given by the user. Actually, $\rho, U$ are the mean density and the mean velocity in the field, and $\epsilon$ is expected to lie in $[0, 1]$ (that is, wall boundaries cannot be imposed by setting $\epsilon = 0$). Nonetheless, the limit case where $\epsilon = 0$, will be discussed in the last section which is devoted to numerical results. The flow is assumed highly unsteady and governed by both convective terms. The free flow on the LHS of the interface is governed by:

$$ \frac{\partial \rho}{\partial t} + \frac{\partial \rho U}{\partial x} = 0 \qquad (3) $$

$$ \frac{\partial \rho U}{\partial t} + \frac{\partial \rho U^2 + P(\rho)}{\partial x} = \frac{\partial \mu \frac{\partial U}{\partial x}}{\partial x} \qquad (4) $$

and the flow in the porous medium agrees with:

$$ \frac{\partial \epsilon}{\partial t} = 0 \qquad (5) $$

$$ \frac{\partial \epsilon \rho}{\partial t} + \frac{\partial \epsilon \rho U}{\partial x} = 0 \qquad (6) $$

$$ \frac{\partial \epsilon \rho U}{\partial t} + \frac{\partial \epsilon \rho U^2}{\partial x} + \epsilon \frac{\partial P(\rho)}{\partial x} = \frac{\partial \epsilon \mu \frac{\partial U}{\partial x}}{\partial x} \qquad (7) $$

The entropy inequality associated with the latter model (5), (6), (7) is classical. An admissible entropy-entropy flux is: $(\eta, f_{\eta})$ where:

$$ \eta = \epsilon \rho (U^2/2 + \phi(\rho)) \qquad (8) $$

$$ f_{\eta} = \epsilon \rho U (U^2/2 + \phi(\rho) + P(\rho)/\rho) \qquad (9) $$

noting:

$$ \phi(\rho) = \int_{0}^{\rho} \frac{P(a)}{a^2} da \qquad (10) $$
For regular solutions of ((5)-(7)), the governing equation for the entropy reads:

$$ \frac{\partial \eta}{\partial t} + \frac{\partial f_{\eta}}{\partial x} - \frac{\partial \epsilon \mu u \frac{\partial U}{\partial x}}{\partial x} = - \epsilon \mu \left( \frac{\partial U}{\partial x} \right)^2 \quad (11) $$

In the limit of a vanishing viscosity, this provides :

$$ \frac{\partial \eta}{\partial t} + \frac{\partial f_{\eta}}{\partial x} \le 0 \quad (12) $$

We from now on will neglect viscous contributions everywhere. Thus, unlike in more complex cases involving different flows, there is no ambiguity here concerning the suitability of the latter model ((5)-(7)) to describe the motion of flows on each side of the interface, taking $\epsilon = 1$ for the free flow. For conveniency, we introduce the homogeneous flux $F(W) = (0, \epsilon\rho U, \epsilon\rho U^2)$. We recall that the function $P(\rho)$ is monotone increasing with: $P(0) = 0$. Numerical tests will be performed using the standard law : $P(\rho) = P_0(\rho/\rho_0)^\beta$.

From now on, we assume that $\mu = 0$ and that (12) holds.

One still needs to provide "interface conditions" to couple both systems.

## 3 Main properties of the non viscous system

### Property 1

The left hand side of system ((5)-(7)) is hyperbolic. Eigenvalues read: $\lambda_1 = 0$, $\lambda_2 = U - c$, $\lambda_3 = U + c$, setting $c = (P'(\rho))^{1/2}$. The field associated with $\lambda_1$ is Linearly Degenerated and the fields $\lambda_2$ and $\lambda_3$ are Genuinely Non Linear.

The proof is obvious. For regular solutions with a continuous porosity, a straightforward counterpart of ((5)-(7)) may be written in a 'weakly' conservative form as:

$$ \frac{\partial \epsilon}{\partial t} = 0 \quad (13) $$

$$ \frac{\partial \epsilon \rho}{\partial t} + \frac{\partial \epsilon \rho U}{\partial x} = 0 \quad (14) $$

$$ \frac{\partial U}{\partial t} + \frac{\partial \frac{U^2}{2}}{\partial x} + \frac{\partial \psi(\rho)}{\partial x} = 0 \quad (15) $$

while noting: $\psi(\rho) = \int_0^\rho \frac{c^2(a)}{a} da$.

We recall that even when $\epsilon = \epsilon_0$ is uniform, the jump conditions associated with ((13), (15)) are not equivalent to those derived from the LHS of ((5), (7)).
## 3.1 Regular waves and Riemann invariants

### Property 2

(i) - The Riemann invariants in the stationary wave are defined as

$$I_1^1 = \rho \epsilon U \qquad (16)$$

$$I_1^2 = \frac{U^2}{2} + \psi(\rho) \qquad (17)$$

(ii) - The Riemann invariants in the 2 wave (resp. the 3 wave) associated with $\lambda_2$ (resp. $\lambda_3$) are:

$$I_2^1 = \epsilon \qquad (18)$$

$$I_2^2 = U + \int_0^\rho \frac{c(a)}{a} da \qquad (19)$$

and:

$$I_3^1 = \epsilon \qquad (20)$$

$$I_3^2 = U - \int_0^\rho \frac{c(a)}{a} da \qquad (21)$$

We note that the entropy flux $f_\eta$ corresponds to the product of Riemann invariants in the 1-wave.

## 3.2 Interface conditions between media

We may now propose two natural candidates for the -steady- interface conditions between the free medium and the porous medium.

* A direct consequence of (16)-(17) is that an obvious candidate to couple both media is to set:

$$[I_1^1]_l^r = [\rho \epsilon U]_l^r = 0 \qquad (22)$$

$$[I_1^2]_l^r = [\frac{U^2}{2} + \psi(\rho)]_l^r = 0 \qquad (23)$$

through the sharp interface. This strong candidate does not account for any head loss.

* A second "weak" candidate immediately follows, which may account for the fact that the entropy dissipation should agree with (12), that is:

$$[I_1^1]_l^r = 0 \qquad (24)$$

$$[f_\eta]_l^r \le 0 \qquad (25)$$

This second one may account for the superposition of a genuinely non linear wave with the 1-contact wave. Both candidates incorporate the same condition corresponding to the continuity of the mass flow rate $I_1^1 = \rho \epsilon U$.
### 3.3 Shock waves and jump conditions

**Case 1 : $\epsilon$ is continuous**

Obviously, if some shock is located at some point where the porosity is regular, jump conditions simply write:

$$[\epsilon] = 0 \qquad (26)$$

$$-\sigma[\rho] + [\rho U] = 0 \qquad (27)$$

$$-\sigma[\rho U] + [\rho U^2 + P(\rho)] = 0 \qquad (28)$$

Thus the entropy inequality (12) provides some counterpart of $-\sigma[\eta] + [f_\eta] < 0$ which reads:

$$(P_l + P_r - 2\rho_l \rho_r [\phi]_l^r / [\rho]_l^r) [U]_l^r < 0 \qquad (29)$$

which, for standard EOS $P(\rho)$, turns out to : $[U]_l^r < 0$.

**Case 2 : $\epsilon$ is discontinuous (sharp interface)**

If some shock wave happens to locate precisely at $x_0$ where $\epsilon$ is discontinuous, we may formally set for some $\hat{\epsilon}$:

$$-\sigma[\epsilon\rho] + [\epsilon\rho U] = 0 \qquad (30)$$

$$-\sigma[\epsilon\rho U] + [\epsilon\rho U^2] + \hat{\epsilon}[P(\rho)] = 0 \qquad (31)$$

For a stationary shock wave ($\sigma = 0$), and for given sign of the velocity $u_l > 0$ on the left side of the discontinuity, the entropy inequality will yield:

$$(\epsilon\rho U)_l [u^2/2 + \phi(\rho) + P(\rho)/\rho]_l^r < 0 \qquad (32)$$

Hence, inserting $\sigma = 0$ and ((30)-(31)) in the latter ((32)) leads to:

$$\hat{E} = (-\hat{\epsilon}(\frac{1}{\epsilon_l \rho_l} + \frac{1}{\epsilon_r \rho_r})[P(\rho)]_l^r + 2[\psi(\rho)]_l^r) < 0 \qquad (33)$$

The whole will not depend on the values of $\epsilon_l$ and $\epsilon_r$ as soon as either:

$$\hat{\epsilon}(\frac{1}{\epsilon_l \rho_l} + \frac{1}{\epsilon_r \rho_r}) = \frac{1}{\rho_l} + \frac{1}{\rho_r} \qquad (34)$$

which simply corresponds to the average :

$$\hat{\epsilon} = \hat{\epsilon}^{\text{max}}(W_l, W_r) = \frac{\epsilon_l m_r + \epsilon_r m_l}{m_l + m_r} = \epsilon_l \epsilon_r \frac{\rho_r + \rho_l}{m_l + m_r} \qquad (35)$$

noting $m = \epsilon\rho$, or if :

$$\hat{\epsilon} = \hat{\epsilon}^0(W_l, W_r) = \frac{2m_l m_r [\psi(\rho)]_l^r}{m_l + m_r [P(\rho)]_l^r} = \epsilon_l \epsilon_r \frac{2\rho_r \rho_l [\psi(\rho)]_l^r}{m_l + m_r [P(\rho)]_l^r} \qquad (36)$$
when $\rho_l$ and $\rho_r$ are distinct, and : $\hat{\epsilon}^0(W_l, W_r) = \frac{2\epsilon_l\epsilon_r}{\epsilon_l+\epsilon_r}$ otherwise. If (35) is used, $\hat{E}$ reaches it lower bound, otherwise if (36) is used, $\hat{E} = 0$. Both provide positive formulas for the mean porosity at the interface $\hat{\epsilon}$. Meaningful jump solutions must also agree with:

$$m_l m_r ([U]_l^r)^2 = \hat{\epsilon}[P(\rho)]_l^r [m]_l^r > 0 \quad (37)$$

The formula (35) is the only one which ensures that: ($\epsilon_l = \epsilon_r$ implies: $\hat{\epsilon}(W_l, W_r) = \epsilon_r$).

Thus it is the only one which guarantees a continuous transition to natural jump relations (26) given above. Consequently, the second one (36) will be disregarded from now on in this study.

Whenever (35) or (36) is applied, a straightforward consequence is the following:

### Property 3

Using either (36) or (35), and assuming that no fluid is present on the right side only (say $\epsilon_r = 0$), then the flow exactly behaves as if some wall boundary condition were imposed.

The interface value of porosity (35) will be used in the approximate interface Riemann solver discussed in the section below.

# 4 An approximate Godunov scheme

We now discuss some possible way to discretize the LHS of ((5), (7)).

## 4.1 Finite Volume procedure

The explicit Finite Volume procedure, which takes its grounds on Greenberg and Leroux's work (see [16, 17, 18, 19, 20, 24], ) is defined as follows. It aims at accounting for discontinuous values of porosity, which cannot be achieved when using classical methods. For a given mesh size : $h_i = x_{i+1/2} - x_{i-1/2}$, and a time step $\Delta t$ in agreement with the classical CFL condition, we compute $W_i^{n+1}$ for given data $W_i^n$ as:

$$h_i((\epsilon\rho)_i^{n+1} - (\epsilon\rho)_i^n) + \Delta t(Q_{i+1/2}^{-} - Q_{i-1/2}^{+}) = 0 \quad (38)$$

$$h_i((Q)_i^{n+1} - (Q)_i^n) + \Delta t \left( \frac{(Q_{i+1/2}^{-})^2}{(\rho\epsilon)_{i+1/2}^{-}} - \frac{(Q_{i-1/2}^{+})^2}{(\rho\epsilon)_{i-1/2}^{+}} \right) = 0 \quad (39)$$

$$+ \Delta t \frac{(\epsilon_{i+1/2}^{-} + \epsilon_{i-1/2}^{+})}{2} (P(\rho_{i+1/2}^{-}) - P(\rho_{i-1/2}^{+})) = 0 \quad (40)$$
where $Q = \epsilon\rho U$. We moreover introduce piecewise constant approximation of the porosity within each cell. Thus we get $\epsilon_{i+1/2}^- = \epsilon_i$ and similarly : $\epsilon_{i+1/2}^+ = \epsilon_{i+1}$.

It only remains to compute the states with subscript +, - on each side of the steady interface $x/t = 0$.

## 4.2 An interface solver preserving 0-Riemann invariants

We consider herein the approximate Godunov scheme which has been introduced in [6], and is a sequel of the approximate Godunov scheme [10]. The only thing we need to introduce is the non-conservative variable $Y$. The main idea here is to preserve in the linearised approximate Riemann solver the properties of the exact Riemann solver through linearly degenerate fields, unlike in the classical Roe's approach which privileges the capture of shocks, and thus introduces the notion of consistency with respect to the integral form of the conservation laws. Hence, one attempts here to maintain Riemann invariants of the true 0-wave through the approximate discrete 0-contact wave. This may be achieved with help of the variable $Y = (\epsilon, Q, \phi)$, noting conveniently $\phi = I_1^2$. This has been preferred to the former choice $Z = (\epsilon, 2c, U)$ that was made in [9], which was devoted to the accurate approximation of strong rarefaction waves including vacuum areas. Actually, one can easily check that this choice $Z = (\epsilon, 2c, U)$ inhibits to recover the exact preservation of 0-Riemann invariants in the linearized interface solver.

The linearisation of the LHS of ((5), (7)) will thus be based on :

$$ \frac{\partial Y}{\partial t} + B(Y) \frac{\partial Y}{\partial x} = 0 \qquad (41) $$

where:

$$ B(Y) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & u & \epsilon\rho \\ 0 & \frac{\psi'(\rho)}{\epsilon} & u \end{pmatrix} $$

System (41) may be symmetrized, when multiplied on the left by the following matrix:

$$ S(Y) = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & (\frac{\epsilon\rho}{c})^2 \end{pmatrix} $$

The right eigenvectors of the matrix $B(Y)$ associated with the eigenvalues $\lambda_1, \lambda_2, \lambda_3$ are:

$$ r_1(Y) = (1, 0, 0) \qquad (42) $$

$$ r_2(Y) = (0, \epsilon\rho, -c) \qquad (43) $$

$$ r_3(Y) = (0, \epsilon\rho, c) \qquad (44) $$
We note that $r_1, r_2, r_3$ always span $IR^3$ (unless $\epsilon = 0$ or if the vacuum occurs).

If $\Omega(Y)$ stands for the matrix of right eigenvectors of $B(Y)$, then:

$$ (\Omega(Y))^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{1}{2\epsilon\rho} & -\frac{1}{2c} \\ 0 & \frac{1}{2\epsilon\rho} & \frac{1}{2c} \end{pmatrix} $$

The approximate Godunov solver [6] requires to compute the exact solution of the linearised problem:

$$ \frac{\partial Y}{\partial t} + B(\hat{Y}(Y_L, Y_R)) \frac{\partial Y}{\partial x} = 0 \qquad (45) $$

while setting here:

$$ \hat{Y}(Y_L, Y_R) = (Y_L + Y_R)/2 $$

If $\tilde{\lambda}_k$ (respectively $\tilde{r}_k$) stands for the k-th eigenvalue (respectively the k-th right eigenvector) of $B(\hat{Y}(Y_L, Y_R))$ and if we set:

$$ Y^+ = Y(x/t = 0+) \qquad (46) $$

$$ Y^- = Y(x/t = 0-) \qquad (47) $$

we obtain the final form of intermediate states:

If $\tilde{\lambda}_3 < 0$, then:

$$ Y^+ = Y_R \qquad (48) $$

$$ Y^- = Y_R - \alpha_1 \tilde{r}_1 \qquad (49) $$

If $\tilde{\lambda}_2 < 0$ and $\tilde{\lambda}_3 > 0$, then:

$$ Y^+ = Y_R - \alpha_3 \tilde{r}_3 \qquad (50) $$

$$ Y^- = Y_L + \alpha_2 \tilde{r}_2 \qquad (51) $$

If $\tilde{\lambda}_2 > 0$, then:

$$ Y^+ = Y_L + \alpha_1 \tilde{r}_1 \qquad (52) $$

$$ Y^- = Y_L \qquad (53) $$

This must be complemented with:

If $\tilde{\lambda}_2 = 0$, then:

$$ Y^+ = Y_R - \alpha_3 \tilde{r}_3 \qquad (54) $$

$$ Y^- = Y_L \qquad (55) $$

If $\tilde{\lambda}_3 = 0$, then:

$$ Y^+ = Y_R \qquad (56) $$

$$ Y^- = Y_L + \alpha_2 \tilde{r}_2 \qquad (57) $$
We recall that coefficients $\alpha_k$ are computed as:

$$ (\alpha_1, \alpha_2, \alpha_3)^t = (\Omega(\hat{Y}))^{-1}(Y_R - Y_L) \quad (58) $$

We have the following property, the proof of which is obvious, owing to the form of the first right-eigenvector $\tilde{r}_1$:

### Property 4

This interface solver is in agreement with (16, 17) since it ensures that the numerical intermediate states computed on each side of the contact discontinuity associated with $\lambda_1$ comply with:

$$ Q_{i+1/2}^{-} = Q_{i+1/2}^{+} \quad (59) $$

$$ (I_1^2)_{i+1/2}^{-} = (I_1^2)_{i+1/2}^{+} \quad (60) $$

(recall that $I_1^2 = \frac{U^2}{2} + \psi(\rho)$).

The whole scheme (38), (39) thus ensures the conservation of mass. It also locally ensures the conservation of momentum provided that $\epsilon_i = \epsilon_{i+1}$.

Of course, we have in general: $(\rho)_{i+1/2}^{-} \neq (\rho)_{i+1/2}^{+}$, and also: $u_{i+1/2}^{-} \neq u_{i+1/2}^{+}$, unless: $\epsilon_i = \epsilon_{i+1}$.

# 5 Numerical results

The cell scheme is (38), (39), and we use the interface Riemann solver detailed above. We also focus in a subsection on the effect of "smoothing" of the interface (see figures (1-8)).

The state law is:

$$ P(\rho) = P_0 \rho^\beta, \text{ with } \beta = 3, \text{ and } P_0 = 10^5. $$

Owing to (3) and (6), we expect that through the stationary interface:

$$ Q_l = Q_r \quad (61) $$

if $l, r$ subscripts refer to states on both sides of $x/t = 0$.

Owing to (12) we also expect that either $[I_1^2]_l^r = 0$, or:

$$ [Q I_1^2]_l^r \le 0 \quad (62) $$

which turns to be: $[I_1^2]_l^r \le 0$ when the flow is right going ($u_l > 0$), and $0 \le [I_1^2]_l^r$ when the flow is going to the left ($u_l < 0$).

The suitability of the numerical answer may thus be examined on the basis of this very simple entropy principle.
## 5.1 First test case

We start with the following initial data. The porosity agrees with:

$$ \epsilon(x < 0.5) = 1. \qquad (63) $$

$$ \epsilon(x > 0.5) = 0.5 \qquad (64) $$

and the initial conditions are:

$$ \rho(x < 0.4, t = 0) = 2., U(x < 0.4, t = 0) = 0. \qquad (65) $$

$$ \rho(x > 0.4, t = 0) = 1., U(x > 0.4, t = 0) = 0. \qquad (66) $$

The mesh contains 4000 nodes, and the CFL number has been set to 1. A shock wave (respectively a rarefaction wave) starting from $x = 0.4$ travels to the right (respectively to the left), and hits the porosity discontinuity located at $x = 0.5$. Part of the shock wave is reflected, and part is transmitted through the right. We may rely on the approximations on each side of $x = 0.5$, and we thus focus on the behaviour of $\epsilon\rho U$ which we expect to be continuous (at least on sufficiently fine meshes), and of $I_1^2$. When the flow around the porous interface is steady, the computed approximations of the momentum around $x = 0.5$ are almost uniform, and the discrete form of the entropy inequality is preserved.

## 5.2 Second test case

The porosity distribution is exactly the same. The initial condition for the density and the velocity are:

$$ \rho(x < 0.4, t = 0) = 1., U(x < 0.4, t = 0) = 0. \qquad (67) $$

$$ \rho(x > 0.4, t = 0) = 2., U(x > 0.4, t = 0) = 0. \qquad (68) $$

The mesh and the CFL number are the same as in the previous case. A shock wave (respectively a rarefaction wave) starting from $x = 0.4$ travels to the left (respectively to the right). The head of the 2 rarefaction wave hits the porosity discontinuity in $x = 0.5$. $\epsilon\rho U$ remains continuous through the latter interface, and the jump of $I_1^2$ is in agreement with the entropy inequality.

## 5.3 Third test case

The porosity profile is now:

$$ \epsilon(x < 0.5) = 1. \qquad (69) $$

$$ \epsilon(x > 0.5) = 0.05 \qquad (70) $$

The density and velocity initial condition are exactly those of the first test case, and the mesh and the CFL number are also the same as in the first case. The amplitude of the incoming right-going shock is thus exactly the same, but the contrast in porosity is much higher (and probably meaningless when focusing on flows in nuclear power plants). We may check that the behaviour of steady cell values of momentum $\epsilon\rho U$ are continuous around $x = 0.5$, and that the second Riemann invariant $I_1^2$ correctly varies, since we expect that $[I_1^2] < 0$ since $u_l(x = 0.5, T) > 0$. This correct behaviour remains when the mesh is refined.
## 5.4 Fourth test case

The porosity profile is still:

$$ \epsilon(x < 0.5) = 1. \qquad (71) $$

$$ \epsilon(x > 0.5) = 0.05 \qquad (72) $$

The density and velocity initial condition are exactly those of the second test case, and the mesh and the CFL number are also the same as in the second test case. We still notice here that $[Q] = 0$ and $[I_1^2] > 0$ through the porous interface, which is in agreement with the entropy condition since $u_l(x = 0.5, T) < 0$.

## 5.5 Effect of sharp or smooth representation of the free-porous interface

The porosity profile is again:

$$ \epsilon(x < 0.5) = 1. \qquad (73) $$

$$ \epsilon(x > 0.5) = 0.05 \qquad (74) $$

and the initial condition are :

$$ \rho(x < 0.4, t = 0) = 2., U(x < 0.4, t = 0) = 0. \qquad (75) $$

$$ \rho(x > 0.4, t = 0) = 1., U(x > 0.4, t = 0) = 0. \qquad (76) $$

We wish to examine here the effect of smoothing of the interface. This is achieved inserting 10 points in the $\epsilon$ profile through the interface. The profile may either be linear or parabolic. The plot also includes the Heavyside distribution (sharp interface) (see figures (1) to (4)).

### Remark 1

We emphasize that we still use the same approximate Godunov scheme with variable $Y = (\epsilon, Q, \phi)$ for all interfaces, whatever the regularization of the porous profile is.

### Remark 2

One might consider the counterpart of this numerical approach from a continuous point of view, by regularizing $\epsilon$ on a shrinking domain, and then pass to the limit as has been achieved for generalized shallow-water equations in [7] (see [19] also). It seems doubtful anyway that such a regularization of $\epsilon$ will enable to retrieve the relation (35) (see [28]).

A similar test is performed setting:

$$ \epsilon(x > 0.5) = 0.05 \qquad (77) $$

$$ \epsilon(x < 0.5) = 1. \qquad (78) $$

and the initial conditions are :

$$ \rho(x < 0.4, t = 0) = 1., U(x < 0.4, t = 0) = 0. \qquad (79) $$

$$ \rho(x > 0.4, t = 0) = 2., U(x > 0.4, t = 0) = 0. \qquad (80) $$
A series of four figures (see figures (5) to (8)) provides the numerical approximations obtained using a fine mesh with 3800 cells.

## 5.6 Comparison with a pure reflexion on a wall boundary

One aims here at comparing approximated solutions achieved with the interface solver, and those obtained when computing the pure reflexion on a wall of the same incoming shock wave. The mesh is rather fine here (3800 cells). Results are displayed on figures (9), (10), (11), (12). The wall boundary -in the case $\epsilon_R = 0.$- is located at $x = 0.5$ (respectively the coupling interface -in the case $\epsilon_R = 0.05-$).

# 6 Conclusion

This work aims at highlighting some of the difficulties occuring when computing any system of non linear PDE describing fluid mechanics when taking porosity into account. While focusing on a simple two-equation model, it appears that steady discontinuous patterns may appear, even when computing simple regular flows such as rarefaction waves traveling in pipes (this actually exactly stands for the Loss Of Coolant Accident in nuclear power plants), as soon as these pipes contain some sudden enlargement or contraction.

Two physical constraints should be fulfilled through the interface, that is:

$$ (\epsilon\rho U)_l = (\epsilon\rho U)_r \qquad (81) $$

and :

$$ \epsilon_l \rho_l U_l \left[ \frac{U^2}{2} + \int_0^\rho \frac{c^2(a)}{a} da \right]_l^r \le 0 \qquad (82) $$

A cell scheme which relies on Leroux ideas has been examined. An *interface* Riemann solver has been proposed, which enforces the continuity of both Riemann invariants of the steady wave. On the basis of many numerical experiments, it appears that the latter is suitable, at least when the contrast of porosity through the interface is not too large ( for nuclear applications, this assumption is meaningful). Moreover, it occurs that there is no influence of the interface path when using this approach. We wish to emphasize that another *interface* solver based on the same approach, but using the state variable $Y' = (\epsilon, Q, \rho\epsilon)$ instead of $Y = (\epsilon, Q, \phi)$ which was considered herein, may violate the entropy inequality (12) in some cases, even when the contrast $\epsilon_r/\epsilon_l$ between media is rather close to 1. (see [21] for more details.

One might wonder whether a higher order scheme would be suitable to cope with the interface. Actually, our opinion is that this is not desirable at the coupling interface for stability reasons, whereas it is indeed a good and standard ingredient apart from the coupling interface, which enables to increase the local accuracy. In order to improve the method, it now remains to examine whether one may prescribe
a physical head loss provided by the user.

Acknowledgments

This work has received financial support from the NEPTUNE project, which is spon-
sored by EDF (Electricité de France), CEA (Commissariat à l'Energie Atomique),
IRSN (Institut de Radioprotection et de Sûreté Nucléaire) and AREVA-NP.

References

[1] A. AMBROSO, C. CHALONS, F. COQUEL, E. GODLEWSKI, F. LAGOUTIERE, P.A. RAVIART, N. SEGUIN, Homogeneous two-phase flows: coupling by Finite Volume methods, *Proceedings of FVCA IV*, Marrakech, July 4-8, 2005.

[2] A. AMBROSO, C. CHALONS, F. COQUEL, E. GODLEWSKI, F. LAGOUTIERE, P.A. RAVIART, N. SEGUIN, Coupling of multiphase flows, *Proceedings of NURETH11*, Avignon, October 2-6, 2005.

[3] F. BACHMANN, Analysis of a scalar conservation law with a flux function with discontinuous coefficients, *Advances in Differential Equations*, vol. 11-12, pp.1317-1338, 2004.

[4] F. BACHMANN AND J. VOVELLE, Existence and uniqueness of entropy solution of a scalar conservation law with a flux function involving discontinuous coefficients, *Communications in Partial Differential Equations*, to appear in 2006.

[5] F. BACHMANN, Finite Volume schemes for a non linear hyperbolic conservation law with a flux function involving discontinuous coefficients, *Int. J. of Finite Volumes*, vol. 3 (1), pdf file available on http://averoes.math.univ-paris13.fr/IJFV/, 2006.

[6] T. BUFFARD, T. GALLOUËT AND J.M. HÉRARD, A sequel to a rough Godunov scheme. Application to real gas flows, *Computers and Fluids*, 2000, vol. 29-7, pp. 813–847.

[7] A. CHINNAYYA, A. LEROUX AND N.S. SEGUIN, A well balanced numerical scheme for the approximation of the shallow-water equations with topography: the resonance phenomenon, *Int. J. of Finite Volumes*, vol. 1(1), pdf file available on http://averoes.math.univ-paris13.fr/IJFV/, 2004.

[8] B. EINFELDT, C.D. MUNZ, P.L. ROE AND B. SJÖGREEN, On Godunov-type methods near low densities, *J. Comp. Phys.*, 1991, vol. 92-2, pp. 273–295.

[9] T. GALLOUËT, J.M. HÉRARD AND N.S. SEGUIN, Some approximate Godunov schemes to compute shallow-water equations with topography, *Computers and Fluids*, 2003, vol. 32, pp.479–513.
[10] T. GALLOUËT, AND J.M. MASELLA, A rough Godunov scheme, *Comptes Rendus Académie des Sciences de Paris*, 1996, vol. I-323, pp. 77-84.

[11] P. GOATIN AND P.G. LE FLOCH, The Riemann problem for a class of resonant hyperbolic systems of balance laws, *Ann. Inst. H. Poincaré, Non linear Analysis*, vol. 21, pp. 881-902, 2004.

[12] E. GODLEWSKI, K.C. LE THANH AND P.A. RAVIART, The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. The case of systems, *Math. Mod. Num. Anal.*, vol. 39(4), pp.649-692 (2005).

[13] E. GODLEWSKI AND P.A. RAVIART, *Numerical approximation of hyperbolic systems of conservation laws*, Springer Verlag, 1996.

[14] E. GODLEWSKI AND P.A. RAVIART, The numerical interface coupling of non-linear hyperbolic systems of conservation laws. The scalar case, *Numerische Mathematik*, vol. 97, pp. 81-130 (2004).

[15] S.K. GODUNOV, A difference method for numerical calculation of discontinuous equations of hydrodynamics, *Mat. Sb.*, 1959, pp. 271-300. In Russian.

[16] L. GOSSE, A well balanced flux splitting scheme designed for hyperbolic systems of conservation laws with source terms, *Computers and Mathematics with Applications*, 2000, vol. 39, pp. 135-159.

[17] L. GOSSE, A well balanced flux splitting scheme using non conservative products designed for hyperbolic systems of conservation laws with source terms, *Math. Mod and Meth. in Appl. Sciences*, 2001, vol. 70, pp. 339-365.

[18] L. GOSSE AND A.Y. LEROUX, Un schéma équilibre adapté aux lois de conservation scalaires non homogènes, *Comptes Rendus Académie des Sciences de Paris*, 1996, vol. I-323, pp. 543-546.

[19] J.M. GREENBERG AND A.Y. LEROUX, A well balanced scheme for the numerical processing of source terms in hyperbolic equations, *SIAM J. of Num. Anal.*, 1996, vol. 33, pp. 1-16.

[20] J.M. GREENBERG, A.Y. LEROUX, R. BARAILLE AND A. NOUSSAIR, Analysis and approximation of conservation laws with source terms, *SIAM J. of Num. Anal.*, 1997, vol. 34, pp. 1980-2007.

[21] J.M. HERARD, *Naive schemes to couple isentropic flows between free and porous medium*, internal EDF report HI-81/04/008/A, 2004 (unpublished).

[22] J.M. HERARD AND O. HURISSE, *Coupling two and one dimensional model through a thin interface*, AIAA paper 2005-4718, 2005.

[23] E. ISAACSON AND B. TEMPLE, Convergence of the Godunov 2x2 method for a general resonant nonlinear balance law, *SIAM Journal of Applied Math.*, vol. 55, pp. 625-640, 1995.
[24] A.Y. LEROUX, *Discrétisation de termes sources raides dans les problèmes hyperboliques*, *Cours CEA-EDF-INRIA*, 1998, INRIA Rocquencourt, France.

[25] C. KLINGENBERG AND N.H. RISEBRO, Convex conservation laws with discontinuous coefficients, *Communications in Partial Differential Equations*, vol. 20, pp. 1959-1990, 1995.

[26] K.H. KARLSEN AND N.H. RISEBRO, Convergence of finite difference schemes for viscous and inviscid conservation laws with rough coefficients, *Math. Model. and Numer. Anal.*, vol. 35, pp. 239-269, 2001.

[27] K.H. KARLSEN AND N.H. RISEBRO AND J.D. TOWERS, Upwind difference approximations for degenerate parabolic convection-diffusion equations with a discontinuous coefficient, *IMA Journal Numer. Anal.*, vol. 22, pp. 623-664, 2002.

[28] N. SEGUIN, *Private communication*, 02/02/2006.

[29] N. SEGUIN AND J. VOVELLE, Analysis and approximation of a scalar conservation law with a flux function with discontinuous coefficients, *Math. Mod. and Meth. in Applied Sci.*, vol. 13, pp. 221-250, 2003.

[30] A. VASSEUR, Well posedness of scalar conservation laws with singular sources, *Methods Appl. Anal.*, vol. 9, pp. 291-312, 2002.
<figure><img src="image_2.png" /><figcaption>Figure 1: Effect of smoothing on the shock wave coming from an open medium and entering a porous medium $\epsilon_M = 0.05$. Linear path - Heavyside path- quadratic path- Riemann invariant $I_1^2 = \frac{U^2}{2} + \psi(\rho)$</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Effect of smoothing on the shock wave coming from an open medium and entering a porous medium $\epsilon_M = 0.05$. Linear path - Heavyside path- quadratic path - Momentum $\epsilon\rho U$</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Effect of smoothing on the shock wave coming from an open medium and entering a porous medium $\epsilon_M = 0.05$. Linear path - Heavyside path- quadratic path- Velocity</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Effect of smoothing on the shock wave coming from an open medium and entering a porous medium $\epsilon_M = 0.05$. Linear path - Heavyside path- quadratic path- Density</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Effect of smoothing on the rarefaction wave coming from an open medium and entering a porous medium $\epsilonpsilon_{M} = 0.05$. Linear path - Heavyside path-quadratic path- Riemann invariant $I_{1}^{2} = \frac{U^{2}}{2} + \psi(\rho)$</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Effect of smoothing on the rarefaction wave coming from an open medium and entering a porous medium $\epsilon_M = 0.05$. Linear path - Heavyside path- quadratic path- Momentum $\epsilon\rho U$</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Effect of smoothing on the rarefaction wave coming from an open medium and entering a porous medium $\epsilon_M = 0.05$. Linear path - Heavyside path- quadratic path- Velocity</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Effect of smoothing on the rarefaction wave coming from an open medium and entering a porous medium $\epsilon_M = 0.05$. Linear path - Heavyside path- quadratic path- Density</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Shock wave in an open medium hitting a porous medium $\epsilon_M = 0.05$ or a wall - Riemann invariant $I_1^2 = \frac{U^2}{2} + \psi(\rho)$</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Shock wave in an open medium hitting a porous medium $\epsilon_M = 0.05$ or a wall - Riemann invariant $I_1^1 = \rho\epsilon U$</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Shock wave in an open medium hitting a porous medium $\epsilon_M = 0.05$ or a wall - Density</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Shock wave in an open medium hitting a porous medium $\epsilon_M = 0.05$ or a wall - Velocity</figcaption></figure>