# Simulation de l'Accident de Dimensionnement du Confinement dans un réacteur de IVe génération avec une approche parallèle hybride dans EUROPLEXUS

Vincent Faucher, Thierry Gautier

► To cite this version:

Vincent Faucher, Thierry Gautier. Simulation de l'Accident de Dimensionnement du Confinement dans un réacteur de IVe génération avec une approche parallèle hybride dans EUROPLEXUS. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717086

HAL Id: hal-01717086

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

Submitted on 25 Feb 2018

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 CC0 1.0 - Universal - International License
# Simulation de l'Accident de Dimensionnement du Confinement dans un réacteur de IVe génération avec une approche parallèle hybride dans EUROPLEXUS

V. Faucher<fnref n="1" />*, T. Gautier<fnref n="2" />

**Résumé** — Cette contribution les derniers développements méthodologiques mis en œuvre dans le programme EUROPLEXUS pour la simulation de l'Accident de Dimensionnement du Confinement dans un réacteur de IVe génération. Cette simulation requiert des méthodes avancées en interaction fluide-structure avec interfaces et grands déplacements, ainsi qu'un environnement d'exécution parallèle efficace et extensible.

**Mots clés** — Interaction fluide-structure, transitoire accidentel, calcul parallèle

## 1 Eléments de contexte

L'Accident de Dimensionnement du Confinement est l'accident de référence pour les réacteurs nucléaire de VIe génération, pour lesquels la tenue aux accidents graves doit être démontrée dès la conception.

Il s'agit d'un transitoire complexe impliquant l'interaction de plusieurs fluides et interfaces contenus dans une cuve métallique avec des structures internes immergées sujettes à de grands déplacements. La géométrie considérée est un premier remontage du prototype ASTRID conçu au CEA et présenté sur la Figure 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Modèle de calcul pour la simulation de l'Accident de Dimensionnement du Confinement</figcaption></figure>

Les éléments méthodologiques nécessaires à la simulation d'une telle situation accidentelle, en vue de son intégration dans un processus de conception dans une unité d'ingénierie, sont les suivants :

<footnote n="1">1 CEA, DEN, DANS, DM2S, SEMT, DYN, vincent.faucher@cea.fr</footnote>

<footnote n="2">2 INRIA, EPI MOAIS, Laboratoire d'Informatique de Grenoble, thierry.gautier@imag.fr</footnote>

<footnote>* Auteur correspondant</footnote>
* ✓ prise en compte robuste et précise des écoulements avec interface en régime transitoire brutal,

* ✓ interaction fluide-structure en grands déplacements,

* ✓ résolution parallèle efficace et extensible.

La présente contribution présente la stratégie utilisée dans le programme EUROPLEXUS (EPX) pour répondre à ces impératifs, en insistant particulièrement sur la stratégie de résolution parallèle hybride mise en œuvre en collaboration avec l'INRIA dans le cadre du projet RePDyn (ANR-COSI-09-001).

# 2 Méthodologie mise en œuvre

## 2.1 Ecoulement avec interfaces et interaction fluide-structure

### 2.1.1 Approche anti-diffusive pour la représentation des interfaces

La représentation des interfaces repose sur une extension de l'approche VOFIRE [1][2][3] pour les écoulements liquide-gaz. Elle s'appuie sur une représentation homogène du fluide, caractérisée par les concentrations locales de chacun des composants, liquide ou gazeux, considérés. Cette approche permet de s'affranchir de la reconstruction des interfaces, pour un gain manifeste en robustesse pour des géométries complexes en représentation tridimensionnelle.

Le suivi des interfaces est assuré par un blocage de la diffusion numérique des concentrations (cf. Figure 2), limitant l'épaisseur des zones de mélange à la dimension d'une maille fluide.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Diffusion numérique des concentrations avec un problème modèle d'explosion sous eau</figcaption></figure>

Le principe anti-diffusif est illustré en représentation mono-dimensionnelle par le transport à vitesse constante de la fonction de Heavyside avec un schéma Volumes Finis (cf. Figure 3).

<figure><img src="image_4.png" /><figcaption>Fig. 3. Transport de la fonction de Heavyside dans un champ de vitesse uniforme</figcaption></figure>
A un instant donné de l'échelle d'intégration temporelle, de pas de temps $\Delta t$, le front de la fonction se trouve dans la cellule j. La question porte alors sur le choix de la valeur discrète $h_{j+1/2}$ de la fonction h à considérer au niveau de la frontière $[j+1/2]$ de la cellule j pour écrire le flux de h à travers la frontière, intervenant dans le bilan de la fonction h dans la cellule entre l'instant courant et le suivant.

Cette situation présente deux cas limite.

1. Décentrement total amont (upwind) : $h_{j+1/2} = h_j$

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

Le flux est non-nul dès que $h_j$ est non-nul, si bien qu'au niveau discret, la fonction se propage dans la cellule $j+1$ avant que le front analytique ait atteint la frontière $[j+1/2]$, ce qui est à l'origine de la diffusion numérique.

2. Décentrement total aval (downwind dans la suite du texte) : $h_{j+1/2} = h_{j+1}$

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

Le flux est toujours nul. Tant que le front théorique est localisé dans la cellule j, la diffusion numérique est bloquée. La solution devient fausse lorsque le front passe dans la cellule $j+1$, car la fonction ne peut être transportée à travers la frontière $[j+1/2]$. Il en résulte une augmentation instable de h dans la cellule j.

Ces deux cas limite introduisent le concept de base établi par Després et Lagoutière : la recherche d'un schéma avec décentrement aval, associé à des contraintes lui permettant de rester stable lorsqu'une discontinuité traverse la frontière d'une cellule (Downwind Scheme with Constraints).

Les contraintes à considérer sont les suivantes :

$$ \text{Consistance pour le flux : } m_{j+1/2} = \min(h_j, h_{j+1}) \le h_{j+1/2} \le \max(h_j, h_{j+1}) = M_{j+1/2} \quad [\text{Eq. 1a}] $$

$$ \text{Stabilité pour l'advection : } m_{j-1/2} = \min(h_{j-1}, h_j) \le h_j^* \le \max(h_{j-1}, h_j) = M_{j-1/2} \quad [\text{Eq. 1b}] $$

$$ \text{avec } h^* = h_j + \mathbf{u} \frac{\Delta t}{\Delta x} (h_{j-1/2} - h_{j+1/2}) $$

En prenant en compte l'encadrement [Eq. 1a] pour le flux $h_{j-1/2}$, une condition suffisante pour la stabilité est donnée par les deux inégalités :

$$ \begin{cases} m_{j-1/2} \le h_j + \mathbf{u} \frac{\Delta t}{\Delta x} (m_{j-1/2} - h_{j+1/2}) \\ h_j + \mathbf{u} \frac{\Delta t}{\Delta x} (M_{j-1/2} - h_{j+1/2}) \le M_{j-1/2} \end{cases} \quad [\text{Eq. 2}] $$

On en déduit deux intervalles de confiance pour le $h_{j+1/2}$ :
$$ m_{j+1/2} \le h_{j+1/2} \le M_{j+1/2} $$
$$ M_{j-1/2} + \frac{\Delta x}{\mathbf{u}\Delta t} (M_{j-1/2} - h_j) \le h_{j+1/2} \le m_{j-1/2} + \frac{\Delta x}{\mathbf{u}\Delta t} (m_{j-1/2} - h_j) \quad [\text{Eq. 3}] $$

L'intersection de ces deux intervalles est forcément non-vide, car le choix *upwind* pour le flux est toujours stable. Le traitement anti-diffusif consiste alors à choisir pour $h_{j+1/2}$ la valeur la plus proche dans cette intersection de la valeur *downwind* $h_{j+1}$. Le résultat sur le problème modèle de la Figure 2 est donné sur la Figure 4.

<figure><img src="image_7.png" /><figcaption>Fig. 4. Traitement anti-diffusif sur les concentrations pour localiser les interfaces</figcaption></figure>

### 2.1.2 Interaction fluide-structure

Pour l'interaction fluide-structure, deux approches sont utilisées simultanément : une approche avec maillage coïncidents pour l'interaction avec la cuve sur la frontière du domaine fluide et une approche avec maillage découpés pour les structures internes immergées. Dans tous les cas, les grands déplacements de structure pris en compte.

La gestion des connexions cinématiques entre fluide et structure exclue tout paramètre numérique arbitraire :

* les forces d'interaction sont obtenues à l'aide de multiplicateurs de Lagrange dans le cas des maillages conformes, selon une procédure de dualisation classique des contraintes cinématiques,

* une approche directe maître-esclave est utilisée pour les structures immergées.

En considérant la configuration géométrique donnée par la Figure 5, les équations d'interaction entre des nœuds fluide et des nœuds de structure immergée sont :
<figure><img src="image_8.png" /><figcaption>Fig. 5. Interaction entre un fluide et une structure immergée</figcaption></figure>

$$ \left( \mathbf{u}_F - \sum_{i \in D} \alpha_i \mathbf{u}_{S_i} \right) \cdot \mathbf{n}_S = 0 \qquad [\text{Eq. 4}] $$

où $\mathbf{u}_F$ est la vitesse nodale du fluide et $\mathbf{u}_S$ celle de la structure,

$D$ désigne le domaine d'influence, comprenant soit 1 nœud (domaine nodal), soit 2 nœuds (domaine linéique), soit 3 ou 4 nœuds (domaine surfacique),

$\{\alpha_i\}$ est un jeu de coefficients de pondération déterminés géométriquement, donnant localement la vitesse de la structure en fonction des vitesses nodales pour les nœuds du domaine d'influence.

L'approche maître-esclave consiste à calculer la force normale agissant sur la structure sans recourir à la résolution d'un système linéaire portant sur des multiplicateurs. Elle est obtenue en considérant le nœud fluide comme esclave de la structure dans la direction normale, si bien qu'on ajoute dans les équations d'équilibre des nœuds de structure du domaine d'influence les forces normales agissant sur le fluide et la masse associée au nœud fluide, en prenant en compte la pondération introduite dans l'équation [Eq. 4].

Une fois l'accélération calculée et la vitesse mise à jour pour les nœuds de structure, la vitesse normale du fluide est déduite directement de celle de la structure. En l'absence d'autre contrainte cinématique portant sur les nœuds concernés, les équations correspondantes sont les suivantes :

$$ (a) \quad \mathbf{f}_{S_i} = \mathbf{f}_{S_i}^{\text{ext}} - \mathbf{f}_{S_i}^{\text{int}} + \alpha_i ([\mathbf{f}_{F}^{\text{ext}} - \mathbf{f}_{F}^{\text{int}}] \cdot \mathbf{n}_S) \mathbf{n}_S $$

$$ (b) \quad m_{S_i}^{\text{tot}} = m_{S_i} + \alpha_i m_F $$

$$ (c) \quad \mathbf{g}_{S_i} = \frac{\mathbf{f}_{S_i}}{m_{S_i}^{\text{tot}}}; \mathbf{g}_F = \frac{\mathbf{f}_F^{\text{ext}} - \mathbf{f}_F^{\text{int}}}{m_F} \qquad [\text{Eq. 5}] $$

$$ (d) \quad \mathbf{u}_{S_i}^{+} = \mathbf{u}_{S_i} + \Delta \mathbf{t} \mathbf{g}_{S_i}; \tilde{\mathbf{u}}_{F}^{+} = \mathbf{u}_{F} + \Delta \mathbf{t} \mathbf{g}_{F} $$

$$ (e) \quad \mathbf{u}_{F}^{+} = \tilde{\mathbf{u}}_{F}^{+} - (\tilde{\mathbf{u}}_{F}^{+} \cdot \mathbf{n}_{S}) \mathbf{n}_{S} + \sum_{i \in D} \alpha_{i} (\mathbf{u}_{S_{i}}^{+} \cdot \mathbf{n}_{S}) \mathbf{n}_{S} $$

où $\mathbf{f}_k^{\text{ext}}$ et $\mathbf{f}_k^{\text{int}}$ sont respectivement les forces nodales externes et internes pour l'entité $k$ (structure S ou fluide F),
$m_k$ est la masse nodale pour l'entité $k$,

$\mathbf{u}_k$ et $\mathbf{u}_k^+$ sont respectivement les vitesses nodales pour l'entité $k$ à l'instant courant et à l'instant suivant dans le processus d'intégration temporelle.

## 2.1. Résolution parallèle hybride

### 2.1.3 Solveur à mémoire distribué optimisé pour la gestion des liaisons cinématiques

Pour permettre une exécution parallèle sur un grand nombre de cœurs, EPX est doté d'un solveur parallèle à mémoire distribué reposant sur une décomposition de domaine automatique. Dans le cadre du projet RePDyn, ce solveur, classiquement conçu pour l'accélération des tâches élémentaires (calcul local de forces internes et de flux de quantités eulériennes) a été optimisé pour la prise en compte des nombreuses contraintes cinématiques qui caractérisent le problème posé de dynamique avec couplage fluide-structure.

L'effort porte principalement sur les liaisons traitées par l'intermédiaire de multiplicateurs de Lagrange (contact unilatéral entre éléments de structure notamment pour le cas considéré au § 1). Cette approche introduit la solution nécessaire d'un système de la forme [Eq. 6], les termes de la matrice faisant intervenir lors de son écriture de nombreux couplages entre les sous-domaines et donc un volume de communications important entre les processus MPI.

$$ \mathbf{H}\mathbf{A} = \mathbf{B} \qquad [\text{Eq. 6}] $$

La résolution directe de ce système est remplacée par une résolution itérative évitant l'assemblage du système, le remplaçant par l'évaluation du résidu $\mathbf{H}\mathbf{A}_k - \mathbf{B}$, où $\mathbf{A}_k$ est le $k$-ième itéré de l'algorithme. Cette évaluation peut être réalisée avec des opérations majoritairement locales sur les sous-domaines via la procédure suivante :

1. évaluation locale des forces correspondant au vecteur de multiplicateurs $\mathbf{A}_k$ :

$$ \mathbf{F}_k = {}^t\mathbf{C}\mathbf{A}_k \qquad [\text{Eq. 7a}] $$

où $\mathbf{C}$ est la matrice des liaisons cinématiques écrites sur un sous-domaine donné, faisant intervenir des inconnues locales et éventuellement des inconnues appartenant à d'autres sous-domaines (inconnues dites remote),

2. sommation des forces sur les inconnues remote : seule étape de communication, nettement moins volumineuse que la construction explicite de l'opérateur $\mathbf{H}$,

3. évaluation locale de la contribution de chaque sous-domaine au résidu :

$$ \mathbf{R}_k = \mathbf{C}\mathbf{M}^{-1}\mathbf{F}_k - \mathbf{B} \qquad [\text{Eq. 7b}] $$

Une telle approche itérative est accompagnée nécessairement d'un préconditionnement pour accélérer la convergence, conçu pour être construit uniquement à partir d'opérations locales à chaque sous-domaine. La nature du préconditionneur utilisé est déduite de la structure de l'opérateur $\mathbf{H}$ (cf. Figure 6), formé de blocs diagonaux correspondant à des groupes de liaisons couplées, chacun de ces blocs contenant une partie correspondant aux couplages entre liaisons écrites sur un même sous-domaine (représentée en bleu et évaluée localement sur un sous-domaine) et une partie correspondant aux couplages entre liaisons écrites sur des sous-domaines différents (représentée en rouge et à l'origine des communications inter-domaines).
<figure><img src="image_9.png" /><figcaption>Fig. 6. Structure de l'opérateur de calcul des multiplicateurs de Lagrange</figcaption></figure>

Si la décomposition de domaine est correcte, les blocs rouges sont de taille modérée et le préconditionneur utilisé consiste à ignorer ces blocs. L'efficacité de cette stratégie (dite solveur GPCG dans EPX) est démontrée sur un cas difficile de contact entre rotor et stator, où on montre qu'elle permet d'abaisser et de préserver le coût par cycle en temps (cf. Figure 7), couplée à une stratégie de reconstruction périodique de la décomposition de domaine assurant la qualité du préconditionneur.

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

Contact rotor-stator avec arrachement d'ailettes (modèle ONERA). Grande variation de la structure du problème de contact en cours de calcul.

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

Impact de la réactualisation de la décomposition de domaine et du solveur GPCG pour les liaisons sur la réduction et la maîtrise du coût par cycle en cours de calcul pour 16 sous-domaines.

<caption>Fig. 7. Validation de l'optimisation du solveur parallèle à mémoire distribuée pour la gestion des liaisons cinématiques</caption>

### 2.1.4 Parallélisme à mémoire partagée complémentaire et approche hybride

Pour répondre à l'évolution actuelle des architectures pour les supercalculateurs, composés de grappes de nœuds à mémoire partagée, une approche dédiée à ces derniers est intégrée à EPX en complément de l'approche distribuée.

Dans le cadre de l'approche hybride proposée, le parallélisme à mémoire partagée est mis en œuvre via la librairie KAAPI, développée au laboratoire INRIA/LIG [5], assurant l'équilibre de la charge entre les *threads* par vol de travail. Cet équilibre, indispensable pour les systèmes dynamiques considérés dans le présent document, est démontré sur la Figure 8 sur un cas pénalisant de crash de missile sur un support rigide.
<figure><img src="image_12.png" /></figure>

Visualisation du missile à l'instant initial et après 1000 cycles en temps : apparition localisée de déformation plastique et de contacts

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

Equilibre parfait du coût par threads pour les boucles les plus consommatrices (83% du temps total), avec une distribution non-uniforme des itérations grâce au vol de travail.

<caption>Fig. 8. Validation de l'optimisation du solveur parallèle à mémoire distribuée pour la gestion des liaisons cinématiques</caption>

Cette stratégie a permis d'obtenir un parallélisme extensible de 32 à 1024 cœurs pour des systèmes fluide-structure fortement couplés tels que le calcul d'ADC faisant l'objet de ce document et intégrant toutes les stratégies présentées. Le parallélisme résultant a été testé dans le cadre d'un PRACE Preparatory Access (PA 0624) sur le supercalculateur TGCC/Curie.

Références

[1] B. Després, F. Lagoutière, Numerical resolution of a two-component compressible fluid model with interfaces, Progress in Computational. Fluid Dynamics, Vol 7(6), Pages 295-310, 2007.

[2] S. Kokh, F. Lagoutière, An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of a five equation model, J. Comput. Phys., 229:2773-2809, 2010.

[3] V. Faucher, S. Kokh, Extended Vofire algorithm for fast transient fluid-structure dynamics with liquid-gas flows and interfaces, J. Fluids and Struc., accepted January 2013.

[4] V. Faucher, projet RePDyn, Rapport d'avancement T0+18, http://www.repdyn.fr, 2011.

[5] computations on cluster of multi-processor, PASCO'07 (2007)