Monte Carlo methods for the adjoint transport equation: application to radiation shielding problems

The Monte Carlo method is the reference approach for simulating the transport of neutrons and photons, particularly in the field of radiation shielding, due to the very low number of approximations that it introduces. The usual Monte Carlo strategy is based on the sampling of a large number of particle histories, which start from a source, follow the physical laws of collision available in nuclear data libraries and explore the geometry of the system : the contributions of the particles to the response of interest (e.g. a count rate in a detector), averaged over all simulated stories, estimate the value predicted by the Boltzmann equation. If the detector region is "small", statistical convergence of the standard Monte Carlo approach becomes very difficult, because only an extremely limited number of stories will be able to contribute. It then becomes advantageous to use Monte Carlo methods for the solution of the adjoint transport equation: the histories of the particles are sampled from the detector backwards, and the collection region is the source of the starting problem (which is typically assumed to be “large” relative to the detector). This approach, simple in principle, offers the possibility of considerably reducing the statistical uncertainty. However, the adjoint Monte Carlo methods present scientific obstacles that are both practical and conceptual: how to sample the physical laws of collision “backwards”? How to control the numerical stability of adjoint simulations? In this thesis, we will explore different strategies in order to provide answers to these questions, in view of applying these methods to radiation shielding problems. The practical implications of this work could open up very encouraging perspectives for the new TRIPOLI-5® simulation code.

CFD development and modeling applied to thermal-hydraulics of hydrogen storage in salt caverns

A PhD thesis is available at LMSF lab of CEA in collaboration with Storengy, a world specialist in natural gas storage in salt caverns. Measurements carried out in the cavity showed that gas is in convective motion in the upper part of the cavity and is not necessarily in thermodynamic equilibrium with the brine at the bottom of the cavity, leading to gas stratification phenomena. The different flow regimes (convective or not) will strongly influence, on the one hand, mass exchanges between the gas and the brine and therefore the evolution of the gas composition (in moisture and other components) at the cavity exit and, on the other hand, thermal exchanges between the gas and the rock mass surrounding the cavity. In this context, CFD-based prediction tools are highly beneficial for understanding these phenomena and will contribute to a better interpretation of the physical measurements made in the cavity, to the design improvment of surface installations and to monitoring storage facilities, particularly for hydrogen storage. In this doctoral project, the aim is to develop a thermal-hydraulics model based on TrioCFD software for gas storage in realistically-shaped cavities and under cavity operating conditions (injection and withdrawal phases). To this end, the operation of storage salt cavities will be modeled, initially for a real geometry and in single-phase flow, then in two-phase flow, taking into account mass exchanges between the brine and the gas in the cavity.

High-Performance Computing (HPC) resolution of "point-saddle" problems arising from the mechanics of contact between deformable structures

In the field of structural mechanics, simulated systems often involve deformable structures that may come into contact. In numerical models, this generally translates into kinematic constraints on the unknown of the problem (i.e. the displacement field), which are dealt with by the introduction of so-called dual unknowns that ensure the non-interpenetration of contacting structures. This leads to the resolution of so-called "saddle-point" linear systems, for which the matrix is "indefinite" (it has positive and negative eigenvalues) and "sparse" (the vast majority of terms in this matrix are zero).

In the context of high-performance parallel computing, we're turning to "iterative" methods for solving linear systems, which, unlike "direct" methods, can perform well for highly refined numerical models when using a very large number of parallel computing processors. But for this to happen, they need to be carefully designed and/or adapted to the problem at hand.

While iterative methods for solving "positive definite" linear systems (which are obtained in the absence of kinematic constraints) are relatively well mastered, solving linear point-saddle systems remains a major difficulty [1]. A relatively abundant literature proposes iterative methods adapted to the treatment of the "Stokes problem", emblematic of incompressible fluid mechanics. But the case of point-saddle problems arising from contact constraints between deformable structures is still a relatively open problem.

The proposed thesis consists in proposing iterative methods adapted to the resolution of linear "saddle-point" systems arising from contact problems between deformable structures, in order to efficiently handle large-scale numerical models. The target linear systems have a size of several hundred million unknowns, distributed over several thousand processes, and cannot currently be solved efficiently, either by direct methods, or by "basic" preconditioned iterative methods. In particular, we will validate the approach proposed by Nataf and Tournier [2] and adapt it to cases where the constraints do not act on all the primal unknowns.

The work carried out can be applied to numerous industrial problems, particularly in the nuclear industry. One example is the case of fuel pellets, which expand under the effect of temperature and the generation of fission products, and come into contact with the metal cladding of the fuel rod, which can lead to cladding failure [3].

This thesis is in collaboration with the LIP6 laboratory (Sorbonne-université).

An internship can be arranged in preparation for thesis work, depending on the candidate's wishes.

[1] Benzi, M., Golub, G. H., & Liesen, J. (2005). Numerical solution of saddle point problems. Acta numerica, 14, 1-137. (
[2] Nataf, F., & Tournier, P. H. (2023). A GenEO Domain Decomposition method for Saddle Point problems. Comptes Rendus. Mécanique, 351(S1), 1-18. (
[3] Michel, B., Nonon, C., Sercombe, J., Michel, F., & Marelle, V. (2013). Simulation of pellet-cladding interaction with the pleiades fuel performance software environment. Nuclear Technology, 182(2), 124-137. (

Implicit/explicit transition for numerical simulation of Fluid-Structure Interaction problems treated by immersed boundary techniques

In many industrial sectors, rapid transient phenomena are involved in accident scenarios. An example in the nuclear industry is the Loss of Primary Coolant Accident, in which an expansion wave propagates through the primary circuit of a Pressurized Water Reactor, potentially vaporizing the primary fluid and causing structural damage. Nowadays, the simulation of these fast transient phenomena relies mainly on "explicit" time integration algorithms, as they enable robust and efficient treatment of these problems, which are generally highly non-linear. Unfortunately, because of the stability constraints imposed on time steps, these approaches struggle to calculate steady-state regimes. Faced with this difficulty, in many cases, the kinematic quantities and internal stresses of the steady state of the system under consideration at the time of occurrence of the simulated transient phenomenon are neglected.

Furthermore, the applications in question involve solid structures interacting with the fluid, undergoing large-scale deformation and possibly fragmenting. A immersed boundary technique known as MBM (Mediating Body Method [1]) recently developed at the CEA enables structures with complex geometries and/or undergoing large deformations to be processed efficiently and robustly. However, this coupling between fluid and solid structure has only been considered in the context of "fast" transient phenomena treated by "explicit" time integrators.

The final objective of the proposed thesis is to carry out a nominal regime calculation followed by a transient calculation in a context of fluid/immersed-structure interaction. The transient phase of the calculation is necessarily based on "explicit" time integration and involves the MBM fluid/structure interaction technique. In order to minimize numerical disturbances during the transition between nominal and transient regimes, the calculation of the nominal regime should be based on the same numerical model as the transient calculation, and therefore also rely on an adaptation of the MBM method.

Recent work defined an efficient and robust strategy for calculating steady states for compressible flows, based on "implicit" time integration. However, although generic, this approach has so far only been tested in the case of perfect gases, and in the absence of viscosity.

On the basis of this initial work, the main technical challenges of this thesis are 1) to validate and possibly adapt the methodology for more complex fluids (in particular water), 2) to introduce and adapt the MBM method for fluid-structure interaction in this steady-state calculation strategy, 3) to introduce fluid viscosity, in particular within the framework of the MBM method initially developed for non-viscous fluids. At the end of this work, implicit/explicit transition demonstration calculations with fluid-structure interaction will be implemented and analyzed.

An internship can be arranged in preparation for thesis work, depending on the candidate's wishes.

[1] Jamond, O., & Beccantini, A. (2019). An embedded boundary method for an inviscid compressible flow coupled to deformable thin structures: The mediating body method. International Journal for Numerical Methods in Engineering, 119(5), 305-333.

Seismic analysis of the soil-structure interface: physical and numerical modelling of global tilting and local detachment

The effect of soil-structure interaction, currently not taken into account in the seismic design of civil engineering structures and their foundations in professional practice, could influence the design of load-bearing structures. Soil-structure interaction effects are linked to inertial interaction (forces in the soil-structure system) and kinematic interaction (influence of the soil-foundation contact surface) (Semblat and Pecker, 2009). A more precise analysis of these two aspects requires three-dimensional (3D) numerical modeling of the soil-foundation-structure system and its temporal response, the definition of relevant constitutive laws for materials in the nonlinear plastic regime and the characterization of their mechanical properties. This makes it possible to consider directly the reduction in soil bearing capacity resulting from loss in soil strength, and the modification over time of the seismic action at the base of the structure. In 3D soil-structure interaction models, the connection between the soil and the embedded part of the footing structure is generally considered to be rigid, and the effects of friction and up-lift are neglected.
A series of tests on the CEA's Azalée shaking table in October 1999 (CAMUS IV, Combescure and Chaudat, 2000) demonstrated the existence of complex phenomena of structural detachment from the ground during seismic shaking. This involved a 1/3 scale model of the structure resting on a sandbox, anchored to the shaking table. Tests revealed up lifting at the foundation level, leading to energy dissipation, as well as significant residual settlement and rotation. Other studies have also highlighted the significant impact of rocking and consequent up-lift at the soil-foundation interface on the seismic response of the structure (Abboud, 2017; Chatzigogos, 2007; Gajan et al., 2021; Gazetas and Apostolou, 2004), as well as the loss of elasticity and nonlinear behavior of the soil, which increases permanent settlement (Pelekis et al., 2021). However, few studies in the literature evaluate the effect of the roughness of the soil-foundation interface and propose contact laws to model settlement and up lifting during rocking of the structure under seismic action.
In the context of interaction effects, understanding the parameters influencing the behavior of the soil-foundation interface and modeling the contact surface remains a challenge. A combined experimental and numerical approach will be developed in the proposed thesis.
The main aim of this thesis is to enable the transition from the modeling of local effects (friction, up lifting) to the simulation of the structure's global response (rocking, settlement, sliding). This is achieved by identifying the experimentally measurable physical parameters that manage the phenomenon locally and, at the same time, the global dynamic parameters altered by interaction effects (change in effective height).
On the one hand, an experimental campaign will be conducted on Vesuve, a single-axis vibrating table. The experimental model will consist of a rigid box containing the reference soil and a structure placed on the surface. The behavior of the system will be monitored using pressure sensors, LVDTs, flexiforce, accelerometers, etc. In addition, a numerical modeling method will be proposed and validated by comparison with experimental results. Finally, a numerical strategy will be proposed for different study cases. The output parameters obtained by the numerical simulations will be correlated with the measured parameters in order to optimize their calibration on the one hand, and to validate the numerical approach on the other.

Numerical simulation of turbulence models on distorted meshes

Turbulence plays an important role in many industrial applications (flow, heat transfer, chemical reactions). Since Direct Simulation (DNS) is often an excessive cost in computing time, Reynolds Models (RANS) are then used in CFD (computational fluid dynamics) codes. The best known, which was published in the 70s, is the k - epsilon model.
It results in two additional non-linear equations coupled to the Navier-Stokes equations, describing the transport, for one, of turbulent kinetic energy (k) and, for the other, of its dissipation rate (epsilon). ). A very important property to check is the positivity of the parameters k and epsilon which is necessary for the system of equations modeling the turbulence to remain stable. It is therefore crucial that the discretization of these models preserves the monotony. The equations being of convection-diffusion type, it is well known that with classical linear schemes (finite elements, finite volumes, etc ...), the numerical solutions are likely to oscillate on distorted meshes. The negative values of the parameters k and epsilon are then at the origin of the stop of the simulation.
We are interested in nonlinear methods allowing to obtain compact stencils. For diffusion operators, they rely on nonlinear combinations of fluxes on either side of each edge. These approaches have proved their efficiency, especially for the suppression of oscillations on very distorted meshes. We can also take the ideas proposed in the literature where it is for example described nonlinear corrections applying on classical linear schemes. The idea would be to apply this type of method on the diffusive operators appearing in the k-epsilon models. In this context it will also be interesting to transform classical schemes of literature approaching gradients into nonlinear two-point fluxes. Fundamental questions need to be considered in the case of general meshes about the consistency and coercivity of the schemes studied.
During this thesis, we will take the time to solve the basic problems of these methods (first and second year), both on the theoretical aspects and on the computer implementation. This can be done in Castem, TrioCFD or Trust development environments. We will then focus on regular analytical solutions and application cases representative of the community.

Greedy method for the model order reduction in neutronics : application of the reduced basis method

We are interested in a methodology that perform a computation in a very short amount of time while preserving the accuracy. A reduced basis approach could meet this requirement.
In the framework of the reduced basis methods [1,3], we devise an approximation space associated to a parameter-dependent partial differential equation. The construction of this approximation space includes a phase of exploration of the space of parameters where it is important to quantify the error between the solution obtained from the approximation space (in construction) and the solution obtained from a standard (fine) discretization.
This crucial step allows to certify the construction of the reduced basis.
Recently, some research work in the laboratory have provided a posteriori error estimator in the context of neutronics [4].
In this context [2], we are interested in generalized non-symmetric eigenvalue problems. Typically, we consider a linear Boltzmann operator of the form:
Find (u, v) such that Lu = Hu + v Fu,
where Lu is the advection operator, Hu is the scattering operator that modelize the collisions of the neutrons, Fu is the fission operator and the unknown u represents the neutron flux. The equation is also called the neutron transport equation. The fact the operator is not symetric comes form the scattering operator.
A first implementation of the reduced basis method based on the Proper Orthogonal Decomposition has been made for the neutron diffusion model in the APOLLO3® code [5]. Reduced basis methods have been studied for the neutron diffusion model [6-8] and the neutron transport model [9-14] with a varying degree of intrusivity.
The objective of this thesis is to contribute to the construction of greedy reduced basis methods for a model of neutronics, especially on assembling the reduced problem and the computation of an a posteriori estimator based on an affine decomposition of the operator. In a second step, many possibilities may be investigated :
- The extension of the reduced basis method to the simplified transport;
- The extension of the reduced basis method to the transport model;
- The application to the loading pattern of a research reactor.

[1] Y. Maday, O. Mula, A generalized empirical interpolation method: application of reduced basis techniques to data assimilation. Analysis and Numerics of Partial Differential Equations, XIII:221-231,2013.
[2] O. Mula, Some contributions towards the parallel simulation of time dependent neutron transport and the integration of observed data in real time, Chapter 1, 2014.
[3] G. Rozza, D. Huynh, and A. Patera, “Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations,” Archives of Computational Methods in Engineering, vol. 15, no. 3, pp. 1–47, 2008.
[4] Y. Conjungo Taumhas, G. Dusson, V. Ehrlacher, T. Lelièvre, F. Madiot. Reduced basis method for non-symmetric eigenvalue problems: application to the multigroup neutron diffusion equations. 2023. ?HAL cea-04156959?
[5] Y. Conjungo Taumhas, F. Madiot, T. Lelièvre, V. Ehrlacher, and G. Dusson. An Application of Reduced Basis Methods to Core Computation in APOLLO3®. M&C 2023
[6] Sartori, A. Cammi, L. Luzzi, M. E. Ricotti, and G. Rozza. Reduced order methods: applications to nuclear reactor core spatial dynamics.15566, in ICAPP 2015 Proceedings, 2015.
[7] S. Lorenzi, An adjoint proper orthogonal decomposition method for a neutronics reduced order model, Annals of Nuclear Energy, 114 (2018), pp. 245–
[8] P. German and J. C. Ragusa, Reduced-order modeling of parameterized multi-group diffusion k-eigenvalue problems, Annals of Nuclear Energy, 134
(2019), pp. 144–157
[9] I Halvic, JC Ragusa. Non-intrusive model order reduction for parametric radiation transport simulations. Journal of Computational Physics 492 (2023), 112385
[10] P Behne, J Vermaak, J Ragusa. Parametric Model-Order Reduction for Radiation Transport Simulations Based on an Affine Decomposition of the Operators. Nuclear Science and Engineering 197 (2), 233-261 (2023)
[11] P Behne, J Vermaak, JC Ragusa. Minimally-invasive parametric model-order reduction for sweep-based radiation transport. Journal of Computational Physics 469, 111525
[12] Z Peng, Y Chen, Y Cheng, F Li. A reduced basis method for radiative transfer equation. Arxiv preprint (2021).
[13] Sun, Y., Yang, J., Wang, Y., Li, Z., & Ma, Y. (2020). A POD reduced-order model for resolving the neutron transport problems of nuclear reactor. Annals of Nuclear Energy, 149, 107799.
[14] Wei, C., Di, Y., Junjie, Z., Chunyu, Z., Helin, G., Bangyang, X., ... & Lianjie, W. (2021). Study of non-intrusive model order reduction of neutron transport problems. Annals of Nuclear Energy, 162, 108495.

Toward robust simulations of uncertain industrial nonlinear dynamical systems

Vibrations are encountered in many components of nuclear reactors subject to fluid flows (steam generator for example). Resulting impacts and friction contacts can lead to wear of the materials which must be controlled for improved maintenance and to guarantee the lifespan of these industrial components.
Thus, it is necessary to know the solicitation experienced by the structure, define a reliable modeling of the system and implement efficient and accurate numerical methods to obtain the non-linear dynamic responses of the system.
In order to improve the robustness and performance of the numerical methods currently used at Framatome, CEA and EDF, this phD thesis have the following objectives:
- implement an efficient temporal integration algorithm suited for the case of systems with contact and friction, in particular to accurately reproduce “stick-and-slip” regimes and then the wear power of the surfaces in contact.
- demonstrate the ability to simulate both a single tube and a bundle of nearly 5000 interacting tubes,
- study the applicability of fluid-elastic coupling models identified on a single straight tube to a bundle of multi-supported 3D tubes,
- identify the relative influence of physical and numerical parameters by probabilistic approaches.

A Master 2 internship is planned before the phD as an introduction to this work.
Profile required: Final year engineer or Master 2 equivalent. Specialized in mechanics, attracted by numerical methods (both computer development and practical use).

Detection and diagnosis of anomalies in civil engineering structures using machine learning and numerical simulation

Monitoring reinforced concrete structures is of particular importance when it comes to identifying potential anomalies (cracking or excessive deformation, for example) in relation to nominal operation. These anomalies can have consequences both for the overall behavior (strength, etc.) and the functionality (tightness, etc.) of the structure. To meet this challenge (fault detection and prediction of consequences), a strong coupling between measurement data and simulations is essential. Current methodology relies mainly on initial instrumentation of the structure, based on expert opinion or feedback, but the data is not processed and analyzed using numerical calculation codes. The subject of the proposed thesis falls within the framework of a methodological breakthrough, through the combination of machine learning and numerical simulation tools for the detection and diagnosis of anomalies on civil engineering structures, in order to develop intelligent and adaptive instrumentation for monitoring the life of the structure. The methodology revolves around the following axes: processing of measurement data by machine learning, leading to the identification of potentially faulty zones, reconstruction of suitable boundary conditions around the previously detected anomaly by metamodeling, and identification of the fault and its consequences by numerical simulation. The thesis will be carried out jointly by two CEA laboratories: LM2S, specializing in structural mechanics, and LIAD, a unit specializing in Artificial Intelligence and Data Science.
The desired candidate (M2 level) should have an appetite for advanced numerical methods (including machine learning), as well as knowledge of mechanics and/or civil engineering. By the end of the thesis, the candidate will have developed knowledge and skills in numerical simulation, data assimilation and mechanics that can be effectively exploited in both industrial and academic environments.
The thesis may be combined with a preliminary M2 internship.

New condensation model in stratified flow at CFD and macroscopic scale by two-phase upscaling

In the context of safety of Pressurized Water Reactor (PWR), the Primary Coolant Loss Accident (LOCA) is of great importance. The LOCA is a hypothetical accident caused by a breach in the primary circuit. This leads to a pressure decrease in the primary circuit and a loss of water inventory in this circuit. Its resulting in heating of the fuel rods, which must remain limited so that damage to the fuel does not reduce cooling of the reactor core and prevents meltdown.

To remedy this situation, safety injection is activated to inject cold water, in the form of a jet, into the horizontal cold branch, which is totally or partially dewatered by the presence of pressurized steam. A stratified flow appears in the cold branch, with significant condensation phenomena in the vicinity of the jet and at the free surface in stratified flow zones. Numerous experimental and numerical works have been carried out on interfacial transfers at the free surface on rectangular and cylindrical cross-sections. CFD simulations of condensation at the free surface are carried out with the Neptune_CFD code, used by FRAMATOME, EDF and CEA. Currently, three models for heat transfer at the free surface are available in Neptune_CFD. These models have been established from a reduced number of simulations (DNS, LES and RANS) on rectangular configurations that remain far from the configuration of interest. Flows in a rectangular section tend to be parallel, whereas flows in a cylindrical section are three-dimensional.

The aim of this thesis is to improve the modeling of free surface condensation in a cylindrical cross-section configuration. Initially, a bibliographic study will be carried out on the free surface flow map, as well as on experimental works devoted to characterizing of interfacial area, mean interfacial velocity, turbulence terms in the vicinity of the free surface and heat transfer. In parallel, a new model will be developed in relation to the various improvement elements identified, and the associated validation carried out. Work is also planned to upscale two-phase CFD simulations to a macroscopic CATHARE approach. This up-scaling method will be based on Tanguy Herry's thesis work.