Formulation continue en vitesse du contact frottement.
Application à la simulation du temps de chute d'une
grappe de commande dans un cœur de réacteur nucléaire

Marc Kham, Mohamed Torkhani, Patrick Massin, Hachmi Ben Dhia

► To cite this version:

Marc Kham, Mohamed Torkhani, Patrick Massin, Hachmi Ben Dhia. Formulation continue en vitesse du contact frottement. Application à la simulation du temps de chute d'une grappe de commande dans un cœur de réacteur nucléaire. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422173

HAL Id: hal-01422173

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

Submitted on 23 Dec 2016

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
# Formulation continue en vitesse du contact frottement. Application à la simulation du temps de chute d'une grappe de commande dans un cœur de réacteur nucléaire.

M. Kham¹, M. Torkhani¹, P. Massin², H. Ben Dhia³

¹Électricité de France
1, avenue du Général de Gaulle, 92141 Clamart Cedex, France
{marc.kham, mohamed.torkhani}@edf.fr

²LAboratoire de Mécanique des Structures Industrielles Durables
UMR EDF-CNRS 2832
1, avenue du Général de Gaulle, 92141 Clamart Cedex, France
patrick.massin@edf.fr

³École Centrale Paris, MSSMAT
Grande Voie des Vignes, 92295 Châtenay-Malabry Cedex, France
bendhia@mssmat.ecp.fr

**Résumé** – On propose une formulation continue augmentée en vitesse [4] du contact-frottement pour la dynamique. On montre une application à la modélisation de la chute de cette grappe de commande dans le système de guidage [1]. On montre une très bonne correspondance avec les résultats expérimentaux aussi bien pour les évaluations des forces d'insertion que pour le temps de chute. On montre aussi un progrès notoire réalisé par rapport à une formulation Lagrangienne ponctuelle du contact.

**Mots clés** – contact, frottement, formulation continue augmentée, dynamique, theta-méthode.

## 1 Introduction

Les grappes de commande servent à piloter la puissance des réacteurs des centrales REP, et surtout à les arrêter en cas de besoin par chute rapide sous poids propre. Le temps de chute maximal autorisé pour la sécurité de l'arrêt doit impérativement être respecté sous peine de non-redémarrage de la tranche. Ces grappes traversent successivement un système de mécanisme de levée (mécanisme de commande), un guide de grappe au niveau du plénum supérieur du cœur (guidage discontinu) et les tubes-guide de l'assemblage combustible (guidage continu) (FIG.1).

En début de vie, un assemblage combustible est pratiquement droit, et le frottement entre la grappe et le guide peut être négligé. Seule la force hydraulique s'oppose alors à la chute de la grappe. Mais au cours de sa vie, l'assemblage combustible subit des chargements (mécaniques ou thermiques) extrêmement violents, qui peuvent entraîner une

<figure><img src="image_2.png" /><figcaption>**Figure 1-** Grappe de commande dans son assemblage combustible.</figcaption></figure>
La maîtrise du temps de chute de la grappe sous poids propre est un enjeu important pour EDF. A cette fin, il est crucial d'avoir une modélisation fine des phénomènes de contact/frottement se produisant entre la grappe et le guidage au cours de la chute. Dans cette optique se posent deux axes de recherche : ① un modèle performant de contact/frottement ; ② son association à un schéma d'intégration temporelle conservatif. Dans cet article, on propose une formulation Lagrangienne augmentée en vitesse du contact associée à un θ-schéma [3, 5, 6], ce que l'on décrira brièvement dans le paragraphe §2. Dans le paragraphe §3, en considérant un modèle simple de rotation de bloc, on cherchera à discuter l'intérêt de la méthode par rapport à des méthodes plus classiques. Dans le paragraphe §4, on présentera les résultats de simulation pour la chute d'une grappe de commande. Enfin, on terminera au paragraphe §5 par les perspectives de développement et d'amélioration en cours de la méthode.

# 2 Méthode d'intégration temporelle du contact-frottement

## 2.1 Formulation augmentée en vitesse du contact

Les lois de contact en vitesse sont connues sous le nom de Signorini-Moreau. Elles s'expriment comme suit (pour tout point x de la surface de contact) :

$$ 
\left\{
\begin{array}{ll}
\text{si} & d_n < 0 \quad \Rightarrow \quad \lambda = 0 \\
\text{sinon} & [[v_n]] < 0; \lambda \le 0 \text{ et } [[v_n]] \times \lambda = 0
\end{array}
\right.
\qquad (1)
$$

où $d_n$ est le jeu normal au point x entre les deux surfaces de contact, $[[v_n]]$ le saut de vitesse normale correspondant et $\lambda$ le multiplicateur de contact (pression) associé.

Ces lois traduisent un contrôle du contact, à la fois au niveau des déplacements et des vitesses. En introduisant deux champs de signe $S_u$ et $S_v$, les inéquations de Signorini-Moreau (1) peuvent être réécrites sous forme augmentée :

$$ 
\left\{
\begin{array}{l}
\lambda = S_u S_v (\lambda - \rho_n . [[v_n]]) \\
S_u = \chi (-d_n) \\
S_v = \chi (\lambda - \rho_n . [[v_n]])
\end{array}
\right.
\qquad (2)
$$

où $\chi$ est la fonction indicatrice de $\mathbb{R}^n$ définie par $\chi(x) = \begin{cases} 1 & x \le 0 \\ 0 & x > 0 \end{cases}$ et $\rho_n$ un réel strictement positif.

## 2.2 Formulation augmentée du frottement

Pour représenter les phénomènes de frottement, on adopte la loi de Coulomb qui s'écrit, pour tout x de la surface de contact :

$$ 
\left\{
\begin{array}{l}
\|\mathbf{r}_{\tau}(\mathbf{x}, t)\| < \mu |\lambda(\mathbf{x}, t)| \quad \Leftrightarrow \quad \mathbf{v}_{\tau}(\mathbf{x}, t) = 0 \\
\|\mathbf{r}_{\tau}(\mathbf{x}, t)\| = \mu |\lambda(\mathbf{x}, t)| \quad \Leftrightarrow \quad \exists \alpha \in {}^{*+} > 0 / \mathbf{v}_{\tau}(\mathbf{x}, t) = \alpha \mathbf{r}_{\tau}(\mathbf{x}, t)
\end{array}
\right.
\qquad (3)
$$

où $\mu$ est appelé coefficient de frottement de Coulomb, $\mathbf{v}_\tau$ vitesse relative tangente et $\mathbf{r}_\tau$ le multiplicateur (ou pression) de frottement associé.

On montre que la loi de Coulomb (3) peut être réécrite de manière augmentée [6] :

$$ 
\begin{array}{l}
\mathbf{r}_{\tau} = \mu \lambda \mathbf{\Lambda} \\
\mathbf{\Lambda} - \mathbf{P}_{B(0,1)}(\mathbf{g}_{\tau}) = \mathbf{0} \\
\mathbf{g}_{\tau} = \mathbf{\Lambda} + \rho_{\tau} \mathbf{v}_{\tau}
\end{array}
\qquad (4)
$$

où $\rho_\tau$ est un réel strictement positif, homogène à l'inverse d'une vitesse, et $\mathbf{P}_{B(0,1)}$ l'opérateur
de projection sur la boule unité. $\Lambda$ est un vecteur appelé semi-multiplicateur de frottement, et $\mathbf{g}_\tau$ un vecteur appelé semi-multiplicateur de frottement augmenté.

La loi de Coulomb pour le frottement (4) est complétée par une équation de type exclusion :

$$ (1 - \chi(g_n)) \times \Lambda = 0 \qquad (5) $$

Après manipulation des équations 4 et 5, on arrive à cette forme ultime des conditions de frottement de Coulomb :

$$ \frac{\mu \lambda \Delta t}{\rho_{\tau}} \chi(g_n) \left\{ \Lambda - \mathbf{P}_{B(0,1)} \left( \Lambda + \rho_{\tau} \mathbf{v}_{\tau} \right) \right\} + (1 - \chi(g_n)) \Lambda = 0 \qquad (6) $$

## 2.3 Formulation en vitesse des équations de contact-frottement

Le principe des puissances virtuelles appliquées à deux solides hyperélastiques 1 et 2 pouvant éventuellement entrer en contact-frottement, s'écrit [3, 4] :

Trouver les champs $(\mathbf{u}_k^1, \mathbf{u}_k^2, \mathbf{v}_k^1, \mathbf{v}_k^2, \lambda_k) \in \mathbf{CA}_v^1 \times \mathbf{CA}_v^2 \times \mathbf{CA}_v^1 \times \mathbf{CA}_v^2 \times H$ tels que :
$\forall (\mathbf{u}_k^{*1}, \mathbf{u}_k^{*2}, \mathbf{v}_k^{*1}, \mathbf{v}_k^{*2}, \lambda_k^{*}) \in \mathbf{CA}_v^1 \times \mathbf{CA}_v^2 \times \mathbf{CA}_v^1 \times \mathbf{CA}_v^2 \times H$

$$ \sum_{i=1}^{2} \int_{\Gamma_c} \rho_i^x \frac{\partial \mathbf{v}^i}{\partial t} \cdot \mathbf{v}^{*i} d\Omega + \sum_{i=1}^{2} \int_{\Gamma_c} \mathbf{F}^i \cdot \mathbf{S}^i : \nabla_x (\mathbf{v}^{*i}) d\Omega \\ - \int_{\Gamma_c} \mathbf{S}_u \cdot \mathbf{S}_v g_n [[\mathbf{v}_n^*]] d\Gamma_c - \int_{\Gamma_c} \mu \mathbf{S}_u \cdot \mathbf{S}_v \lambda \mathbf{P}_{B(0,1)}(\mathbf{g}_\tau) [[\mathbf{v}_\tau^*]] d\Gamma_c = 0 $$

$$ \int_{\Gamma_c} \frac{-1}{\rho_n} \{\lambda - \mathbf{S}_u \cdot \mathbf{S}_v g_n\} \lambda^* d\Gamma_c = 0 $$

$$ \int_{\Gamma_c} \frac{-\mu \mathbf{S}_u \cdot \mathbf{S}_v \lambda}{\rho_{\tau}} \{\Lambda - \mathbf{P}_{B(0,1)}(\mathbf{g}_{\tau})\} \Lambda^* d\Gamma_c + \int_{\Gamma_c} (1 - \mathbf{S}_u \cdot \mathbf{S}_v) \Lambda \cdot \Lambda^* d\Gamma_c = 0 $$

$$ g_n = \lambda - \rho_n [[\mathbf{v}_n]] $$

$$ \mathbf{g}_{\tau} = \Lambda + \rho_{\tau} (\mathbf{1} - \mathbf{n} \times \mathbf{n}) \cdot \Delta[[\mathbf{v}]] $$

$$ \mathbf{S}_u = \chi(-d_n) $$

$$ \mathbf{S}_v = \chi(\lambda - \rho_n [[\mathbf{v}_n]]) $$

$$ \mathbf{u}^i(t) = \mathbf{u}^i(t_0) + \int_{t_0}^{t} \mathbf{v}^i(\tau)d\tau \qquad (7) $$

où $\mathbf{CA}_v^i$ désigne l'espace des champs de vitesse cinématiquement admissibles pour le solide i. Cette formulation du contact sera dite *continue*, par opposition aux méthodes ponctuelles traditionnellement utilisées dans les codes de calcul.

Cette équation est résolue par un $\theta$-schéma d'ordre 1, défini par :

$$ \mathbf{u}^{+} = \mathbf{u}^{-} + \Delta t(1 - \theta)\mathbf{v}^{-} + \Delta t\theta\mathbf{v}^{+} \qquad (8) $$

$$ \mathbf{a}^{+} = \frac{\mathbf{v}^{+} - \mathbf{v}^{-}}{\Delta t} $$
## 3 Analyse du balancement d'un bloc

On se propose dans ce paragraphe de comparer la formulation *continue en vitesse* du contact à d'autres méthodes plus classiques. Pour cela on considère le cas du balancement d'un bloc rectangulaire supposé *rigide* reposant sur un demi-espace également supposé *rigide* (FIG.2). Ses dimensions sont 36×80 cm pour une masse de 417.6 Kg. Le bloc est légèrement incliné d'un angle initial de 10⁻² degrés, puis relâché. Soumis uniquement à la pesanteur, on suppose que le mouvement du bloc consiste en une série de rotations autour de ses deux coins (modèle dit de Housner [7]).

<figure><img src="image_3.png" /><figcaption>**Figure 2-** modèle de balancement de bloc : le bloc est soulevé d'un angle initial de 10⁻² degrés, puis relâché.</figcaption></figure>

On compare la formulation *continue en vitesse* du contact décrite précédemment (VELOCITY), à une formulation avec des *ressorts discrets* attachés aux coins du bloc (DISCRETE), à une formulation par *multiplicateurs de Lagrange* affectés au coins du bloc (LAGRANGE) et à une formulation *continue classique* contraignant uniquement le jeu en déplacement (GAP). Sur la Figure 3, on représente les solutions calculées à l'aide du code d'éléments finis d'EDF (Code_Aster) en termes d'évolution de l'énergie cinétique du bloc au cours du temps. Les valeurs analytiques des maxima attendus sont représentées par des points. Il apparaît que la solution donnée par la formulation *continue en vitesse* est la plus satisfaisante.

<figure><img src="image_4.png" /><figcaption>**Figure 3 :** modèle de balancement de bloc : le bloc est soulevé d'un angle initial de 10⁻² degrés, puis relâché.</figcaption></figure>
# 4 Modélisation du temps de chute d'une grappe de commande

La grappe de commande avec ses 24 crayons est modélisée par une *poutre équivalente*. De la même façon, les éléments du guidage (FIG.1) constituant l'assemblage combustible sont représentés par un alignement d'éléments *poutres* pour dessiner un couloir de guidage de la grappe (FIG.4) [1]. La partie inférieure du guidage, formant le tube-guide de l'assemblage combustible, est modélisée dans une *configuration déformée* typique (en forme de S) relevée sur les assemblages usagés.

Sur la Figure 5, on présente les résultats de simulations réalisées avec Code\_Aster, en comparant la formulation *continue en vitesse* aux formulations *Lagrangienne discrète* et *continue en déplacement*. Les solutions numériques sont ainsi comparées à la courbe expérimentale de temps de chute (évolution de la vitesse de chute en fonction du temps).

**Figure 4:** Modèle simplifié de la grappe de commande et des éléments de guidage de l'assemblage combustible.

Les différentes méthodes de contact permettent d'estimer de manière relativement satisfaisante à la fois l'évolution de la vitesse de chute et le temps de chute, égal à environ 1.52s. En particulier, les formulations continues en vitesse et en déplacement sont très proches l'une de l'autre, et coïncident parfaitement avec la courbe expérimentale. Comparativement, la formulation *Lagrangienne discrète*, moins « élaborée » que les précédentes, se révèle ici moins précise.

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

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

**Figure 5:** Évolution de la vitesse verticale de chute du crayon en fonction du temps : comparaison des formulations continues en déplacement (bleu), en vitesse (rouge) du contact avec une formulation *Lagrangienne ponctuelle* du contact (rose), et avec des résultats d'essai (noir).
# 5 Conclusions et perspectives

Dans cet article, nous avons présenté une formulation augmentée en vitesse du contact associée à un schéma d'intégration d'ordre 1 (θ-schéma). Cette formulation permet un contrôle de la vitesse de contact et le schéma d'ordre 1 d'assurer une meilleure stabilité de l'intégration temporelle. Pour tester l'apport d'une telle formulation, on a dans un premier temps considérer le cas-test d'un balancement de bloc, pour lequel la supériorité de la méthode introduite est mise en évidence. Enfin, on a réalisé une simulation du temps de chute d'une grappe de commande dans un assemblage combustible en configuration déformée, et montré une meilleure performance de cette méthode par rapport à un modèle de contact Lagrangien discret.

Mais la formulation augmentée du contact reste aujourd'hui relativement coûteuse : ainsi, la simulation (par Code_Aster) du temps de chute (§4) requiert 100 heures de calculs, contre environ 20 heures quand on utilise un modèle Lagrangien discret pour le contact. Cet inconvénient reste problématique pour les besoins d'industrialisation d'EDF. L'objectif principal porte donc sur l'amélioration de la performance de la méthode dans Code_Aster, avec quelques pistes déjà identifiées : résorption de la boucle de point-fixe sur la géométrie (grands déplacements) [4], utilisation d'un solveur itératif (Gradient Conjugué) au lieu d'un solveur direct etc...

# Références

[1] H. Andriambobolona, Bosselut D., Massin P. *Methodology for a numerical simulation of an insertion or a drop of the rod cluster control assembly in a PWR, Nuclear Engineering and Design*, Vol. **237**, 600-606, 2007.

[2] H. Ben Dhia, Zarroug M. *Hybrid frictional contact particles-in elements. Revue Européenne des Éléments finis*, Vol. **11**, n° 2-3-4, pp. 417-430, 2002.

[3] H. Ben Dhia *Level -sets fields, placement and velocity-based formulation of contact-impact problems*, Int. J. Numer. Methods Eng., Vol. **69**, 2711-2735, 2007.

[4] H. Ben Dhia, Kham M., Massin P., Torkhani M. *Un cadre général pour le traitement des problèmes d'interface*, Giens, 2009.

[5] C. Zammali *Contribution à la modélisation mécanique et numérique des problèmes de contact-impact*, thèse de doctorat École Centrale Paris, 2005.

[6] M. Kham *Éléments de contact dérivés d'une formulation hybride continue*, document de référence Code_Aster, R5.03.52, 2007.

[7] G.W. Housner *The behavior of inverted pendulum structures during earthquakes*. Bull. of the Seismological Society of America, **53**(2), pp.403-417, 1963.