# Méthodes de domaine fictif et de raffinement adapté de maillage pour la mécanique des fluides numérique des composants nucléaires dans le contexte du projet NEPTUNE

Michel Belliard

► To cite this version:

Michel Belliard. Méthodes de domaine fictif et de raffinement adapté de maillage pour la mécanique des fluides numérique des composants nucléaires dans le contexte du projet NEPTUNE. 19e Congrès Français de Mécanique, Aug 2009, Marseille, France. papier: 507. hal-01053751

HAL Id: hal-01053751

https://hal.science/hal-01053751v1

Submitted on 1 Aug 2014

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>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Méthodes de domaine fictif et de raffinement adapté de maillage pour la mécanique des fluides numérique des composants nucléaires dans le contexte du projet NEPTUNE

M. BELLIARD

CEA Cadarache, DER/SSTH/LMDL, Bât. 238, 13108 Saint Paul-lez-Durance Cedex
Tél. 04 42 25 23 17, E-mail : michel.belliard@cea.fr

## Résumé :

Dans le cadre du projet NEPTUNE, on présente la simulations du bilan d'énergie d'un générateur de vapeur de centrale nucléaire (équation de convection - diffusion - réaction) par des méthodes de domaine fictif. Le domaine physique est inclus dans un domaine fictif plus grand et plus simple géométriquement. La dissociation de la géométrie des deux domaines présente de sérieux avantages : utilisation d'un maillage Cartésien (méthodes numériques performantes), traitement des frontières mobiles, ...

Plusieurs types de solveurs domaine fictif ont été développés pour des CL immergées générales (Dirichlet, Neumann ou Robin), en éléments finis $Q_1$ ou en volumes finis avec saut de flux et/ou de solution. Des techniques de raffinement de maillage permettent de réduire l'erreur de modélisation. Le solveur domaine fictif est couplé à un code industriel pour la fourniture des relations de fermeture. Les gains en temps CPU et mémoires sont soulignés ainsi que les potentialités pour le calcul parallèle.

## Abstract :

In the framework of the NEPTUNE project, we present the energy balance equation computation for a nuclear power plant steam generator (convection - diffusion - reaction equation) using a fictitious domain method. The physical domain is immersed in a bigger and geometrically simpler fictitious domain. As the shape of the computation domain is no longer correlated with the physical domain shape, it is possible to use performing numerical techniques for Cartesian meshes or for moving interfaces without remeshing the computation domain.

Various fictitious domain solvers are presented for general boundary conditions (Dirichlet, Neumann or Robin), using $Q_1$ finites elements or finites volumes with flux and/or unknown jumps. Adaptive mesh refinement techniques allow the reduction of the model error. The fictitious domain solver is coupled with an industrial code for the supply of the closing relations. CPU time and memory saving are presented and perspectives concerning parallel computations or Navier-Stokes extension are sketched.

**Mots clefs :** domaine fictif, volumes finis, saut de flux et de solution, éléments finis, pénalisation, générateur de vapeur

## 1 Introduction

Ces travaux ont été réalisés dans le cadre du projet NEPTUNE, financé par le CEA (Commissariat à l'Energie Atomique), EDF, l'IRSN (Institut de Radioprotection et de Sûreté Nucléaire) et AREVA-NP [1]. C'est un projet de codéveloppement en thermohydraulique nucléaire diphasique. Il a pour principal objectif le développement à horizon 2010-2015 de la nouvelle génération d'outils pour les différentes échelles de simulation (échelles locale, composant et système). La plate-forme NEPTUNE doit permettre à terme la réalisation de couplages inter-échelles et inter-disciplinaires. Le projet regroupe des activités de développement logiciel, de R&D en modèles physiques et méthodes numériques et de validation expérimentale). Le travail présenté ici s'inscrit dans les activités de R&D en méthodes numériques et concerne la résolution de l'équation du bilan d'énergie d'un générateur de vapeur (GV) de centrale nucléaire par des méthodes de domaine fictif. Cette équation est de type convection - diffusion - réaction.

La méthode des domaines fictifs a été initialement introduite par V. K. Saul'ev [2] au début des années 1960. Considérons le problème suivant, en dimension *d*, défini sur un domaine réel $\tilde{\Omega}$ (cf. Figure 1) : trouver $\tilde{u} \in$
$H^1(\tilde{\Omega})$ tel que

$$
\left\{
\begin{array}{ll}
-div(\tilde{\mathbf{a}} \cdot \nabla \tilde{u}) + div(\tilde{\mathbf{v}} \tilde{u}) + \tilde{b}(\tilde{u} - \tilde{u}_e) = \tilde{f} & \text{dans } \tilde{\Omega}, \\
\text{conditions aux limites} & \text{sur } \partial\tilde{\Omega} (:= \tilde{\Gamma} \cup \Sigma),
\end{array}
\right.
\qquad (1)
$$

avec $\tilde{\mathbf{a}} \in (L^\infty(\tilde{\Omega}))^{d \times d}$, $\tilde{\mathbf{v}} \in (L^\infty(\tilde{\Omega}))^d$, $\tilde{b} \in L^\infty(\tilde{\Omega})$, $\tilde{u}_e \in H^1(\tilde{\Omega})$ et $\tilde{f} \in L^2(\tilde{\Omega})$. Le tenseur de diffusion $\tilde{\mathbf{a}}$
et le coefficient de réaction $\tilde{b}$ vérifient l'hypothèse classique d'ellipticité ($\tilde{u}$ existe et est unique). La méthode
des domaines fictifs consiste à immerger le domaine réel dans un domaine fictif $\Omega$ de forme géométrique plus
simple. Le problème fictif se définit par : trouver $u \in H^1(\Omega)$ tel que

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

$$
\left\{
\begin{array}{ll}
-\mathrm{div}(\mathbf{a} \cdot \nabla u) + \mathrm{div}(\mathbf{v}u) + b(u - u_e) = f & \text{dans } \Omega, \\
\text{Conditions aux limites réelles pour } u & \text{sur } \partial\Omega \cap \partial\tilde{\Omega} := \tilde{\Gamma}), \\
\text{Conditions aux limites appropriées pour } u & \text{sur } \partial\Omega \setminus \partial\tilde{\Omega} := \Sigma),
\end{array}
\right.
\qquad (2)
$$

FIG. 1 – Définition du concept de la méthode des domaines fictifs (exemple d’un quart de disque immergé dans un carré).

où **a** ∈ (L^∞(Ω))^{d×d}, **v** ∈ (L^∞(Ω))^d, f ∈ L^2(Ω), u_e ∈ H^1(Ω) et b ∈ (L^∞(Ω)) tels que **a**, **b** vérifient l'hy-
pothèse classique d'ellipticité dans Ω (u existe et est unique). On a donc une frontière immergée Σ associée
aux conditions aux limites réelles (*Embedded Boundary conditions*, EBC).

Nous nous focalisons ici sur des conditions aux limites de type Robin $(-\tilde{\mathbf{a}} \cdot \nabla \tilde{u}) \cdot \mathbf{n} = \alpha_R(\tilde{u} - u_R) + g_R$, $\mathbf{n}$ vecteur unité normal extérieur) car elles nous permettent d'inclure de façon naturelle les conditions aux limites de type Neumann $(-\tilde{\mathbf{a}} \cdot \nabla \tilde{u}) \cdot \mathbf{n} = g_R$ ou de type Dirichlet $(\tilde{u} = u_R)$ par une technique de pénalisation.

Contrairement aux méthodes classiques, qui calquent la géométrie de la grille de calcul sur celle du domaine réel, l'avantage principal de l'approche en domaine fictif est de dissocier le maillage de calcul de la forme géométrique du domaine réel : mise en œuvre de méthodes performantes de résolution numérique en maillage Cartésien (précision, rapidité et faible coût en mémoire), facilité de gestion des interfaces mobiles (modification de la géométrie de la frontière immergée au lieu d'un remaillage du domaine de calcul), possibilité d'utiliser des méthodes se prêtant bien aux maillages Cartésien (analyse multirésolution par exemple), ... Ceci est particulièrement motivant pour la simulation de l'interaction fluide-structure ; les tubes des GV étant forte- ment sollicités par l'écoulement diphasique turbulent.

Dans la Section 2, l'équation du bilan d'énergie d'un GV est détaillée. Plusieurs méthodes de domaine fictif sont présentées pour des conditions limites immergées générales dans la Section 3, en éléments finis (EF, éléments $Q_1$) ou en volumes finis (VF, méthode avec saut de flux et/ou de solution [4]), utilisant un traitement en volume ou en surface. Elles permettent de modifier aussi peu que possible le schéma numérique initial d'un code donné. Pour augmenter la précision des calculs, on fait appel à des méthodes de raffinement local de maillage qui font l'objet de la Section 4(a). On rapporte dans la Section 5 des cas d'applications numériques académiques ou industriels (GV). Enfin on donne quelques remarques en guise de conclusion.

2 Bilan d'énergie d'un générateur de vapeur

La Figure 2(a) présente le schéma de fonctionnement d'un GV. Le calcul concerne l'évaporateur (partie où le fluide secondaire est en contact avec le fluide primaire). De fait, les simulations industrielles des générateurs de vapeur utilisent des techniques d'homogénéisation (prise de moyenne, cf. Figure 2(b)) et décrivent le fluide diphasique secondaire comme un mélange équivalent (liquide plus vapeur) se déplaçant dans un milieu modélisé comme poreux. La porosité $\beta$ est le rapport entre le volume du mélange et celui de la cellule d'homogénéisation. Deux modélisations basées sur ces techniques sont fréquemment rencontrées dans le contexte des simulations industrielles : *Homogeneous Equilibrium Model* et *Homogeneous Relaxation Model* [5]. L'expression de l'équation du bilan d'énergie pour le mélange, où $H$ est l'enthalpie spécifique ($kJ.kg^{-1}$), est de la forme suivante (code CEA GENEPI [6]) :

$$
\beta \rho \partial_t H + \beta \mathbf{G} \cdot \nabla H - \operatorname{div}(\beta \chi_T \nabla H) = Q - \operatorname{div}(\beta x(1-x) \rho \mathcal{L} \mathbf{V}_R), \quad (3)
$$

où

- $\rho(H, P)$ est la densité ($kg.m^{-3}$) fonction de la pression $P$ ($Pa$) et de $H$,

- $\chi_T = a|\mathbf{G}|L$ est le coefficient de diffusion turbulente ($kg.m^{-1}.s^{-1}$) avec $a$ le coefficient de Schlichting [7],
  **G** le flux de masse ($kg.m^{-2}.s^{-1}$), $L$ une longueur typique de tourbillon et $|\cdot|$ la norme Euclidienne,
- $Q(H, P)$ est le terme source de chaleur ($W.m^{-3}$) donné par la résolution de l'équation du bilan d'énergie du fluide primaire,

- $V_R(H, P)$ est la vitesse relative (vitesse de la vapeur moins vitesse du liquide, $m.s^{-1}$) [8],

- $x(H, P)$ est le titre statique ($\equiv \frac{H-H_{ls}}{L}$) avec $H_{ls}(P)$ l'enthalpie du liquide à saturation ($KJ.kg^{-1}$) et $L(P)$ la chaleur latente ($J.kg^{-1}$). Le terme $\beta x(1-x)\rho LV_R$ est appelé terme de dérive.

<figure><img src="image_3.png" /><figcaption>FIG. 2 – Schéma d’un générateur de vapeur et cellule d’homogénéisation.</figcaption></figure>

Dans le code GENEPI, la discrétisation spatiale de la formulation faible de l'équation (3) est établie par des EF $Q_1(H)$ et $\beta$ et par des éléments finis $Q_0(P, \rho, G, \chi_T)$ et $Q$. La discrétisation en temps est faite par la méthode de Crank-Nicholson. La correction BTD [9] est appliquée pour empêcher de fausses oscillations introduites par la discrétisation centrée du terme de convection. A chaque pas de temps, le système linéaire est lissé par une vingtaine d'itérations CGS préconditionnées ILLU.

## 3 Méthodes de domaine fictif

La confrontation des équations (1) et (3), écrites sous forme conservatrice, permet d'identifier :

$$ \tilde{u} = H^n, \quad \tilde{a} = \beta\chi_T, \quad \tilde{b} = \frac{\beta\rho}{\delta t}, \quad \tilde{u}_e = H^{n-1}, \quad \tilde{v} = \beta G $$

où $\delta t$ est le pas de temps et $n$ est l'indice en temps. Il est nécessaire de définir plus précisement les coefficients de l'équation (2). Après discrétisation spatiale, nous prenons naturellement pour $\tilde{\Omega}_h$ :

$$ a|_{\tilde{\Omega}_h} = \tilde{a}|_{\tilde{\Omega}_h}, \quad v|_{\tilde{\Omega}_h} = \tilde{v}|_{\tilde{\Omega}_h}, \quad b|_{\tilde{\Omega}_h} = \tilde{b}|_{\tilde{\Omega}_h}, \quad f|_{\tilde{\Omega}_h} = \tilde{f}|_{\tilde{\Omega}_h}. $$

De même pour le domaine externe discrétisé $\Omega_{e,h}$, nous choisissons :

$$ a|_{\Omega_{e,h}} = \mathbf{Id}, \quad v|_{\Omega_{e,h}} = \mathbf{0}, \quad b|_{\Omega_{e,h}} = \mathbf{0}, \quad f|_{\Omega_{e,h}} = \mathbf{0}. $$

On prend des conditions aux limites de type Robin homogène ($\alpha_R = u_R = g_R = 0$) sur $\partial\Omega_h \setminus \partial\tilde{\Omega}_h$.

Les diverses méthodes de domaines fictifs se différencient par la façon d'imposer les EBC dans l'équation (2). Typiquement, on peut modéliser la frontière immergée $\Sigma$ soit par une frontière approchée épaisse $\omega_{h,\Sigma}$ (volumule, cellules traversées par $\Sigma$), soit par une frontière approchée fine $\Sigma_h$ (surface). L'étude du problème de transmission sous forme faible entre $\tilde{\Omega}$ et $\Omega_e$ fait apparaître une intégrale de saut de flux surfacique sur $\Sigma$ (I. Ramière et al. [10]) qui s'écrit $\int_{\Sigma} \{\alpha_R(\tilde{u}-u_R)+g_R+(\tilde{\mathbf{v}}.\mathbf{n})\tilde{u}\}$ ds si les flux diffusif et convectif sont nuls côté $\Omega_e$. On égale cette intégrale à une intégrale de volume $\int_{\omega_{h,\Sigma}} (\alpha_R(u_h-u_R)+g_R+(\mathbf{v}.\mathbf{n})u_h) \frac{dx}{\epsilon_h}$ ou de surface $\int_{\Sigma_h} (\alpha_R(u_h-u_R)+g_R+(\mathbf{v}.\mathbf{n})u_h) \frac{ds}{\epsilon_h}$ sur la frontière approchée, suivant le type de modélisation retenue. $\epsilon_h$ est un coefficient introduit pour assurer au mieux cette égalité, n'ayant pas la même valeur suivant le modèle de frontière immergée.
## 3.1 Modélisation épaisse de la frontière immergée

Dans le cas d'une frontière approchée épaisse, nous imposons l'EBC par des termes de réaction et de source en volume, les termes b et f dans (2) [10] :

$$a|_{\omega_{h,\Sigma}} = \tilde{a}|_{\omega_{h,\Sigma}}, \quad v|_{\omega_{h,\Sigma}} = \tilde{v}|_{\omega_{h,\Sigma}}, \quad b|_{\omega_{h,\Sigma}} = \tilde{b}|_{\omega_{h,\Sigma}} + \frac{\alpha_R + \mathbf{v} \cdot \mathbf{n}}{\epsilon_h}, \quad f|_{\omega_{h,\Sigma}} = \tilde{f}|_{\omega_{h,\Sigma}} + \frac{\alpha_R u_R - g_R}{\epsilon_h},$$

Ce modèle est utilisable aussi bien en volumes qu'en éléments finis. Soit $0 < \eta << 1$ un paramètre de pénalité réel qui tend vers zéro. En prenant $\frac{\alpha_R}{\epsilon_h} = \frac{1}{\eta}$ et $\frac{g_R}{\epsilon_h} = \frac{\mathbf{v} \cdot \mathbf{n}}{\epsilon_h} = 0$, nous pouvons imposer des conditions aux limites de type Dirichlet ($\tilde{u} = u_R$) par une pénalisation $L^2$ de l'interface approchée (en supposant tous les autres coefficients négligeables par rapport à $\frac{1}{\eta}$).

## 3.2 Modélisation fine de la frontière immergée

La méthode des volumes finis se prête bien au cas d'une frontière approchée fine. La frontière approchée $\Sigma_h$ suit les faces des volumes de contrôle. L'intégrale de saut de flux du problème de transmission est prise en compte directement dans l'expression du flux de diffusion aux faces de $\Sigma_h$ [11].

On utilise un modèle récent de fracture proposé par Ph. Angot [4], qui lie les sauts de flux de diffusion $(\phi(u).\mathbf{n} = -(a\nabla u).\mathbf{n})$ et de solution sur $\Sigma$. Au niveau continu, pour une fonction $\psi$ (flux ou solution) dans $H^1(\tilde{\Omega} \cup \Omega_e)$, on choisit de noter $\psi_{\Sigma}^{-}$ et $\psi_{\Sigma}^{+}$ les traces de $\psi|_{\tilde{\Omega}}$ et $\psi|_{\Omega_e}$ sur chaque face de $\Sigma$ respectivement, $\tilde{\psi}|_{\Sigma} = (\psi_{\Sigma}^{+} + \psi_{\Sigma}^{-})/2$ la moyenne arithmétique des traces de $\psi$, et $[[\psi]]_{\Sigma} = (\psi_{\Sigma}^{+} - \psi_{\Sigma}^{-})$ le saut des traces de $\psi$ sur $\Sigma$ orienté par $\mathbf{n}$. A la définition du problème fictif (2), on rajoute les relations suivantes :

$$[(a\nabla u).\mathbf{n}]_{\Sigma} = \alpha \bar{u}|_{\Sigma} - q \quad \text{sur } \Sigma, \quad (4)$$

$$\frac{(a\nabla u) \cdot \mathbf{n}|_{\Sigma}}{\alpha \bar{u}|_{\Sigma}} = \beta [[u]]_{\Sigma} - g \quad \text{sur } \Sigma, \quad (5)$$

où $g$ et $q$ sont données dans $L^2(\Sigma)$ et les coefficients de transfert $\alpha, \beta \ge 0$ sur $\Sigma$ sont des fonctions mesurables et bornées vérifiant les hypothèses classiques d'ellipticité. Quand $\alpha = g = q = 0$ et $\beta \to \infty$, nous sommes dans le cas du problème de transmission classique avec $u \in H^1(\Omega)$ (i.e. continuité du flux et de la solution). Au niveau discret, on choisi $(\phi_h(u_h) \cdot \mathbf{n})_{\Sigma_h}^{-}$ égale au flux de Robin, avec $\phi_h(u_h)$ le flux numérique. Dans la méthode "pas de contrôle extérieur", le choix $\alpha = 4\beta$ garantit que $(\phi_h \cdot \mathbf{n})_{\Sigma_h}^{-}$ ne dépend que de $(u_h)_{\Sigma_h}^{-}$ et non de $(\phi_h \cdot \mathbf{n})_{\Sigma_h}^{+}$ et de $(u_h)_{\Sigma_h}^{+}$. On prend : $\alpha = 4\beta = 2\alpha_R$ et $g - \frac{q}{2} = g_R - \alpha_R u_R$.

# 4 Raffinement adapté du maillage

En utilisant un schéma d'éléments finis $Q_1$, I. Ramière a montré que l'approche décrite ci-dessus conduit à une méthode en $O(h)$ [10, 12]. Cela est dû à la position approximative de la frontière dans le maillage cartésien. Des techniques de raffinement de maillage (Adaptive Mesh Refinement, AMR), même si elles n'augmentent pas l'ordre de convergence, permettent de réduire la valeur de l'erreur. La méthode Local defect correction de [13] est mise en œuvre en EF [10] et la méthode Flux Interface Correction de Ph. Angot et al. [14] modifiée par I. Ramière et al. (méthode FIC-EBC) est mise en œuvre en VF [11]. Le raffinement de maillage est adapté de façon à raffiner les cellules interceptées par $\Sigma$. Les deux méthodes ont en commun l'utilisation d'un algorithme multi-grilles entre la grille de base et celles de raffinement, cf. Figure 3. Pour la méthode FIC-EBC, un patch de raffinement est crée autour de chaque cellule interceptée par $\Sigma$, cf. Figure 4(a). Le bilan de flux d'une cellule grossière $K$ est corrigé en tenant compte des flux du patch sur les faces intérieures $\sigma$ de $K$, cf. Figure 4(b). Pour un niveau de raffinement, le calcul sur chaque patch est indépendant et pourrait être exécuté en parallèle. Tou- tefois, un préconditionnement, sur chaque grille, par décomposition de domaine permet d'accélérer la conver- gence de l'algorithme, mais nécessite des échanges entre processeurs. Il est également important de bien gérer l'échange d'informations vers la grille grossière de base pour éviter un goulot d'étranglement (décomposition de domaine de la grille de base et regroupement des patches par sous-domaine).

# 5 Applications numériques

## 5.1 Cas tests académiques

Dans [10, 11] de nombreux cas tests numériques, pour lesquels la solution analytique d'une équation de convec- tion - diffusion - réaction est connue, sont présentés en EF et VF avec AMR. Rappelons qu'ils confirment que la méthode est d'ordre $O(h_f)$ où $h_f$ est le pas d'espace de la grille la plus fine, cf. Figure 5(a). Ils mettent également en évidence le gain conséquent en temps CPU que procure l'AMR, cf. Figure 5(b). La confirmation de ces résultats sur des simulations d'intérêts industriels fait l'objet de la Section suivante.
<figure><img src="image_4.png" /><figcaption>FIG. 3 – Algorithme multi-grilles pour AMR.</figcaption></figure>

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

<figure><img src="image_6.png" /><figcaption>FIG. 4 – Méthode AMR FIC-EBC. Patch de raffinement et correction des flux.</figcaption></figure>

## 5.2 Un cas test industriel

La validation sur un cas test industriel est la simulation 3D en domaine fictif du bilan d'énergie de CLOTAIRE, maquette de GV du CEA. Les EBC sont mixtes (Dirichlet et Neumann). Le solveur est couplé au code GENEPI pour la fourniture des relations de fermeture physique. La convergence en maillage, sans AMR, est étudiée avec des maillages du domaine fictif en forme de parallélépipède, comportant de 7 200 à 460 800 cellules, et englobant le domaine physique (i.e. l'évaporateur du GV), cf. Figure 6(a). Que la discrétisation spatiale soit en EF ou en VF, l'erreur introduite par la méthode converge bien en norme $L^2$ en $O(h)$ en moyenne vis-à-vis d'une solution numérique à 648 388 cellules en maillage adapté à la géométrie, cf. Figure 6(b). L'erreur relative est réduite à $10^{-3}$ environ dès que le maillage a un nombre raisonnable de cellules. L'écart relatif mesuré sur des quantités physiques, telle l'énergie échangée, ne dépasse pas $5 \cdot 10^{-3}$.

Globalement, la méthode de domaine fictif en EF est légèrement plus précise que celles en VF, mais les écarts sont faibles. Pour une EBC de type Dirichlet, la méthode "pas de contrôle de l'extérieur" est même légèrement plus précise. Il en est de même lorsque le maillage est grossier. On peut également souligner le gain en temps CPU (30%) vis-à-vis du calcul de référence à 648 388 cellules pour le même pas de discrétisation.

## 6 Conclusion

Des méthodes de domaine fictif en EF et en VF ont été testées avec succès sur des cas test académiques et industriels dans le contexte des générateurs de vapeur. La convergence en maillage de la méthode est illustrée. Le raffinement du maillage Cartésien accroît de manière importante la précision des calculs. A pas d'espace similaires, les résultats de calcul sont aussi précis que ceux obtenus par de façon classique avec un maillage adapté à la géométrie du GV. On peut même dire que pour un maillage 10 fois plus grossier que celui de référence, l'erreur relative de calcul est déjà suffisamment petite pour nos applications (environ $10^{-3}$).

## Références

[1] Guelfi A., Bestion D., Boucker M., Boudier P., Fillion P., Grandotto M., Hérard J.-M., Hervieu E., and Péturaud P. NEPTUNE - A New Software Platform for Advanced Nuclear Thermal Hydraulics. Nuclear Science and Engineering, 156, 281-324, July 2007.

[2] Saul'ev V. On the solution of some boundary value problems on high performance computers by fictitious domain method. Siberian Math. Journal, 4(4), 912-925, 1963 [in Russian].

[3] Marchuk G. Methods of Numerical Mathematics. Springer-Verlag, New York, USA, 1982 (1st ed. 1975).
<figure><img src="image_7.png" /></figure>

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

FIG. 5 – Cas tests académiques, EBC Robin : convergence de l’erreur relative en norme $L^2$ en fonction du pas d’espace et du temps CPU. Modélisation fine de la frontière immergée. Méthode de domaine fictif en VF “pas de contrôle de l’extérieur” [11].

[4] Angot P. A model of fracture for elliptic problems with flux and solution jumps. C.R. Acad. Sci. Paris, Ser I Math., 337 (6) 425-430, 2003.

[5] Clerc S. Numerical simulation of the Homogeneous Equilibrium Model for Two-phase Flows. Journal of Computational Physics, 161(1), 354–375, 2000.

[6] Grandotto M. and Obry P. Calculs des écoulements diphasiques dans les échangeurs par une méthode aux éléments finis. Revue Européenne des Éléments Finis, 5(1), 53–74, 1996 [in French].

[7] Schlichting H. Boundary Layer Theory. Mac Graw Hill, New York, USA, 1968.

[8] Zuber N. and Findlay J. Average Volumetric Concentration in Two-Phase Flow Systems. J. Heat Transfer, 87(4), 453–468, 1965.

[9] Gresho P., Chan S., Lee R. L., and Upson C. A modified Finite Element Method for Solving the Time-Dependent Incompressible Navier Stokes Equations (Part 1 : Theory). Int. J. Num. Methods in Fluids, 4(6), 557–598, 1984.

[10] Ramière I., Angot P., and Belliard M. A fictitious domain approach with spread interface for elliptic problems with general boundary conditions. Computer Methods in Applied Mechanics and Engineering, 196(4-6), 766–781, 2007.

[11] Ramière I., Angot P., and Belliard M. A general fictitious domain method with immersed jumps and non-conforming structured mesh. Journal of Computational Physics, 225(2), 1347–1387, 2007.

[12] Ramière I. Convergence analysis of the Q1-finite element method for elliptic problems with non-boundary-fitted meshes. International Journal for Numerical Methods in Engineering, 75(9), 1007–1052, 2008.

[13] Hackbusch W. Local Defect Correction Method and Domain Decomposition Techniques, volume 5 of Computing Suppl., pages 89–113. Springer-Verlag (Wien), 1984.

[14] Angot P., Caltagirone J. P., and Khadra K. Une méthode adaptative de raffinement local : la correction du flux à l'interface. C. R. Acad. Sci. Paris, Sér. I, 315, 739–745, 1992.
<figure><img src="image_9.png" /></figure>

(a) Coupe horizontale du maillage du domaine fictif (en rouge). Le maillage bleu représente un maillage adapté à la géométrie du GV.

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

(b) Erreur relative en norme $L^2$ le long de profils verticaux d'enthalpie du mélange en fonction du pas d'espace $h$.

FIG. 6 – Simulation du bilan d'énergie de CLOTAIRE par la méthode de domaine fictif en VF “pas de contrôle de l'extérieur” [11].