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

HAL
open science

# From domain decomposition to model reduction for Large nonlinear structures

Bertrand Leturcq, Patrick Le Tallec

► To cite this version:

Bertrand Leturcq, Patrick Le Tallec. From domain decomposition to model reduction for Large nonlinear structures. Comptes Rendus. Mécanique, 2023, 351 (S1), pp.1-17. 10.5802/crmeca.168 . cea-04149398

HAL Id: cea-04149398

https://cea.hal.science/cea-04149398

Submitted on 3 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Mécanique

Bertrand Leturcq and Patrick Le Tallec

**From Domain Decomposition to Model Reduction for Large Nonlinear
Structures**

Published online: 15 May 2023

https://doi.org/10.5802/crmeca.168

**Part of Special Issue:** The scientific legacy of Roland Glowinski

**Guest editors:** Gregoire Allaire (CMAP, Ecole Polytechnique, Institut Polytechnique de Paris, Palaiseau, France), Jean-Michel Coron (Laboratoire Jacques-Louis Lions, Sorbonne Université) and Vivette Girault (Laboratoire Jacques-Louis Lions, Sorbonne Université)

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

BY

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
<figure><img src="image_4.png" /></figure>

The scientific legacy of Roland Glowinski / L'héritage scientifique de Roland Glowinski

# From Domain Decomposition to Model Reduction for Large Nonlinear Structures

## De la décomposition de domaine à la réduction de modèle pour les grandes structures non linéaires

Bertrand Leturcq$^a$ and Patrick Le Tallec$^b$

$^a$ Université Paris-Saclay, CEA, Service d'études mécaniques et thermiques, 91191 Gif-sur-Yvette, France

$^b$ M3DISIM/LMS, École Polytechnique, CNRS, Institut Polytechnique de Paris, 91128 Palaiseau, France

E-mails: bertrand.leturcq@cea.fr (B.Leturcq), patrick.letallec@polytechnique.edu (P.LeTallec)

In memoriam of Roland Glowinski (1937-2022), great scientist, leader and friend.

**Abstract.** The numerical simulation of multiscale and multiphysics problems requires efficient tools for spatial localization and model reduction. A general strategy combining Domain Decomposition and Nonuniform Transformation Field Analysis (NTFA) is proposed herein for the simulation of nuclear fuel assemblies at the scale of a full nuclear reactor. The model at subdomain level solves the full elastic problem but with a reduced nonlinear loading, based on simplified boundary conditions, reduced creep flow rules, projected sign preserving contact conditions, and a NTFA like reduced friction law to get the evolution of each slipping mode. With this loading reduction, the local solution can be explicitly obtained from a small set of precomputed elementary elastic solutions. The numerical tests indicate that considerable cost reduction (a factor of 50 to 1000) can be achieved while preserving engineering accuracy.

**Résumé.** La simulation numérique de problèmes multi-échelles et multi-physiques nécessite des outils performants de localisation spatiale et de réduction de modèle. Une stratégie générale combinant la décomposition de domaine et l'analyse de champ de transformation non uniforme (NTFA) est proposée ici pour la simulation d'assemblages de combustible nucléaire à l'échelle d'un réacteur nucléaire complet. Le modèle au niveau du sous-domaine résout le problème élastique complet mais avec un chargement non linéaire réduit, basé sur des conditions aux limites simplifiées, des règles d'écoulement de fluage réduites, des conditions de contact préservant le signe projeté et une loi de frottement réduit de type NTFA pour obtenir l'évolution de chaque mode de glissement. Avec cette réduction de chargement, la solution locale peut être obtenue explicitement à partir d'un petit ensemble de solutions élastiques élémentaires précalculées. Les tests numériques indiquent qu'une réduction considérable des coûts (un facteur de 50 à 1000) peut être obtenue tout en préservant la précision du calcul.

**Keywords.** Model reduction, Domain decomposition, Nuclear fuel assemblies, creep, contact.

**Mots-clés.** Réduction de modèle, Décomposition de domaines, Assemblage de combustible, fluage, contact.

**Funding.** This work has been partly finalized during a stay of the second author at the Oden Institute of University of Texas at Austin.
# 1. Introduction

Multiscale and multiphysics simulation becomes a necessity. On one hand, there is often a strong coupling between the global response of a system of interest and local dimensioning nonlinear phenomena such as nonlinear creep, multiscale friction, change of phases, molecules unfolding, liquid crystals microstructure evolution and so on. On the other hand, large amount of experimental data originating from different sources and scales must be used in the construction, validation and optimisation of the different components. Such issues occur for example in large scale industrial systems as encountered in nuclear industry [1], in the development and use of microstructured soft materials or architectures in transport industry, or in the development of patient specific medical models [2].

These problems have been central in the career of Roland Glowinski and his approach was based on basic fundamentals which proved to be successful in many important applications. He has introduced powerful localisation techniques based on operator [3] or on domain splitting [4], he has outlined the importance of properly specifying and controlling the boundary and interface conditions [5], and of respecting and exploiting the constitutive non linearities [6, 7].

Following these principles, the simulation of real life multiscale and multiphysics problems combines two basic ingredients. The first one concerns the data transfer strategy between the global model and the local problems. The notion of imposed averaged deformation as proposed in standard homogeneization techniques is not adequate when local models deal with molecular or fibril networks as in [8] or in [2]. Additional work is therefore needed to better understand and expand the microsphere models which are used in those directions. Averaged deformation is also not sufficient when the scales are not well separated. Domain decomposition based strategies can then be a possible solution, especially in situations like in [9, 10] where the construction of subdomains is very straightforward.

The other issue deals with the necessary model reduction which must be performed at the local scale. An accurate treatment of all nonlinearities and multiphysics microscale phenomena occurring at that local scale is often out of reach. The Nonuniform Transformation Field Analysis (NTFA) strategy of [11] proposes an approximation of the local nonlinear fields by projecting them on an orthogonal basis of precomputed modes. The challenge is then to project the nonlinear constitutive law on these modes without using exact local integration at each time step, which would be of the same complexity as the original multiscale problem. Extending this approach to more complex situations also requires to build sign preserving modal expansion of contact forces or of other positive internal variables.

The particular problem on which we would like to illustrate these fundamentals is the study of the progressive deformation of nuclear fuel assemblies that occurs during a succession of irradiation cycles within a nuclear reactor [12, 13]. The engineering challenge is to limit the fuel assemblies permanent bows which can be at the origin of costly incidents during the pull out or of poor functioning of the reactor's control safety bars. The collective response of the hundreds of fuel assemblies inside the core couples the global effect of the hydraulic loading, the localized creep of the structural elements and the local evolution of the friction force inside each holding spring. This multi-body multi-physics calculations involves something like 4 million contact surfaces and up to sixty thousands fuel rods. The time scale is of the order of one year and corresponds to a reactor loading cycle. The spatial scale ranges from a rod length of 5m to spring fixing sizes of 1mm. A full model would require more than 300 millions degrees of freedom and a considerable effort in the local integration of creep, contact and friction. Model reduction is thus needed. In order to improve the representativeness of the core model without penalizing
the computational cost of the coupled simulations, this model reduction is proposed herein at the scale of each of the 241 fuel assemblies present inside the reactor core. It borrows the concepts of both domain decomposition [14] and of non-linear homogenization as proposed in [11].

The first idea herein is thus to consider a single fuel assembly as the Representative Volume Element (RVE) classically used in multiphase homogenization. In this framework, the solution within a RVE in presence of creep and friction is historically obtained by Transformation Field Analysis (TFA) as introduced by [15] and [16]. In their construction, these authors have selected a limited number of phases within the RVE, each nonlinear phase being supposed to be uniform in space. This crude approximation allows to apply the original constitutive law (creep or friction) to predict the evolution of the constitutive variable inside the phase p at the scale of each nonlinear phase. For a RVE with elasticity tensor C, given history of irreversible strain $\epsilon_{=p}^{in}$ and subjected to a given macro deformation $\tilde{\epsilon}$ imposed in average to the RVE, the resulting stress field to be used in equilibrium equations can be explicitly given at all points by

$$ \underline{\underline{\sigma}} = \underline{\underline{C}} : \underline{\underline{A}} : \underline{\underline{\tilde{\epsilon}}} + \sum_{p} \underline{\underline{C}} : (\underline{\underline{D}}_{p} - \underline{\underline{I}}_{p}) \underline{\underline{\epsilon}}_{=p}^{in} $$

with localisation and influence tensors A and Dp obtained a priori by solving elementary elasticity problems at the scale of the RVE. The NTFA strategy of [11, 17, 18] proposes a more accurate approximation of the inelastic strain fields $\underline{\underline{\epsilon}}_{=p}^{in}$ by replacing the constant fields by a linear expansion on an orthogonal basis of strain modes obtained at a controlled accuracy by Proper Orthogonal Decomposition (POD) [19] based on a collection of precomputed snapshots. A simple but efficient averaging rule is also used for projecting nonlinear viscoplastic constitutive laws on to the selected modes.

Extending this approach to nonlinear industrial structures requires additional ingredients in order to extend the notion of imposed deformation to a complex substructure and to build a sign preserving modal expansion of the contact forces. For the first point, we use a domain decomposition strategy. A substructure like a fuel assembly is a natural subdomain of the full structure, and its motion can be accurately controlled by the rigid body motions of each of its interface with the global structure. On the other hand, the modal expansion of the contact forces uses positive modes obtained through a Nonegative Matrix Factorisation (NMF) of contact forces snapshots as proposed in [20]. The local contact conditions are then transferred to the reduced model through a simple projection of the contact constraints on these positive modes. Last, a modal basis of slip displacements is obtained by POD, to be combined to averaged stick slip criteria to be directly written in terms of slipping modes.

Altogether, the reduced model solves the full elastic problem but with a reduced nonlinear loading, based on simplified boundary conditions, reduced creep laws (as in standard NTFA) to get the evolution of the creep modes, reduced (NMF projected) contact conditions to get the reduced normal forces, and a NTFA like reduced friction law to get the evolution of each slipping mode. With this loading reduction, the local solution can be explicitly obtained from a small set of precomputed elementary elastic solutions.

The paper is organized as follows. The full problem is described in Section 2, where it is decomposed in local subproblems, with adequate reduction of the local nonlinear fields using POD like strategies. The formalism for the construction of reduced modal constitutive laws and its application to nonlinear creep, contact and friction are presented in Section 3. Numerical results are presented in Section 4, before drawing conclusions in Section 5.
## 2. Description and decomposition of the full model

### 2.1. *The mechanical problem*

We consider the long term mechanical evolution of fuel assemblies within the core of a pressurized water nuclear reactor. A typical reactor contains more than 200 fuel assemblies that will bow and deform under the combined action of hydraulic forces, thermal loads, irradiation growth and irradiation creep.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** The pressure vessel of a nuclear reactor core, containing approximately 240 fuel assemblies, taken from [21].</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** The loading environment of the fuel assemblies inside a reactor core [13].</figcaption></figure>

A typical fuel assembly is 5 *m* high, 0.21 *m* wide, contents more than 250 thin fuel rods which are maintained on 10 spacing grids assembled together through 25 guide tubes. The fixings of
the fuel rods onto the spacing grids is achieved by small holding springs working in compression and friction, involving in each fuel assembly 17 000 contact surfaces between the rods and the holding springs, at a cm scale. This holding system controls the lateral motion of the fuel rods while allowing for their irradiation growth. It is also at the origin of a global hysteretic response of the structure.

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

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Bowed fuel assembly [22]. A typical map of fuel assemblies bows inside a nuclear core as reported in [22].</figcaption></figure>

In this framework, it is important to monitor and limit the permanent bows of fuel assemblies. The internal forces which are at the origin of permanent deformation of this otherwise elastic structure include non uniform irradiation growth, irradiation creep, trapped friction forces and local relaxation of the holding springs.

The existence of irradiation creep transforms the elastic constitutive law into

$$ \underline{\underline{\sigma}}(x) = \mathbf{C}(x) : (\underline{\underline{\varepsilon}}(x) - \underline{\underline{\varepsilon}}^{in}(x)) \quad (1) $$

while the evolution of the creep strain is governed by a flow rule of the type

$$ \underline{\underline{\dot{\varepsilon}}}^{in} = \frac{3}{2} \frac{\partial \psi(\sigma_{eq})}{\partial \sigma_{eq}} \frac{\underline{\underline{\sigma}}'}{\sigma_{eq}}. \quad (2) $$

Above **C** denotes the local elasticity tensor, $\sigma_{eq} = \sqrt{\frac{3}{2}\underline{\underline{\sigma}}' : \underline{\underline{\sigma}}'}$ the local equivalent Von Mises stress, $\underline{\underline{\sigma}}'$ the local deviatoric stress, and $\psi(\sigma_{eq})$ the local viscoplastic potential. A typical potential properly describing the creep flow of metals under irradiation is

$$ \psi(T, \sigma_{eq}) = \frac{1}{\tau(x)} \sigma_0(x) \left( \frac{\sigma_{eq}}{\sigma_0(x)} \right)^m \quad (3) $$

with a material dependent exponent $m$ close to two in many cases, and characteristic creep time $\tau$ and yielding stress $\sigma_0$ functions of the local temperature and of the local amount of cumulated irradiation.

As for contact, each holding spring $J$ acts in compression with a negative contact force given by

$$ \tau_N(J) = k_N(J) \min(0, \delta u^N(J) - \varepsilon^N(J)). \quad (4) $$
<figure><img src="image_9.png" /><figcaption>Figure 4. Mechanical description of a fuel assembly [23].</figcaption></figure>

Above $k_N(J)$ denotes the stiffness of the holding spring, $\delta u^N(J)$ the difference of normal displacement between the spring extremities, and $\varepsilon^N(J)$ denotes the positive or negative prestrain of the holding spring which can vary in time due to spring relaxation under irradiation. Introducing the local gap defined by

$$ \text{gap}(J) = \delta u^N(J) - \varepsilon^N(J) - \tau_N(J)/k_N(J), \quad (5) $$

this constitutive law rewrites as a set of complementary inequalities

$$ \text{gap}(J) \ge 0, \quad -\tau_N(J) \ge 0, \quad -\text{gap}(J)\tau_N(J) = 0. \quad (6) $$

As for the friction force at a given fixation point, it is supposed to be governed by a regularized Coulomb friction law relating this friction force to the local sliding velocity $\dot{u}_t$ by

$$ \underline{\tau}_t = -\min \left( \mu_t, \mu \frac{|\tau_N|}{\|\dot{u}_t\|} \right) \dot{u}_t. \quad (7) $$

## 2.2. Subdomain decomposition and loading reduction

A fuel assembly is a natural subdomain of the full structure. Moreover, we know from domain decomposition studies [14] that a subdomain motion can be accurately controlled by the rigid body motions of each of its interface with an error of at most $O(r^2/H^2)$ in energy norm, with $r$ the diameter of the interface, and $H$ the distance appearing when estimating the norm of the lifting operator building an equilibrium solution from an imposed trace, distance which will be here the distance between two consecutive interfaces. The mechanical interfaces of a fuel assembly are the spacing grids, which are subjected to the hydraulic forces applied by the local flow of the pressurized water and to potential contact forces with neighboring assemblies. Conversely these loads are mostly controlled by the average position of the spacing grids. In this framework, the “global deformation” of a fuel assembly is thus defined by the three translational modes of each
of its 11 spacing grids. Its evolution is then characterized by its elastic response to the imposed external and nonlinear internal forces. In view of the large dimension of this loading space, a model reduction is then needed at that level. To account for hydraulic forces, nonlinear creep, contact and friction, we need to introduce in the spirit of NTFA [11]

(1) the field of external forces $\{f_{jg}^{ext}\}_{jg}$ applied in average on each spacing grid $g$ and along each direction $j$;

(2) an orthonormal basis $\{\varepsilon_{i}^{in}\}_{i=1,N_i}$ of small dimension $N_i$ in the space of irreversible creep strain fields on which to develop the local creep strain fields

$$ \underline{\varepsilon}^{in}(x) = \sum_{i} \alpha_{i} \underline{\varepsilon}_{=i}^{in}(x) $$

at each material point $x$ of the structure. This basis is to be constructed by Proper Orthogonal Decomposition [19] on a large collection of snapshots of creep strain fields $\underline{\varepsilon}^{in}(x)$ obtained by performing series of direct simulations of the local problems;

(3) an orthonormal basis $\{\underline{u}_{m}^{s}\}_{m=1,N_m}$ of small dimension $N_m$ on which to develop the slipping velocities at all fixation points $J$ into $\underline{u}_{t}(J) = \sum_{m} \dot{\alpha}_{m} \underline{u}_{m}^{s}(J)$, basis to be also constructed by Proper Orthogonal Decomposition on a large collection of snapshots of slipping velocity fields $\underline{u}_{t}(J)$ obtained in the same direct simulations of the local problems;

(4) a positive decomposition of a matrix of $P$ contact forces snapshots $X$ on $N_k$ positive modes using a Nonnegative Matrix Factorisation $X \approx \hat{\Lambda}H$ obtained by least squares minimization of the approximation error $\|X - \hat{\Lambda}H\|_{2}^{2}$ using alternating directions [24,25]. The result is a modal matrix $\hat{\Lambda}$ of $N_k$ positive columns $\Lambda_k$ (typically $P \gg N_k \approx 10$) on which to decompose the local contact forces $\tau_N(J) \approx -\sum_k f_k \Lambda_k(J)$.

The strain modes $\underline{\varepsilon}_{=j}^{in}(x)$ are unit free, the slipping modes $\underline{u}_{m}^{s}(J)$ will be properly normalized in order to be adimensional, and the local values of the contact modes $\Lambda_k(J)$ are expressed in Newtons. With such a modal decomposition of the external and internal forces, the subdomain problem can be explicitly solved at low cost by solving the weak form of the equilibrium equations tested against any kinematically admissible test function $\mathbf{U}^* = \underline{u}^*(x) \in \mathbf{V}$. In our framework, the space $\mathbf{V}$ of kinematically admissible displacement corresponds to displacement fields which are zero on the part of boundary where the assembly is potentially clamped and which present zero slipping. On the other hand, such fields allow for arbitrary normal overlaps at contact points. They are thus “free contact modes”. The local equations of equilibrium write then

$$ \langle \mathbf{KU}, \mathbf{U}^* \rangle - \sum_{i=1}^{N_i} \alpha_i \langle \mathbf{F}_i^{in}, \mathbf{U}^* \rangle + \sum_{k=1}^{N_k} f_k \langle \mathbf{F}_k^{N}, \mathbf{U}^* \rangle = \sum_{jg} f_{jg}^{ext} \langle \mathbf{F}_{jg}^{ext}, \mathbf{U}^* \rangle, \quad \forall \mathbf{U}^* \in \mathbf{V}. \quad (8) $$

The notation here are rather standard. The fuel assembly elastic stiffness matrix is given by

$$ \langle \mathbf{KU}, \mathbf{U}^* \rangle = \int_V \mathbf{C} \underline{\varepsilon}(\mathbf{U})(x) : \underline{\varepsilon}(\mathbf{U}^*) (x) dV, $$

the term $\langle \mathbf{F}_{jg}^{ext}, \mathbf{U}^* \rangle$ characterizes the virtual power developed by a unit force in the corresponding averaged translational grid motion

$$ \langle \mathbf{F}_{jg}^{ext}, \mathbf{U}^* \rangle := \frac{1}{S_g} \int_{grid-g} \hat{u}_j(x_{ext}) dS $$

with $S_g = \int_{grid-g} dS$ the section of the corresponding grid, the term $\langle \mathbf{F}_i^{in}, \mathbf{U}^* \rangle$ corresponds to the virtual power developed by the internal creep strains (in modal form)

$$ \langle \mathbf{F}_i^{in}, \mathbf{U}^* \rangle = \int_V \mathbf{C} \underline{\varepsilon}_{=i}^{in}(x) : \underline{\varepsilon}(\mathbf{U}^*) (x) dV, $$
the elementary contact loads $\langle \mathbf{F}_k^N, \mathbf{U}^* \rangle$ correspond to the virtual power developed by the contact mode $\Lambda_k$

$$ \langle \mathbf{F}_k^N, \mathbf{U}^* \rangle = \sum_J \Lambda_k(J) \delta(u^*)^N(J). $$

As for the unknown displacement field, for a known slipping motion $\sum_m \alpha_{sm} \mathbf{U}_m^s$, it is by construction given by $\mathbf{U} = \hat{\mathbf{U}} + \sum_m \alpha_{sm} \mathbf{U}_m^s$ with $\hat{\mathbf{U}}$ obtained by solving (8) and $\mathbf{U}_m^s$ the nodal representation of the displacement field obtained by zero extension of the tangential slip mode $\underline{u}_m^s$. Above V denotes either the domain occupied by one assembly, or its volume.

Observe that as usual in such frameworks the solution of (8) can be obtained as an explicit linear combination of $N_i + N_j + N_k + N_m$ precomputed elementary elastic solutions. In more details, we have

$$ \mathbf{U} = \sum_{Jg} \alpha_{Jg} \hat{\mathbf{U}}_{Jg}^{\text{ext}} + \sum_{i} \alpha_{i} \hat{\mathbf{U}}_{i}^{\text{in}} + \sum_{k} f_k \hat{\mathbf{U}}_{k}^{\text{N}} + \sum_{m} \alpha_{sm} (\hat{\mathbf{U}}_{m}^{s} + \mathbf{U}_{m}^{s}) \quad (9) $$

with $\alpha_{Jg} = f_{Jg}^{\text{ext}}$ and under the notation

$$ \hat{\mathbf{U}}_{Jg}^{\text{ext}} = \mathbf{K}^{-1} \mathbf{F}_{Jg}^{\text{ext}}, \quad \hat{\mathbf{U}}_{i}^{\text{in}} = \mathbf{K}^{-1} \mathbf{F}_{i}^{\text{in}}, \quad \hat{\mathbf{U}}_{ik}^{\text{N}} = \mathbf{K}^{-1} \mathbf{F}_{k}^{\text{N}}, \quad \hat{\mathbf{U}}_{m}^{s} = \mathbf{K}^{-1} \mathbf{F}_{m}^{s}. $$

The last term $\mathbf{F}_m^s$ in the above equilibrium equations corresponds to the elastic internal force induced by slipping

$$ \langle \mathbf{F}_m^s, \mathbf{U}^* \rangle = - \langle \mathbf{K} \mathbf{U}_m^s, \mathbf{U}^* \rangle. $$

The expression (9) translates into a similar decomposition of the deviatoric stress field $\underline{\sigma}'$ and of the normal overlap $\delta u^N$.

**Remark 1.** When dealing with a collection of fuel assemblies inside a reactor core, the forces $f_{Jg}^{\text{ext}}$ applied on the grids of the different assemblies include hydraulic forces and potential inter assembly contact forces which are given as a nonlinear function of the relative positions of neighboring grids in the spirit of (4) and the equation $\alpha_{Jg} = f_{Jg}^{\text{ext}}$ becomes an implicit nonlinear equation in $\{\alpha_{kh}\}_{kh}$

$$ \alpha_{Jg} = f_{Jg}^{\text{ext}}(\{\alpha_{kh}\}_{kh}). \quad (10) $$

# 3. Reduced nonlinear laws

## 3.1. Basic principle

The problem at local level will be specified if we can find a proper way to construct a reduced viscoplastic flow rule controlling the evolution of the modal creep coefficients $\alpha_i$ as function of their dissipation work (in Joules)

$$ \sigma_j := \int_V \underline{\sigma} := \underline{\epsilon}_j^{\text{in}} dV, \quad (11) $$

a reduced contact law relating the (adimensionalized) modal contact forces $f_k$ to the penetration work (expressed in Joules)

$$ \delta U_k^N = \sum_J \delta U^N(J) \Lambda_k(J), \quad (12) $$

and a reduced friction law $F_m^s(\{\dot{\alpha}_l\}_l)$ relating the virtual power of the friction forces $S_m^s$ in the slip mode $\mathbf{U}_m^s$ (measured in Watts) to the modal coefficients $\{\dot{\alpha}_m\}_m$ (expressed in meter per second) of the slipping velocity

$$ S_m^s := \sum_J \underline{\tau}_t(J) \cdot \underline{u}_m^s(J) = F_m^s(\{\dot{\alpha}_l\}_l), \quad (13) $$

where $\underline{\tau}_t(J)$ denotes the friction force applied on the $J^{th}$ contact. The problem is to compute at low cost the evolution of these reduced nonlinear loadings while respecting as much as possible the power dissipated by these loadings. This is achieved by Galerkin projection of the nonlinear constitutive laws onto the reduced space of internal loadings. For the contact law, this
projection must respect the sign of the contact forces. For creep and friction, this projection must be complemented by a reduced integration hence by a further approximation, otherwise the projection would have the same complexity as the full problem which would not be efficient.

## 3.2. Reduced contact law

The projection strategy to be used to construct a reduced contact law must respect the sign constraint on the projected contact force. In terms of energy and in presence of contact, the detailed model formally has a saddle point structure

$$ \inf_{\hat{u}(x) \in \mathbf{V}_{\hat{\tau}_{N}(\cdot)} \le 0} \sup_{J} \left\{ W(\hat{u}) + \sum_{J} \left( (\delta u^{N}(J) - \varepsilon^{N}(J)) \hat{\tau}_{N}(J) - \frac{\hat{\tau}_{N}^{2}(J)}{2k_{N}(J)} \right) \right\} $$

with $W(\hat{u})$ denoting the elastic and potential energy of the structure. Indeed the optimality conditions in terms of the normal reaction force $\tau_N$ take the form (4) while the equilibrium equations would correspond to the optimality conditions in displacement. The reduced contact law is then simply obtained as in [20] by reducing the set of admissible contact forces to the negative cone built from the selected modal forces

$$ \mathbf{V}_{\tau} := \left\{ \hat{\tau}_{N} = - \sum_{k} \hat{f}_{k} \Lambda_{k}, \hat{f}_{k} \ge 0, \forall k = 1, N_{k} \right\}. $$

The corresponding optimality conditions in $\{\hat{f}_k\}_k$ now formally become

$$ \sum_{J} \left( \frac{1}{k_N(J)} \left( \sum_{l} f_l \Lambda_l(J) \right) \Lambda_k(J) + (\delta u^N(J) - \varepsilon^N(J)) \Lambda_k(J) \right) - \text{gap}_k = 0, \forall k = 1, N_k, \quad (14) $$

where the Lagrange multipliers $\text{gap}_k$ of the positivity constraint $f_k \ge 0$ satisfy the complementary conditions

$$ \text{gap}_k \ge 0, \quad f_k \ge 0, \quad \text{gap}_k f_k = 0, \quad \forall k = 1, N_k. \quad (15) $$

Observe that, in difference with [20], we do not reduce the displacement space in the lagrangian writing of the contact problem, which guarantees that the classical Inf-Sup stability condition will be satisfied after reduction. The modal contact conditions keep the exact form of the local ones, within the projection of the local gaps $\text{gap}(J) = \delta u^N(J) - \varepsilon^N(J) - \tau_N(J)/k_N(J)$ onto the local modes $\{\Lambda_k\}_k$. As before, the sum $\sum_J \delta u^N(J) \Lambda_k(J)$ is an explicit function of the expansion coefficients $\alpha_{jg}$, $\alpha_i$, $f_k$ and $\alpha_{sm}$ in (9), namely

$$ \sum_{J} \delta u^{N}(J) \Lambda_{k}(J) = \sum_{jg} \alpha_{jg} \sum_{J} \delta^{N} \hat{\mathbf{U}}_{jg}^{\text{ext}}(J) \Lambda_{k}(J) + \sum_{i} \alpha_{i} \sum_{J} \delta^{N} \hat{\mathbf{U}}_{i}^{\text{in}}(J) \Lambda_{k}(J) \\ + \sum_{k} f_{k} \sum_{J} \delta^{N} \hat{\mathbf{U}}_{k}^{\text{N}}(J) \Lambda_{k}(J) + \sum_{m} \alpha_{sm} \sum_{J} (\delta^{N} \hat{\mathbf{U}}_{m}^{s} + \delta^{N} \mathbf{U}_{m}^{s}) (J) \Lambda_{k}(J) \quad (16) $$

where all lengthy summation on $J$ can be precomputed, which enables a fast solution of (14).

## 3.3. Reduced creep law

For the creep part, we need to find the evolution $\sum_i \dot{\alpha}_i \dot{\epsilon}^{in} \equiv \dot{\epsilon}^{in}$ of the projected creep strain which will dissipate the same energy as the one dissipated by the exact creep strain rate

$$ \dot{\epsilon}^{in} = \frac{3}{2} \frac{\partial \psi(\sigma_{eq})}{\partial \sigma_{eq}} \frac{\sigma'}{\sigma_{eq}} $$
in any stress field $\sigma_0 \underline{\underline{\epsilon}}_j^{in}$ defined on our reduced space of nonlinear viscoplastic loading. Above, we have kept the notation of (2) with $\underline{\underline{\sigma}}'$ denoting the local deviatoric stress $\underline{\underline{\sigma}}' = \underline{\underline{\sigma}} - \frac{1}{3}\text{Tr}(\underline{\underline{\sigma}})$ with norm $\underline{\underline{\sigma}}_{eq}$. The reduced flow rule should thus be solution of the projection problem

$$ 
\begin{aligned}
\int_V \sigma_0 \left( \sum_i \dot{\underline{\underline{\alpha}}}_i \underline{\underline{\epsilon}}_i^{in} \right) : \underline{\underline{\epsilon}}_j^{in} dV &= \int_V \sigma_0 \dot{\underline{\underline{\epsilon}}}_i^{in} : \underline{\underline{\epsilon}}_j^{in} dV, \\
&= \int_V \frac{3m}{2\tau(x)} \left( \frac{\sigma_{eq}^2(x)}{\sigma_0^2(x)} \right)^{m/2-1} \underline{\underline{\sigma}}'(x) : \underline{\underline{\epsilon}}_j^{in}(x) dV, \quad \forall j = 1, N_i.
\end{aligned} 
\quad (17) $$

The presence of the nonlinear term $(\frac{\sigma_{eq}^2(x)}{\sigma_0^2(x)})^{m/2-1}$ in the above integral makes it impossible to precompute the right hand side integral independently of the values of the modal coefficients of the internal loading. For situations with uniform coefficients $\tau$ and $\sigma_0$, the idea of [11] is to replace the local value of $\underline{\underline{\sigma}}_{eq}^2(x)$ by the volume average $s_r^2 = \langle \underline{\underline{\sigma}}_{eq}^2 \rangle$ of the square of the local Von Mises equivalent stress as computed after orthogonal projection onto the reduced space of viscoplastic loading. If the basis vectors $\underline{\underline{\epsilon}}_j^{in}$ are properly scaled in the sense that $\int_V \underline{\underline{\epsilon}}_i^{in} : \underline{\underline{\epsilon}}_j^{in} dV = V \delta_{ij}$, we would simply get

$$ s_r^2 = \frac{3}{2V^2} \sum_i \sigma_i^2 \quad \text{with} \quad \sigma_i = \int_V \underline{\underline{\sigma}}'(x) : \underline{\underline{\epsilon}}_i^{in}(x) dV $$

yielding as flow rule

$$ \dot{\underline{\underline{\alpha}}}_i = \frac{3m}{2\tau V} \left( \frac{s_r^2}{\sigma_0^2} \right)^{m/2-1} \frac{1}{\sigma_0} \underline{\underline{\sigma}}_i. $$

The function $s_r^2$ is a quadratic function of the stress field, hence from (9) it reduces to a quadratic form directly computed in terms of modes with no need of an element wise integration of the detailed model at each time step. This results in a considerable gain of computing time and turns out to be very efficient for mildly nonlinear creep.

The nonuniform case is more complex. The correct energy based scalar product is now defined by the mass matrix

$$ M_{ij} = \int_V \sigma_0(x) \underline{\underline{\epsilon}}_i^{in}(x) : \underline{\underline{\epsilon}}_j^{in}(x) dV \quad (18) $$

and the accessible projected quantity is the volume average $s_r^2 = \langle \underline{\underline{\sigma}}_{eq}^2 \rangle$ of the square of the local Von Mises equivalent stress after division by the yield stress $\sigma_0$ and orthogonal projection onto the reduced space of viscoplastic loading. Indeed the coefficients $\alpha'_i$ of the projection of $\frac{1}{\sigma_0} \underline{\underline{\sigma}}_i$ on to the reduced space for the above scalar product are solution of

$$ \int_V \sigma_0 \left( \sum_i \alpha'_i \underline{\underline{\epsilon}}_i^{in} \right) : \underline{\underline{\epsilon}}_j^{in} dV = \int_V \underline{\underline{\sigma}}'(x) : \underline{\underline{\epsilon}}_j^{in}(x) dV = \sigma_j, \forall j = 1, N_i, $$

which gives $\alpha'_i = \sum_j M^{-1}_{ij} \sigma_j$. Within the projection onto the reduced space, we then get

$$ 
\begin{aligned}
\langle \frac{\sigma_{eq}^2}{\sigma_0} \rangle := & \frac{3}{2V} \int_V \frac{1}{\sigma_0} \underline{\underline{\sigma}}' : \underline{\underline{\sigma}}' dV \\
= & \frac{3}{2V} \int_V \sigma_0 \frac{1}{\sigma_0} \underline{\underline{\sigma}}' : \underline{\underline{\sigma}}' dV \\
\approx & \frac{3}{2V} \int_V \sigma_0 \left( \sum_i \alpha'_i \underline{\underline{\epsilon}}_i^{in} \right) : \underline{\underline{\sigma}}' dV \\
= & \frac{3}{2V} \sum_{ij} M^{-1}_{ij} \sigma_j \underline{\underline{\sigma}}_i.
\end{aligned} 
$$
In matrix form, denoting by $\Sigma$ the column vector having $\sigma_j$ as coefficient in line $j$, we get

$$ \left\langle \frac{\sigma_{eq}^2}{\sigma_0} \right\rangle = \frac{3}{2V} \Sigma^t \mathbf{M}^{-1} \Sigma. $$

We can then simplify the right hand side of (17)

$$ \begin{aligned} & \int_V \frac{3m}{2\tau(x)} \left( \frac{\sigma_{eq}^2(x)}{\sigma_0^2(x)} \right)^{m/2-1} \underline{\underline{\sigma}}'(x) : \underline{\underline{\epsilon}}_j^{in}(x) dV \\ & \approx \frac{3m}{2} \left\langle \frac{\sigma_0^{1-m/2}}{\tau} \right\rangle \left( \left\langle \frac{\sigma_{eq}^2}{\sigma_0} \right\rangle \right)^{m/2-1} \int_V \underline{\underline{\sigma}}'(x) : \underline{\underline{\epsilon}}_j^{in}(x) dV \\ & \approx \frac{3m}{2} \left\langle \frac{\sigma_0^{1-m/2}}{\tau} \right\rangle \left( \frac{3}{2V} \Sigma^t \mathbf{M}^{-1} \Sigma \right)^{m/2-1} \underline{\underline{\sigma}}_j. \end{aligned} $$

From (17), the resulting reduced flow rule is then

$$ \dot{\alpha}_i = \frac{3m}{2} \left\langle \frac{\sigma_0^{1-m/2}}{\tau} \right\rangle \left( \frac{3}{2V} \Sigma^t \mathbf{M}^{-1} \Sigma \right)^{m/2-1} (\mathbf{M}^{-1} \Sigma)_i \quad (19) $$

which perfectly respects the variational structure of the original flow rule $\dot{\alpha}_i = \frac{\partial \Psi}{\partial \sigma_i}$ with a global dissipation potential given by

$$ \Psi(\Sigma) = V \left\langle \frac{\sigma_0^{1-m/2}}{\tau} \right\rangle \left( \frac{3}{2V} \Sigma^t \mathbf{M}^{-1} \Sigma \right)^{m/2} . $$

Observe also that in (19), because of the constitutive law (1), the jth line $\sigma_j := \int_V \underline{\underline{\sigma}} : \underline{\underline{\epsilon}}_j^{in} dV$ of $\Sigma$ is an explicit function of the expansion coefficients $\alpha_{jg}, \alpha_i, f_k$ and $\alpha_{sm}$ in (9), namely

$$ \sigma_j = \sum_{jg} \alpha_{jg} \int_V \underline{\underline{C}}_{\underline{\underline{\varepsilon}}} (\underline{\underline{\mathbf{U}}}_{jg}^{\text{ext}}) : \underline{\underline{\varepsilon}}_j^{\text{in}} dV + \sum_i \alpha_i \int_V \underline{\underline{C}}_{\underline{\underline{\varepsilon}}} (\underline{\underline{\mathbf{U}}}_i^{\text{in}}) : \underline{\underline{\varepsilon}}_j^{\text{in}} dV \\ + \sum_k f_k \int_V \underline{\underline{C}}_{\underline{\underline{\varepsilon}}} (\underline{\underline{\mathbf{U}}}_k^{\text{N}}) : \underline{\underline{\varepsilon}}_j^{\text{in}} dV + \sum_m \alpha_{sm} \int_V \underline{\underline{C}}_{\underline{\underline{\varepsilon}}} (\underline{\underline{\mathbf{U}}}_m^s + \underline{\underline{\mathbf{U}}}_m^s) : \underline{\underline{\varepsilon}}_j^{\text{in}} dV. \quad (20) $$

All integrals in this expression can be precomputed which results in considerable cost reduction.

### 3.4. Reduced friction law

The power developed by the friction force in the slipping mode $\underline{\underline{u}}_m^s$, as predicted by the friction law (7), is

$$ S_m^s := \sum_j \tau_t(J) \cdot \underline{\underline{u}}_m^s(J) = - \sum_j \min \left( \mu_t, \mu \frac{|\tau_N(J)|}{\|\dot{\underline{u}}_t(J)\|} \right) \dot{\underline{u}}_t(J) \cdot \underline{\underline{u}}_m^s(J). $$

We will temporarily omit the local minimization. But, even then, we still need to handle at each slipping contact point a division by the local value of the slipping velocity. To get rid of this division, we first consider the case where the slipping velocity is unimodal $\dot{\underline{u}}_t = \dot{\alpha}_l \underline{\underline{u}}_l^s$. The power developed by friction in such a slipping mode is then

$$ S_m^s (\dot{\alpha}_l \underline{\underline{u}}_l^s) = - \sum_j \mu |\tau_N(J)| \frac{1}{\|\underline{\underline{u}}_l^s(J)\|} \underline{\underline{u}}_l^s(J) \cdot \underline{\underline{u}}_m^s(J). $$

For multimodal slipping velocities $\dot{\underline{u}}_t = \sum_l \dot{\alpha}_l \underline{\underline{u}}_l^s$, we then proceed by “spherical interpolation”

$$ S_m^s \left( \sum_l \dot{\alpha}_l \underline{\underline{u}}_l^s \right) \approx - \sum_l \frac{\dot{\alpha}_l}{\| \dot{\alpha}_l \|} S_m^s (\dot{\alpha}_l \underline{\underline{u}}_l^s). $$
Taking now the minimum value (in norm) between the sticking and the slipping case at the modal level and developing the normal force $\tau_N$ along the given normal modes, we arrive to our proposed approximate friction law

$$S_m^s \left( \sum_l \dot{\alpha}_l \underline{u}_l^s \right) \approx F_m^s (\{\dot{\alpha}_l\}_l) := -\min \left( \mu_t \dot{\alpha}_m, \mu \sum_l \frac{\dot{\alpha}_l}{\|\dot{\alpha}\|} \sum_k f_k \left[ \sum_j \Lambda_k(J) \frac{\underline{u}_l^s(J) \cdot \underline{u}_m^s(J)}{\|\underline{u}_l^s(J)\|} \right] \right). \quad (21)$$

The terms between the square brackets can be precomputed which results in considerable savings. On the other hand, by writing the local equilibrium equations in weak form against all slipping modes, we must have

$$S_m^s = \langle \mathbf{KU}, \mathbf{U}_m^s \rangle.$$

After linear expansion of $\mathbf{U}$ by (9), the modal nonlinear friction law becomes

$$\left\langle \mathbf{K} \left( \sum_{jg} \alpha_{jg} \hat{\mathbf{U}}_{jg}^{\text{ext}} + \sum_{i} \alpha_{i} \hat{\mathbf{U}}_{i}^{\text{in}} + \sum_{k} f_{k} \hat{\mathbf{U}}_{k}^{\text{N}} + \sum_{m} \alpha_{sm} (\hat{\mathbf{U}}_{m}^{\text{s}} + \mathbf{U}_{m}^{\text{s}}) \right), \mathbf{U}_{m}^{\text{s}} \right\rangle + \mu_{t} \dot{\alpha}_{m} = R_{m}^{\text{s}} \quad (22)$$

with a nonlinear relaxation of slipping modes given by

$$R_m^s := F_m^s (\{\dot{\alpha}_l\}_l) + \mu_t \dot{\alpha}_m. \quad (23)$$

# 4. Numerical results

## 4.1. Global algorithm

After reduction, the problem to solve at all time $t$ is given by the reduced nonlinear laws ((19), (14)-(15), (22)-(23)) controlling the time evolution of the internal variables $\alpha_i$, $f_k$ and $\alpha_{sm}$ characterizing the creep strains, contact forces, and slipping displacements, respectively, and generating a displacement field and stress field given by (9).

This problem is solved by an implicit time integration algorithm, where at each time step nested prediction corrections algorithms are used to identify the nonlinear components in these laws. In more details, the proposed algorithm at each time step is as follows.

(1) External loop on creep variables

(2) Predict the values of the modal creep variables $\{\alpha_i(t)\}_i$

(a) Internal loop on the slipping relaxation

(b) Predict the value of the slipping relaxation $\{R_m^s(t)\}_m$

(c) Solve ((14)-(15), (22)) with respect to $\{f_k\}_k$ and $\{\alpha_{sm}\}_m$ with given values $\{\alpha_i(t)\}_i$ of creep modal coefficients and of slipping relaxation $\{R_m^s(t)\}_m$ using an active set strategy.

(d) Correct the value of the slipping relaxation $\{R_m^s(t)\}_m$ by (23), and repeat the internal loop until reaching convergence of this field.

(3) Integrate the evolution equation (19) in $\{\alpha_i(t)\}_i$ from time $t-dt$ to time $t$ by an implicit Runge Kutta algorithm with given values of $\{\alpha_{jg} \}_{jg}$, $\{f_k\}_k$ and $\{\alpha_{sm}\}_m$, correct the predicted value of $\{\alpha_i(t)\}_i$, and repeat the external loop until reaching convergence of this field.

This algorithm can be applied to situations involving many interacting fuel assemblies. In this case, the grid forces become a nonlinear function of the averaged grid displacements $\alpha_{jg}$ and step (2c) complemented with (10) become a global system of equations to be solved in all assemblies for all unknowns $\{\alpha_{jg}\}_{jg}$, $\{f_k\}_k$ and $\{\alpha_{sm}\}_m$.
## 4.2. Results for a fuel assembly creep

A total of 88 long-term creep simulations were performed offline with the fuel assembly detailed model under realistic neutron flux and for different force loadings applied constantly on grids [26]. A total of 176 snapshots were collected for creep strain and as many for the mean grids displacements. Two truncated POD then generate the 18 macroscopic modes and the 18 inelastic tensor modes, partly displayed in Figure 5. Now considering a central core position for the fuel assembly, the modal creep characteristics are determined from the corresponding neutron flux map. The grid displacement are evaluated with both models and compared in Figure 6 for 7 arbitrary load cases, which represent in plane loadings on grids. The inelastic deformation fields are also recombined for the load case number 6 and compared, in Figure 7, to those of the reference model. In this case, the ROM shows an excellent mapping although a slight deficit in axial creep intensity. Creep bending of the rods at the passage of each grid presents greater values and correlates very well. We obtain a speed-up of around 50 for this Fuel Assembly (FA) creep simulation. The reason for this modest speed-up is that the FA detailed reference model is already highly optimized and uses exclusively generalized finite elements, such as pipes and shells. Obviously, the same geometry represented by a classical mesh of linear or quadratic brick elements would display a much better time ratio in favor of the reduced model.

<figure><img src="image_10.png" /><figcaption>**Figure 5.** The 18 macroscopic FA displacement modes (3D) and the 5 first inelastic tensor modes.</figcaption></figure>

## 4.3. Results for a 3D contact problem with friction

A first test considers a contact problem between two elastic slender three-dimensional structures of 100 meters high, of 15 meters and 10 meters deep respectively, and of 5 meters wide [26]. They are embedded at their base (Figure 8) and are initially in contact on their interface, where friction is controlled by a Coulomb's law with coefficient $\mu_t = 0.5$. The displacement units in this section are also in meters (m). Two forces F1 and F2 apply respectively at mid-height and at the top of the left structure in order to impose a given in plane or out of plane displacement at these points. The elastic solution induced by these six elementary imposed displacements without contact conditions between the two structures define six boundary modes modes $\tilde{U}_{fg}^{ext}$. In addition, 10 positive modes of contact forces are extracted by NMF from 200 simulations previously performed offline. Similarly 6 interface slipping modes are extracted by POD and are
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Displacements comparison between the ROM and the detailed FA model for 7 load cases. ROM in red, reference in blue. Magnification x200.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Creep comparison for the load case number 6. Left: axial creep strain. Right: creep curvature in the rods (rad/m).</figcaption></figure>

represented on Figure 8 (b), with the 3D extension of the first slipping mode being represented on Figure 8 (c).

Figures 9 and 10 display the results of the full and of the reduced order model for an imposed displacement $\underline{u}_1 = (0.0, 0.0, 0.0)$, $\underline{u}_2 = (1.0, 0.0, 0.0)$ and $\underline{u}_1 = (0.661, 0.47, 0.825)$, $\underline{u}_2 = (-0.237, -0.975, 1.393)$, respectively. The first loading induces a large opening between the structures while the second one generates significant transversal slipping. The penetration is negligible and the displacements are superposed while the slipping and contact forces are well represented.
<figure><img src="image_13.png" /><figcaption>**Figure 8.** 3D-contact problem with friction, from left to right a) face and side view of the structure, b) the 6 slipping modes, c) 3D extension of the first slipping mode.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Loading 1: $\underline{u}_1$ = (0.0, 0.0, 0.0), $\underline{u}_2$ = (1.0, 0.0, 0.0). From left to right, superposition of recombined displacements, interface contact forces and interface slipping from the reduced model (red) and the reference model (blue) [26].</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Loading 2: $\underline{u}_1$ = (0.661, 0.47, 0.825), $\underline{u}_2$ = (-0.237, -0.975, 1.393). From left to right, superposition of recombined displacements, interface contact forces and interface slipping from the reduced model (red) and the reference model (blue) [26].</figcaption></figure>
4.4. *Results with a fuel assembly with contact and friction*

In order to extract the elementary loading modes defined in Section 2.2, a series of 60 calculations are performed offline with the detailed 3D model of a fresh fuel assembly. In these calculations, a constant friction coefficient is set at $\mu_t = 0.5$ and creep is not allowed. Different sets of forces are progressively imposed onto the grids with realistic values, so that the deformed shape of the fuel assembly stays in a range compatible with power plant experimental observations. From the results obtained, a truncated POD defines 11 modes $u_{jg}^{ext}$, 21 sliding modes $u_m^s$ and a NMF decomposition builds 10 positive modes $\Lambda_k$ for internal contact forces. The resulting reduced order model is then used to simulate a load case where grid 6 is subjected to imposed displacements following 2 cycles of increasing amplitude up to +/-20mm. The grid 6 force-displacement curve of the reduced contact-friction model is then compared with that of the detailed model in Figure 11. It can be seen that the reduced order model reproduces very well the stiffness and hysteresis due to the sliding of the rods in the grids. Concerning the computation time, the speed-up establishes to 200.

<figure><img src="image_16.png" /><figcaption>**Figure 11.** Bending test, force-displacement comparison of the reduced order model (red) with the 3D-reference model (blue) [26].</figcaption></figure>

**5. Conclusion**

The numerical tests indicate that considerable cost reduction (a factor of 50 to 1000 for a single assembly) can be achieved by the present strategy of domain based model reduction while preserving engineering accuracy.

The proposed strategy is quite general. It embeds an error control in the modal construction of the boundary displacements and of the creep modes. The open problems concern on one hand the error analysis of the NTFA averaged material law and on the other hand the combined control of the modal expansion of the contact forces and of the slipping modes in order to improve robustness and accuracy.

**Conflicts of interest**

The authors have no conflict of interest to declare.
## References

[1] Á. Horváth, B. Dressel, "On numerical simulation of fuel assembly bow in pressurized water reactors", Nucl. Eng. Des. 365 (2013), p. 814-825.

[2] G. C., J. Diaz, P. Le Tallec, J. M. Allain, "Multiscale mechanical model based on patient-specific geometry: application to early keratoconus development", J. Mech. Behav. Biomed. Mater. 129 (2022), article no. 105121.

[3] R. Glowinski, P. Le Tallec, Augmented Lagrangians and Operator-Splitting Methods in Nonlinear Mechanics, SIAM Studies in Applied Mathematics, vol. 9, Society for Industrial and Applied Mathematics, 1989.

[4] R. Glowinski, T.-W. Pan, T. I. Hesla, D. D. Joseph, J. F. Périaux, "A fictitious domain approach for the direct simulation of incompressible fluid flow past moving rigid bodies: Application to particulate flow", J. Comput. Phys. 169 (2001), no. 2, p. 363-426.

[5] R. Glowinski, O. Pironneau, "Numerical methods for the first biharmonic equation and for the two-dimensional Stokes problem", SIAM Rev. 21 (1979), no. 2, p. 167-212.

[6] A. Caboussat, R. Glowinski, D. C. Sorensen, "A least-squares method for the numerical solution of the Dirichlet problem for the elliptic monge-ampère equation in dimension two", ESAIM, Control Optim. Calc. Var. 19 (2013), no. 3, p. 780-810.

[7] R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer Series in Computational Physics, Springer, 1984.

[8] J. Diani, P. Le Tallec, "Variational upscaling for modeling state of strain-dependent behavior and stress-induced crystallization in rubber-like materials", Contin. Mech. Thermodyn. 33 (2021), p. 749-766.

[9] B. Leturcq, P. Le Tallec, S. Pascal, O. Fandeur, J. Pacull, "NTFA inspired model order reduction including contact and friction", in 14th. World Congress on Computational Mechanics (WCCM 14), 2021.

[10] B. Leturcq, P. Le Tallec, S. Pascal, O. Fandeur, N. Lamorte, "A new reduced order model to represent the creep induced fuel assembly bow in PWR cores", Nucl. Eng. Des. 394 (2022), article no. 111828.

[11] J.-C. Michel, P.-M. Suquet, "Nonuniform transformation field analysis", Int. J. Solids Struct. 40 (2003), no. 25, p. 6937-6955.

[12] A. Billerey, "Evolution of fuel rod support under irradiation — Impact on the mechanical behaviour of fuel assemblies", in Structural behaviour of fuel assemblies for water cooled reactors – Proceedings of a technical meeting held in Cadarache, France, 22–26 November 2004, no. IAEA-TECDOC-1454, International Atomic Energy Agency, 2005, p. 101-112.

[13] A. Wanninger, M. Seidl, R. Macián-Juan, "Mechanical analysis of the bow deformation of a row of fuel assemblies in a PWR core", Nucl. Eng. Technol. 50 (2018), no. 2, p. 297-305.

[14] P. Le Tallec, M. Vidrascu, "Efficient solution of mechanical and biomechanical problems by domain decomposition", Numer. Linear Algebra Appl. 6 (1999), no. 7, p. 599-616.

[15] G. J. Dvorak, "Transformation field analysis of inelastic composite materials", Proc. R. Soc. Lond., Ser. A 437 (1992), no. 1900, p. 311-327.

[16] J. Fish, K. Shek, "Finite deformation plasticity for composite structures: Computational models and adaptive strategies", Comput. Methods Appl. Mech. Eng. 172 (1999), no. 1-4, p. 145-174.

[17] J.-C. Michel, P.-M. Suquet, "A model-reduction approach to the micromechanical analysis of polycrystalline materials", Comput. Mech. 57 (2016), no. 3, p. 483-508.

[18] J.-C. Michel, P.-M. Suquet, "Computational analysis of nonlinear composite structures using the nonuniform transformation field analysis", Comput. Methods Appl. Mech. Eng. 193 (2004), no. 48-51, p. 5477-5502.

[19] Y. Liang, H. P. Lee, S. Lim, W. Lin, K. Lee, C. Wu, "Proper Orthogonal Decomposition And Its Applications—Part I: Theory", J. Sound Vib. 252 (2002), no. 3, p. 527-544.

[20] M. Balajewicz, D. Amsallem, C. Farhat, "Projection-based model reduction for contact problems", Int. J. Numer. Methods Eng. 106 (2016), no. 8, p. 644-663.

[21] IRSN, "IRSN documentation", 2021, https://www.irsn.fr/dechets/dechets-radioactifs/PublishingImages/assemblage-combustible.jpg.

[22] A. Franzén, "Evaluation of Fuel assembly Bow Penalty Peaking Factors for Ringhals 3 Based on a Cycle Specific Core Water Gap Distribution", PhD Thesis, Uppsala Universitet, Sweden, 2017, http://www.diva-portal.org/smash/get/diva2:1166550/FULLTEXT01.pdf.

[23] K.-n. Song, S.-h. Lee, "Effect of Weld Properties on the Crush Strength of the PWR Spacer Grid", Sci. Technol. Nucl. Install. 2012 (2012), no. 9-11, article no. 540285 (12 pages).

[24] D. D. Lee, S. Seung, "Learning the Parts of Objects by Non-Negative Matrix Factorization", Nature, London 401 (1999), no. 6755, p. 788-791.

[25] P. Paatero, "Least squares formulation of robust non-negative factor analysis", Chemometr. Intell. Lab. Syst. 37 (1997), no. 1, p. 23-35.

[26] B. Leturcq, "Réduction de modèles thermomécaniques non linéaires d'assemblages combustibles en vue d'un calcul de cœur REP", PhD Thesis, Institut Polytechnique de Paris, Paris, France, 2022.