Development of a numerical model of the POSEIDON irradiator for qualification in Co-60 radiosterilisation

CEA/Saclay research centre has several 60Co irradiation facilities dedicated to gamma irradiation for both CEA and industrials R&D needs in various fields such as electronuclear, defence, electronics, space as well as health applications.

In the more specific field of health, an irradiator is used to radiosterilise, i.e. to neutralise microbiological contaminants (such as bacteria, viruses, microbes, spores) using ionising radiation, for medical devices such as hip prostheses, orthopaedic screws or plates, on behalf of their suppliers. The great advantage of gamma radiation sterilization, compared with other sterilization alternatives (gas or cold immersion in liquid chemicals), is that medical devices do not have to be removed from their sealed pouches; they are processed directly through their packaging.

Radiosterilization of medical devices is a highly demanding process, in line with the requirements of ISO 13485 and ISO 11137. Firstly, the doses delivered must be neither too low, to ensure product sterility, nor too high, to avoid altering their integrity. Secondly, three qualification stages are required to guarantee validation of irradiation sterilization processes. The first two, known as installation and operational qualification, are respectively designed to demonstrate that the sterilization equipment has been installed in accordance with its specifications and is operating correctly, delivering appropriate doses within the limits of defined acceptance criteria. In particular, operational qualification consists in characterizing the irradiator in terms of dose distribution and reproducibility, by considering the volumes to be irradiated filled with a homogeneous material, with envelope densities representative of the products to be treated. Finally, the third qualification stage, known as performance qualification, must demonstrate, specifically for the medical products to be treated, that the equipment operates consistently, in accordance with predetermined criteria, and delivers doses within the specified dose range, thus producing a product that meets the specified requirement for sterility.

Depending on the supplier, irradiation packaging cartons are generally filled with a variety of different medical products, corresponding to a wide range of sizes and weights. The effects on spatial dose distribution of all possible product loading configurations should therefore be examined, including for different fill rates of the cartons on the irradiator's dynamic conveyor. Finally, it should be noted that the qualification processes must be repeated following any modification likely to affect dose or dose distribution, and therefore each time the sources are restocked. These processes are currently carried out exclusively by measurement, using a multitude of dosimeters appropriately distributed within and on the surface of the packages.

The Laboratoire des Rayonnements Appliqués (DRMP/SPC/LABRA), in charge of the POSEIDON gamma irradiator dedicated to the radiosterilization of medical equipment at CEA/Saclay would like to have a digital tool enabling these validation processes to be carried out by simulation. One of the major benefits expected from the use of this tool is to considerably reduce the downtime of the POSEIDON irradiator, imposed by the experimental validation phases.

The aim of the present thesis is to implement a numerical model of the POSEIDON irradiator, enabling the validation phases to be reproduced by simulation, as quickly as possible, while ensuring the accuracy of the results, to the desired precision. This work will be carried out at the DM2S/SERMA/CP2C laboratory (CEA/Saclay) with regular exchanges with the LABRA laboratory. CP2C is specialized, among other things, in radiation protection studies using numerical simulations.

Thus, the subject of the thesis, divided into three stages, will explore an alternative validation approach to that, carried out experimentally:

• The first stage will involve the development of a numerical model of the POSEIDON irradiator, integrating the dynamic nature of radiosterilization treatments. This numerical model will be based on a calculation methodology to be decided during the thesis (Monte-Carlo or deterministic method), with a compromise between the quality of the results obtained and the speed of calculation execution. For this stage, the radiation transport Monte Carlo code TRIPOLI-4® will be used as a reference, with comparisons made using other numerical tools such as MCNP®, PENELOPE, GEANT4, NARMER-1, etc.;

• The second stage will successively involve validation of the selected numerical model by comparison with experimental measurements, to be defined and carried out during the thesis, and its application to the calculation of operational qualification processes and performances for different families of supplier cartons. As regards validation of the calculations, the instrumentation used for gamma dose measurements will be numerically modelled and analysed, taking into account all the physical phenomena involved in absorbed dose (photon and electron doses). The aim is to consolidate calculation/measurement comparisons for experiments carried out during the thésis;

• The final step will be to analyze the contribution of the numerical model in relation to the experimental approach. This computational approach will nevertheless need to be optimized in terms of calculation time, in order to facilitate the sensitivity analyses to be carried out.

During the thesis, various directions of research will be investigated, such as improving the modelling of reflections during photon transport in a closed environment (PAGURE irradiator casemate; use of deep learning techniques for deterministic codes), implementing stochastic geometries to model the contents of the packaging to be irradiated, and improving algorithms to reduce computation times.

Hydrogen explosions within geometrically-tailored porous media : fluid-solid coupling and safety challenges


Hydrogen is a key asset for the energy transition, but it still poses major scientific and safety challenges. Colorless and odorless, hydrogen leaks easily, ignites at low concentrations and temperatures, and can lead to the propagation of rapid deflagrations as well as detonations, a dangerous type of supersonic combustion. Understanding the mechanisms involved in the transition from deflagration (slow flame) to detonation (supersonic flame accompanied by a shock wave) is therefore vital to the safety of hydrogen production facilities (electrolyzers) and the nuclear industry. In the accidental scenario of loss of cooling and core meltdown, oxidation of uranium rod cladding can lead to the release of hydrogen. It was the subsequent explosion that led to the loss of containment and release of radioactive material at Fukushima and Three Mile Island. Hydrogen risk management is therefore one of the major challenges for nuclear safety.

The main mechanism behind the deflagration -> detonation transition is the presence of obstacles along the flame path. These generate vorticity, which increases the surface area of the flame and accelerates the reactive wave. When obstacles are in sufficient number and proportion, a runaway effect and wave reflections can lead to a shock-chemical reaction coupling: detonation is born, propagating at several kilometers per second. Unfortunately, it's impossible to avoid the fact that industrial plants are cluttered with obstacles: pipes, buildings, machines, walkways, structures... and present this type of scenario.

Conversely, a very densely-packed, porous environment can, on the contrary, smother a too-rapid flame and allow the reverse transition from detonation to deflagration, which is less dangerous in nature. For example, detonation can be attenuated by passage through a porous matrix, or when a porous medium is placed along the walls during propagation in a tube. A crucial safety question then arises: under what circumstances does an obstacle accelerate or slow down a flame? Can porous media be designed to stop dangerous flames?


The aim of this thesis is to approach this question from three angles:

1/ on the one hand, via the preparation and execution of experimental tests on the CEA Saclay hydrogen explosion test bench (SSEXHY). These include:
- exploring different geometries for porous media, based on parameterizable topologies. These porous matrices will then be 3D printed via metal additive manufacturing;
- prepare instrumentation for the SSEXHY explosion test bench, featuring a visualization section using a Schlieren technique coupled with an ultra-fast camera capable of several million images per second;
- post-process the results of the shock and pressure sensors and the OH*-filtered photomultipliers.

2/ secondly, via numerical simulations of the DNS or LES type on research calculation codes. For example, we might be interested in :
- the influence of porous obstacle geometry (shape, porosity, hydraulic diameter, etc.) on flame propagation speed and deflagrationdetonation transitions;
- the influence of the 2/3D character of porous materials;
- the proposal of new criteria for selecting mesh refinement levels to capture the phenomena of interest.

3/ Finally, theoretical modelling of the problem from the point of view of volume-averaged equations will be carried out, with the aim of developing simplified, predictive models of the behavior of porous flame arresters.

Assessment of a new model for the investigation of an hypothetical Steam Generator Tube Rupture (SGTR) accident in Lead-cooled Fast Reactors (LFR).

The purpose of this PhD thesis is to implement and assess a new vapor explosion model in the
CEA code Europlexus, to investigate an hypothetical steam generator tube ruptures (SGTR) inside
the primary coolant pool of a lead-cooled fast reactors (LFR).

Advanced 2D and 3D neutron reflector modelling in support to high fidelity resolutions applied to different types of reactors

The modelling of the radial and axial interfaces of the core with the reflectors is source of approximations either for industrial Light Water Reactors (PWR and VVER types) but also for advanced concepts such as the Small Modular Reactors (SMR) and fast spectrum reactors (FSR) with an impact on the neutron flux distribution over the whole core.
In order to support the development of high fidelity resolutions to be applied to different reactor types, advanced 2D (for radial) and 3D (for axial) neutron reflectors models will be investigated,in this PHD. The analysis will start at level of the assembly (called lattice level) via an identification of relevant Figure of Merit (FoM) to be compared against reference calculations. The preliminary analisys will be completed by Monte Carlo reference calculations. The solutions will be then implemented in core calculations based on FEM methods such as the ones available on APOLLO3 codes for industrial-like homogeneous SP1 2 energy groups up to advanced SP3 pin-by-pin multi-groups solutions. This work will be realized in the framework of the neutronic code APOLLO3®. During his/her work the doctorand will be accompained by the technical team, and so he/she will be asked to intervene exclusively into the academic and theoretical part of the development work to support innovative industrial calculations and reference solutions.
Among the possible theoretical improvements of the classical models the following are identified:
• A new equivalence theory adapted to the Mixed Dual (Raviart Thomas) finite element method will be introduced in APOLLO3® Minos solver.
• In the classical diffusion approach, an equivalent (Benoist) diffusion coefficient is used. Unfortunately, the classical theory does not allow to treat (for example) the case when anisotropy is included in the emission density, or also when higher order expansion of the angular flux is done (SPN approximation). We propose the candidate to improve the classical model to include the possibility of defining anisotropic emission density in the core level and preserving the Benoist principle by using a generalization of this theory.

Staggered schemes for the Navier-Stokes equations with general meshes

The simulation of the Navier-Stokes equations requires accurate and robust numerical methods that
take into account diffusion operators, gradient and convection terms. Operational approaches have
shown their effectiveness on simplexes. However, in some models or codes
(TrioCF, Flica5), it may be useful to improve the accuracy of solutions locally using an
error estimator or to take into account general meshes. We are here interested in staggered schemes.
This means that the pressure is calculated at the centre of the mesh and the velocities on the edges
(or faces) of the mesh. This results in methods that are naturally accurate at low Mach numbers .
New schemes have recently been presented in this context and have shown their
robustness and accuracy. However, these discretisations can be very costly in terms of memory and
computation time compared with MAC schemes on regular meshes
We are interested in the "gradient" type methods. Some of them are based on a
variational formulation with pressure unknowns at the mesh centres and velocity vector unknowns on
the edges (or faces) of the cells. This approach has been shown to be effective, particularly in terms of
robustness. It should also be noted that an algorithm with the same degrees of freedom as the
MAC methods has been proposed and gives promising results.
The idea would therefore be to combine these two approaches, namely the "gradient" method with the same degrees of freedom as MAC methods. Initially, the focus will be on recovering MAC schemes on regular meshes. Fundamental
questions need to be examined in the case of general meshes: stability, consistency, conditioning of
the system to be inverted, numerical locking. An attempt may also be made to recover the gains in
accuracy using the methods presented in for discretising pressure gradients.
During the course of the thesis, time will be taken to settle the basic problems of this method (first and
second years), both on the theoretical aspects and on the computer implementation. It may be carried
out in the Castem, TrioCFD, Trust or POLYMAC development environments. The focus will be on
application cases that are representative of the community.

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.