Modeling of a shock on a filled nuclear power plant tank: Comparison of the fluid-structure model with and without ALE

[Hanna-che]

Modeling of a shock on a filled nuclear power plant tank: Comparison of the fluid-structure model with and without ALE

CSMI M2 Project - Hanna CHETOUANE
Company: Avnir Energy
Supervisors: Gabriel GEORGES, Guilherme VIANA, Christophe PRUD'HOMME, Vincent CHABANNES
Time Period: October 2025 - January 2026

Context

This project is part of the Sinudyn project, in collaboration with Avnir Energy, CEMOSIS, Sonohrc Technologies, and INSA Lyon. Its main goal is to develop a high-fidelity model of wave propagation in a tank following an impact, in the context of monitoring and maintenance of nuclear power plant tanks.

The physical model is based on fluid-structure coupling:

• the solid (tank structure) is modeled by the linear dynamic elasticity equation, describing the mechanical response following a hammer impact,
• the fluid is modeled by the linear dynamic wave equation, describing the propagation of pressure waves in the fluid.

Solid model

$$\rho_s \frac{\partial ^2 u_s}{\partial t^2} - \nabla \cdot \sigma(u_s) = f_s \quad \text{in } \Omega_s \times [0, T]$$
With:

• $u_s$: displacement of the solid (vector field, in m)
• $\sigma(u_s) = \lambda(\nabla \cdot u_s) I + 2 \mu \epsilon (u_s)$: stress tensor (Pa)
• $\epsilon (u_s) = \frac{1}{2}(\nabla u_s + \nabla u_s ^T)$: strain tensor (unitless)
• $\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}$ and $\mu = \frac{E}{2(1+\nu)} \quad$ where $E$ is Young's modulus (Pa) and $\nu$ is Poisson's ratio (unitless)
• $\rho_s$: density of the solid (in kg/ $m^3$)
• $f_s$: source term applied to the solid (in N/ $m^3$)

Fluid model

$$\frac{\partial^2 p}{\partial t^2} - c^2 \Delta p = f_f \quad \text{in } \Omega_f \times [0, T]$$
With:

• $p$: pressure field (scalar, in Pa)
• $c$: speed of sound in the fluid (in m/s)

Fluid-structure coupling

Coupling is achieved via boundary conditions at the fluid-structure interface $\Gamma_{fsi}$, where fluid pressure is applied to the solid and solid acceleration is applied to the fluid:
$$\sigma(u_s) \cdot n_s = -p n_s \text{,} \quad \frac{\partial p}{\partial n_f} = -\rho_f \frac{\partial^2 u_s}{\partial t^2} \cdot n_f \text{\quad on } \Gamma_{fsi}$$

Project objectives

Previous work (notably my M1 internship) used an Eulerian formulation with a fixed fluid mesh. The objective of this project is to extend the model by taking into account the displacement of the fluid mesh in the vicinity of the fluid-structure interface, integrating an ALE formulation. This will allow us to compare the Eulerian and ALE formulations in terms of accuracy, numerical stability, and computation time.
The study is conducted in a linear framework, assuming small deformations of the structure.

Adding ALE to the equations

The solid model remains unchanged, but it is necessary to propagate the displacement of the solid within the fluid domain. To do this, we introduce a diffusion equation to calculate the displacement of the fluid mesh:
$$- \nabla \cdot (k \nabla d_f) = 0 \quad \text{in } \Omega_f \times [0, T] \quad \text{with } \quad d_f = u_s \quad \text{on } \Gamma_{fsi} \times [0, T]$$

• $d_f$: mesh displacement (vector field, in m)
• $k$: diffusion parameter

The wave equation is modified by introducing terms related to the mesh velocity:
$$\frac{\partial^2 p}{\partial t^2} + \boldsymbol{2 \dot{d_f} \cdot \nabla (\frac{\partial p}{\partial t}) + (\nabla \cdot \dot{d_f}) \frac{\partial p}{\partial t} } - c^2 \Delta p = f_f \quad \text {in } \Omega_f \times [0, T]$$
with the same coupling conditions at the interface.

Implementation

The implementation was carried out in C++ using the Feel++ library. The equations are first solved monolithically, then coupled via Picard iterations, where:

• Solve the elasticity equation based on fluid pressure
• Solve the diffusion equation based on solid displacement
• Solve the wave equation with ALE based on solid acceleration
• Repeat until convergence: $|u^{n,k}-u^{n,k-1}| < \epsilon$ and $|p^{n,k}-p^{n,k-1}| < \epsilon$

Then at each time step:

• Execute the Picard loop
• Move the fluid mesh: $\Omega_f^n = \Omega_{\text{ref}} + d_f^n$, via the MoveMesh method provided by Feel++
• Update the fluid mesh $\Omega_f$
• Export the results and move to the next time step

Numerical simulations

The simulations were performed in 2D and 3D on simplified geometries of a reference tank. We considered stainless steel for the tank structure with water inside. The hammer impact was modeled by a wavelet and integrated via the fluid source function:
$f(x,t) = A \psi (t)$, with:

• $\psi(t) = -\sin(2 \pi f_c t) \exp(-5(f_c t -2)^2)$
• $A = 10^4$ N
• $f_c = 5000$ Hz

In the simulations considered, the mesh displacement is of the order of $10^{-12}$. It is therefore not directly visible in Paraview and has been amplified using the Warp by Vector filter.
In the following visualization, you can see the displacement of the fluid mesh with the evolution of wave propagation.

The following visualization shows the evolution of wave propagation in the tank with interaction with the structure, in 3D, with the fluid on the left and the solid on the right.

Comparison of the Eulerian and ALE models

Computation time

The ALE formulation entails a significant numerical cost due to the modified fluid solver and the additional resolution of the diffusion equation. The mesh displacement itself is not very costly.

Stability

The number of Picard iterations is similar in both cases, the residuals decrease rapidly, and the linear solvers (GAMG) converge in the same number of iterations. Thus, the integration of ALE does not degrade numerical stability.

Accuracy

For this case, we performed a time convergence study comparing the absolute and relative $L_2$ norms with models with and without ALE, as illustrated in the following graph.

The results show identical convergence for both cases. No measurable gain is observed with ALE, mainly due to the small displacement of the mesh, which shows no noticeable difference on the model.

Conclusion

In the context studied, i.e., for a linear model with small deformations, the introduction of ALE does not provide sufficient accuracy gains to justify the associated numerical overhead. The Eulerian formulation therefore remains sufficient for this type of configuration.
For a more complete validation, it would be necessary to perform a 3D simulation with a sufficiently fine mesh ($\approx$ 3mm) to compare the numerical results with experimental data obtained by Sonorhc.

Useful links

• 📄Project report: Hanna_CHETOUANE_report.pdf
• 🎤 Presentation: Hanna_CHETOUANE_slides.pdf