Wave propagation simulation in the upper core of sodium-cooled fast reactors using a spectral-element method for heterogeneous media

ASTRID project, French sodium cooled nuclear reactor of 4th generation, is under development at the moment by Alternative Energies and Atomic Energy Commission (CEA). In this project, development of monitoring techniques for a nuclear reactor during operation are identified as a measure issue for enlarging the plant safety. Use of ultrasonic measurement techniques (e.g. thermometry, visualization of internal objects) are regarded as powerful inspection tools of sodium cooled fast reactors (SFR) including ASTRID due to opacity of liquid sodium. In side of a sodium cooling circuit, heterogeneity of medium occurs because of complex flow state especially in its operation and then the effects of this heterogeneity on an acoustic propagation is not negligible. Thus, it is necessary to carry out verification experiments for developments of component technologies, while such kind of experiments using liquid sodium may be relatively large-scale experiments. This is why numerical simulation methods are essential for preceding real experiments or filling up the limited number of experimental results. Though various numerical methods have been applied for a wave propagation in liquid sodium, we still do not have a method for verifying on threedimensional heterogeneity. Moreover, in side of a reactor core being a complex acousto-elastic coupled region, it has also been difficult to simulate such problems with conventional methods. The objective of this study is to solve these 2 points by applying three-dimensional spectral element method. In this paper, our initial results on three-dimensional simulation study on heterogeneous medium (the first point) are shown. For heterogeneity of liquid sodium to be considered, four-dimensional temperature field (three spatial and one temporal dimension) calculated by computational fluid dynamics (CFD) with LargeEddy Simulation was applied instead of using conventional method (i.e. Gaussian Random field). This three-dimensional numerical experiment yields that we could verify the effects of heterogeneity of propagation medium on waves in Liquid sodium.

Abstract-ASTRID project, French sodium cooled nuclear reactor of 4th generation, is under development at the moment by Alternative Energies and Atomic Energy Commission (CEA).In this project, development of monitoring techniques for a nuclear reactor during operation are identified as a measure issue for enlarging the plant safety.Use of ultrasonic measurement techniques (e.g.thermometry, visualization of internal objects) are regarded as powerful inspection tools of sodium cooled fast reactors (SFR) including ASTRID due to opacity of liquid sodium.In side of a sodium cooling circuit, heterogeneity of medium occurs because of complex flow state especially in its operation and then the effects of this heterogeneity on an acoustic propagation is not negligible.Thus, it is necessary to carry out verification experiments for developments of component technologies, while such kind of experiments using liquid sodium may be relatively large-scale experiments.This is why numerical simulation methods are essential for preceding real experiments or filling up the limited number of experimental results.Though various numerical methods have been applied for a wave propagation in liquid sodium, we still do not have a method for verifying on threedimensional heterogeneity.Moreover, in side of a reactor core being a complex acousto-elastic coupled region, it has also been difficult to simulate such problems with conventional methods.The objective of this study is to solve these 2 points by applying three-dimensional spectral element method.In this paper, our initial results on three-dimensional simulation study on heterogeneous medium (the first point) are shown.For heterogeneity of liquid sodium to be considered, four-dimensional temperature field (three spatial and one temporal dimension) calculated by computational fluid dynamics (CFD) with Large-Eddy Simulation was applied instead of using conventional method (i.e.Gaussian Random field).This three-dimensional numerical experiment yields that we could verify the effects of heterogeneity of propagation medium on waves in Liquid sodium.Index Terms-Acoustics, Acoustic simulation, Finite element analysis, Spectral-element method, I. INTRODUCTION HE 4th generation nuclear reactors, which use liquid metals for medium of their coolant circuits, is under development as an international project [1,2,3].Inspection and repairing systems for inside of a reactor core are regarded as important issues for enlarging its safety.Based on this viewpoint, improvement of measurement tools' accuracy, sensitivity, frequency of measurements and reliability are required [8,6].Ultrasonic measurement systems have been studied and developed among LIET CEA Cadarache and LMA CNRS as possible methods for liquid sodium with no transparency [5].
It is known that strong flow fields and heterogeneities of temperature and density of gas bubbles caused by flows occurs above core sub-assemblies during operation of reactor from studies [9,25,26].Also the effects of these heterogeneities of gas density and temperature gradient on a wave propagation have been researched [20,7].By these studies, it was revealed that the effects of those heterogeneities are not negligible as for estimating wave propagations inside of reactor core.Because of the use of liquid metals that is a characteristic 4th generation of nuclear reactors, possible scales and the number of times are limited.Thus numerical methods for predicting this phenomenon is demanded.
There exist several studies on the numerical methods for a wave propagation in heterogeneous medium.D. Fiorina simulated the fluctuation of sound celerity caused by heterogeneity of a medium by applying Gaussian random field (GRF).He concluded that it is possible to estimate the amount of possible effects of medium fluctuation on a wave propagation only for limited propagation distance and strength of fluctuation, which was confirmed for two dimensional acoustic problems with a point and plane source and for isotropic fluctuation condition (i.e.no strong flow) [10,18].This method which applying GRF is recognized as a method which required relatively fewer computational cost than methods which use the CFD.B. Lü improved this method and applied for three dimensional problems by applying the fluctuation for the time-of-flight and signal amplitudes calculated by a ray (tracing) method beforehand [19].B. Iooss expanded this application of GRF for anisotropic fluctuation.He proposed to used anisotropic GRF for the fluctuation in twodimensional flow fields.At the same time, he concluded that the fluctuation pattern generated by Gaussian method is far from the real situation because this is a stochastic and idealized method [13].
Additionally, only the ray (tracing) method [23] or the pencil method [11] have been used as numerical methods for a wave propagation.There are limitations on application of ray tracing techniques for heterogeneous medium.There are limitations on application of ray tracing techniques for heterogeneous medium.
The objective of this study is to simulate the effect of heterogeneity on three-dimensional wave propagation by applying the spectral element method, which is one of higher order finite element methods.As the method to simulate the fluctuation of medium, we applied an CFD simulation results, which was calculated by CEA STMF, then we studied the effects of realistic fluctuated fields on the three-dimensional acoustic field from a circular piston plate.In the section 2, the model target of our numerical study is explained.In the section 3, the numerical tools which we used for the wave simulation, the details of CFD and interpolation of medium properties are described.Then we show our first results about the application of three-dimensional wave simulation with CFD.

II. NUMERICAL METHODS FOR A WAVE PROPAGATION
In our study, we simulate a wave propagation with three dimensional spectral-element method (SEM) for the discretization scheme in spatial dimension.We used a numerical code called SPECFEM.SPECFEM is a generic name of numerical codes, belonging CIG (Computational Infrastructure for Geodynamics), for a wave propagation simulation which uses the SEM.Initially this was developed by Dimitri Komatitsch and Jean-Pierre Vilotte at Institut de Physique du Globe (IPGP) in Paris, France from 1995 to 1997 and then by Dimitri Komatitsch and Jeroen Tromp at Harvard University and Caltech, USA for a seismic wave simulation.[16,17,30] On SPECFEM, continuous Galarkin typed SEM with higher (4 th )-order Gauss-Lobatto-Legendre is implemented.The higher order basis function brings better accuracy for surfacewave, acoustic-elastic coupling problems than the other types of finite-element method using lower-order polynomials.The time domain is discretized by Newmark-beta method.

A. PLAJEST experiment / CFD simulation
We selected an experimental and numerical study called PLAJEST as the modeling target of our study since PLAJEST is also targeting the same object of the same condition i.e. the upper-core region of SFR in operation.This study program is collaboratively worked between Japanese Atomic Energy  Agency (JAEA), the U.S. Department of Energy (DOE) and CEA.The configuration of this experiment is shown in Figure 1.The main purpose of PLAJEST is to observe heat conduction at a liquid-solid boundary in a sodium cooling circuit because the thermal fluctuations lead to high frequency thermal fatigue and crack in adjoining structures of SFR.
After the experiment by JAEA [14] [15], thermofluid analysis has been carried out by CEA/STMF [4].STMF used a CFD tool called TrioCFD (known as Trio_U by 2015) for their calculation and LES was employed as a turbulent model.Tetrahedral element with 4 nodes was selected for TrioCFD computation model.Total number of element reached 5,582,706 and characteristic mesh length was set to 1.40 mm.Initial 200 seconds of their calculation is excluded as duration for stabilization then 200 ~ 210 seconds with 0.1 milli-seconds of time step is available as candidates for our wave simulation.3 jets of sodium exist in this setup.Sodium with lower temperature (304.5 ℃) is emitted from central jet and with higher temperature (347.5 ℃) from two outer jets.The average flow velocity is 0.51 m/s for every jet.The simulated temperature field at time = 200.787seconds is also shown in figure 1 (b).Their simulation results achieved good accordance with the experiment results by JAEA in terms of Spectral power density and standard deviation of temperature values.[4]

B. Wave simulation in realistic fluctuated medium
For the modeling of sodium state, we took few hypothesis:  The density of gas bubble is constant everywhere in whole simulated region.


The effect of flow on a wave propagation is vanishingly small, so called frozen fluid.Thus, the heterogeneity of the temperature field is taken into account for our study.
As can be seen in figure 1 (b), the temperature field is spatially (and also temporally) stable in the region of just above the jets' outlet.We call this field as a static temperature field.The pattern of temperature field becomes more complex and temporal fluctuation occurs as the region goes farther from the outlets.We call this field as a fluctuation field.Our former study [21] applied the hypothesis that the temperature field in the region near the outlet may be described by superimposing a static field and fluctuation field.The static field was provided as very simple profiles including parabolic shapes.The fluctuation fields were generated by GRF.However, as described above, GRF method is only validated for the field without strong flows and we are still not sure about to what extent GRF is applicable for the region including flows.It is expected that GRF is almost perfectly applicable for the simulation region being sufficiently far from the outlets and the difference of acoustic simulation results between with GRF and with CFD will be increased when the region becomes closer the outlet.
Hence one of the objectives of this study is to examine the applicability of the GRF approximation for the medium fluctuation by comparing the results from GRF-applied wave simulation and CFD-applied simulation and changing the distance from the sodium outlets.
In this paper, as a previous step, we apply the temperature field calculated by CFD to our three-dimensional wave simulation, then observe the modified acoustic fields and relation between the distance from the outlets.
We use the Sobolev's empirical equations [24] for the relationships between temperature and wave celerity in sodium: = 2723.0− 0.531   (1) where   is the celerity of ultrasonic waves in meters per second and   is sodium temperature in Kelvin degrees.
Density is also temperature dependent:

C. Configuration of a three-dimensional wave simulation
In order to decrease the total computational costs, the target domain for our acoustic simulation is partially extracted from the entire region calculated in the CFD side and other domain which is acoustically not interested is excluded from our simulations.We extracted the rectangular parallelepiped region drawn as the purple box in figure 1.The boundary lengthes for y. z axis (perpendicular to the direction of propagation) are defined to become about three times of the diameter of the source plane.The dimension for x axis is taken enough for acoustic waves' passing through three sodium jets completely.Herewith, we observe the modification of acoustic fields before and after passing through each jet.
We carried out multiple simulations by changing the z coordinate of the simulation region from z = 0.040 meter (just above the jets) to z = 0.340 meter at 0.010 meter interval (z origin is set as the same level of the edges of outlets, which z value means the central coordinate of the simulated box).X and y coordinates are kept in all simulations.
As the source of all simulations, we used circular plane source, which is composed of monopole point sources on a circular plane with a diameter of 0.0254 meter.The intervals of each point sources are same as the element size of SPECFEM mesh (i.e.approximately 8.0 × 10 −4 meter).Each source point emits a 1MHz Ricker wavelet.The maximum amplitude of each emission is multiplied by hamming window function having a maximum point at the center of the circular plane.Figure 2 shows this source plane.The color shows the maximum amplitudes of the emissions at each point normalized by the maximum amplitude at the center of the circle.
The Finite-Element (FE) Mesh for SPECFEM3D is generated by FE meshing software CUBIT (developed by Sandia National Laboratories, United States).Both 1st and 2ndorder hexahedral FE may be imported to SPECFEM.We selected 1st-order elements for our simulation because the targeted model does not have any curved geometry.It is necessary to convert the data format while the mesh data is exported from CUBIT to SPECFEM3D.In order to speed-up this conversion process, we used IOSS (IO Systems) library included the FE analysis supporting software SEACAS (also developed by Sandia National Laboratories, United States).After being imported into the SPECFEM, the GLL points are inserted into each element.Determination of an element size is done following 2 conditions, which are CFL condition (3)  where  is the interval of each time increment.  is the minimum interval of GLL points.
We selected the averaged Courant number  = 0.4 and the wave celerity   = 2416.268m/s (in sodium with lowest temperature value 274.5℃ during the CFD simulation) for the calculation of a mesh size and time step length.In the equation 3, dx is not the mesh size itself, but it is the interval of GLL points.We decided the average mesh size  = 8.05 × 10 −4 meter, time step interval  = 2.3 × 10 −8 second.For the number of time step updates, we set 5000 steps for having enough temporal duration for waves' passing through entire simulation domain.The mesh model used our simulations has 3,250,000 elements and 3,722,406 nodes and total number of computation nodes including GLL points became 215,320,764.

D. Interpolation of temperature field from CFD to wave simulation
The element type of TrioCFD is a tetrahedral FE with the average element size  = 1.400 × 10 −3 meter.The temperature value is defined on the central point of each element, while the element type of SPECFEM3D is a hexahedral with average element size  = 8.050 × 10 −4 .The temperature values need to be defined on the corner nodes and every internal GLL point.Thus the CFD temperature data needs to be interpolated for importing to SPECFEM3D.For this process, we used a software library so called MEDCoupling, which is one of the libraries included in a scientific numerical method supporting tool SALOME.SALOME is the open source software developed by CEA.By using MEDCoupling, interpolation of physical fields between several types of FE or nodes is possible.In our application, the temperature field of TrioCFD is at first interpolated onto the corner nodes of SPECFEM3D elements using conservative interpolation.After being imported into SPECFEM3D, all temperature values on corner nodes are interpolated on the GLL points using tri-linear interpolation.
Figure 3 shows the interpolated temperature field on SPECFEM3D mesh.For visualization purpose, the mesh size is 10 times larger than the mesh actually used.

E. Computational resources
In this study, all of simulations were performed on 2 super computers, CURIE (CEA TGCC) and OCCIGEN (CINES).By the helps of those advanced computer systems and gigantic capacities dedicated for scientific calculations, the computation domain was divided into 256 parts then the parallelized calculations were carried out.The average duration for 1 acoustic simulation takes about 15 ~ 20 minutes excluding mesh generation and interpolation of temperature fields.

IV. RESULTS
Figure 4 shows the positions of the plane source (blue) and receiver surfaces (orange and red).For the visualization purpose, acoustic fields calculated here are indicated as crosssection surfaces.The acoustic fields and the effects from the temperature heterogeneity were observed for each simulation on x-y, y-z receiver planes and y-z receiver planes positioned at x = -0.105,-0.080, -0.035, 0.000, 0.035, 0.080, 0.105 meter.The y-z plane at x = 0.000 meter is the middle point of the simulation domain.The others are placed at front and back of each sodium jet.
Figure 5 shows part of acoustic fields obtained from all of the simulations.Here an acoustic field refers the maximum pressure values at each spatial point.In figure 5, acoustic fields on x-y, y-z and y-z plane (only x = 0.035, 0.105 meter) are indicated for the simulations whose middle z coordinate value of simulation domains are z = 0.04, 0.12, 0.24, 0.34 meter.Additionally, the top row is the results of homogeneous case.In the homogeneous case, the temperature of medium was set as 333.167℃ which is the ambient temperature value of all heterogeneity cases.From this figure, it is found that the acoustic field of z=0.040, which is close to the outlets, and homogeneous case are quite similar.It is because that the temperature boundary at just above the outlets is almost parallel to the wave front.Also at z = 0.340 meter, the mixing of the hot and cold jets matured enough and the magnitude of temperature difference becomes small; therefore, the wave less effected by the heterogeneity of this altitude.However, at z = 0.120 and 0.240 meter, it can be seen that the acoustic fields are bent by strong heterogeneity and the transmission loss towards propagation distance becomes greater.Fig. 3. the interpolated temperature field on SPECFEM3D mesh.For visualization purpose, the mesh size is 10 times larger than the mesh actually used.In figure 6, values of acoustic fields on z axis of each y-z receiver planes are shown.Red lines represent the results of homogeneous cases and blue lines represent all of the results from heterogeneous cases at each z position of simulation domains.Amplitude values are normalized by the maximum amplitude (i.e.pressure) among all simulations.Also, the values on horizontal axis are normalized by the diameter of the source plane.The larger x values become, i.e. the farther the yz receiver plane goes, the more dispersive among the amplitude curves become.How the maximum amplitude of pressure recorded on each y-z receiver plane changes depending the positional changes for z direction of simulation domains is indicated in figure 7.At planes x=-0.105,-0.080 meter, positioning near to the plane source, maximum amplitudes take approximately constant values regardless of distances from outlets.In contrast, at the farthest and second farthest y-z planes at x = 0.105, 0.080 meter, it was found that the amplitude value of simulation domain z = 0.120 meter is approximately 10 % smaller than the other values.At the same time the second peak may be recognized at z = 0.240 meter.
It is considered that the effect of temperature heterogeneity becomes stronger under 2 conditions, that are:  Mixing of hot and cold flow are grown enough so that the shape of temperature boundary is curved and not parallel to the emitted wave front. Temperature difference is kept large.At the altitude z = 0.120 meter, it is about the distance from outlets where the temperature boundary starts to be changed while the temperature difference is relatively still high compared to the region above.
On the other hand, it seems that, at z = 0.240 meter, the temperature difference there is already small and even smaller that z = 0.230 meter.This might be caused by the temperature boundary shape became a key factor of the total effects rather than the temperature difference.

V. CONCLUSION
This study observed the effects of three-dimensional heterogeneity of temperature field on a wave propagation under a high temperature liquid sodium environment by numerical simulation using spectral element method.
The heterogeneous temperature field were imported from CFD Large-eddy simulation which was carried out based on an experiment and validated its result by comparing with the experimental result.
As a result, it was found that there is a certain relationship between a strength of effects of temperature heterogeneity on acoustic fields and a distance from the outlets of sodium jets, i.e. the state of mixing flows.More concretely, the amount of effects on maximum acoustic pressure can vary up to about 10 % at a certain altitude.
Accordingly, we may assume that the magnitude of temperature difference in a medium and the shape of temperature boundary are the key factor which decide the effect of medium heterogeneity on a wave propagation.
One single time step was examined for this study.We will continue this observation for other time shots in order to obtain an index concerning a relation between the strength of effects from medium heterogeneity and the distance from the outlets.
Subsequently, by comparing with a wave propagation with a heterogeneity generated by GRF, we will verify the availability of isotropic assumption on heterogeneity which was taken in GRF for the medium where a sodium flow exists.Fig. 7.The maximum amplitude of pressure recorded on each y-z receiver plane changes depending the positional changes for z direction of simulation domains.The minimum peak may be found at z = 0.12 meter on the lines x = 0.08, 0.105 meter.The amplitude value of simulation domain z = 0.120 meter is approximately 10 % smaller than the other values.The second peak is also found at z = 0.24 meter.

Fig. 1 .
Fig. 1.(a) Geometry of PLAJEST CFD simulation.3 sodium jets outflow from the gaps with 20 mm width.The region for wave propagation simulation is extracted (purple box).Multiple simulations are carried out with changing the z coordinate of extracted region from 40 mm to 340 mm above the outlet at the intervals of 10 mm.(b) A snapshot of CFD result at time 200.0 seconds, x-y cross-sectional plane at y = 0.09 meter.The temperature pattern just above the outlets are simple and stable while it becomes more complex and unstable when z coordinate goes the farther.

Fig. 2 .
Fig. 2. The circular source plane applied for the simulations.The color shows the maximum amplitudes of the emissions at each point normalized by the maximum amplitude at the center of the circle.The plane source is composed of monopole point sources on a circular plane with a diameter of 0.0254 meter.The intervals of each point sources are same as the element size of SPECFEM mesh (i.e.approximately 8.0 × 10 −4 meter).Each source point emits a 1MHz Ricker wavelet.The maximum amplitude of each emission is multiplied by hamming window function having a maximum point at the center of the circular plane.

Fig. 4 .
Fig.4.The positions of the plane source (blue) and receiver surfaces (orange and red).For the visualization purpose, acoustic fields calculated here are indicated as cross-section surfaces.

Fig. 5 .Fig. 6 .
Fig.5.Acoustic fields obtained from all of the simulations.Here an acoustic field refers the maximum pressure values at each spatial point.Acoustic fields on x-y, y-z and y-z plane (only x = 0.035, 0.105 meter) are indicated for the simulations whose middle z coordinate value of simulation domains are z = 0.04, 0.12, 0.24, 0.34 meter.Additionally, the top row is the results of homogeneous case.In the homogeneous case, the temperature of medium was set as 333.167℃ which is the ambient temperature value of all heterogeneity cases.
and the number of elements par one wave length,