David L Clark

John C Gordon

P Jeffrey Hay
email: pjhay@lanl.gov

Richard L Martin

Rinaldo Poli

A DFT Study of Tris(bis(trimethylsilyl)methyl)lanthanum and -samarium

DFT calculations have been carried out on compounds Ln[CH(SiR2R')(SiR3)]3 for Ln = La and Sm and (i) R = R' = Me; (ii) R = H, R' = Me; (iii) R = R' = H. The results are compared with the X-ray structures that are available from the literature for both metals and R = R' = Me. The calculations correctly reproduce the experimental structural features in these complexes exhibiting the peculiar pyramidal coordination geometry. The results show significant increases in the Si-C bond lengths associated with -Si-C agostic interactions, whereas little structural changes are found for -C-H agostic interactions. The latter are in fact repulsive. The simplified model system with R = H and R' = Me that retains one agostic methyl interaction in each alkyl ligand also correctly reproduces the essential geometrical features. The simplest model with only SiH3 groups, while also adopting a pyramidal coordination geometry, no longer accurately describes the real molecule since the -Si-C agostic interactions are replaced by -Si-H interactions. A Mulliken analysis of the electronic structure shows a relatively covalent Ln-C interaction with significant 5d orbital participation in the bonding. A number of calculations with different basis sets on the model Sm[CH(SiH2Me)(SiH3)]3 system shows that the addition of polarization functions (d functions on Si and C, p functions on H, or f functions on Sm) has little or no beneficial effect on the quality of the results.

Introduction

In the late 1980s Power et al. reported two remarkable structures of homoleptic trialkyl complexes of lanthanum and samarium, Ln[CH(SiMe3)2]3 (Ln = La, Sm). 1 These compounds have a formal electron count of only 6 valence electrons if the (4f) n electrons are considered as core electrons. They adopt a highly distorted trigonal pyramidal structure, with C-Ln-C angles close to the tetrahedral value, and unusually short Ln•••Me contacts. Although the potential for the participation of agostic interactions (of the -CH type) was recognized, the molecular geometry was in fact rationalized on purely steric grounds, "the deviation from planarity occurring in order to maximize ligand-metal attractions and minimize ligand-ligand repulsions". 1 The isolation an structural characterization of the tris(amido) samarium(III) derivative Sm[N(SiMe3)2]3, 2 has led us to examine and re-evaluate the structural features of the isolobal and isostructural tris(bis(trimethylsilyl)methyl) complexes. A further impetus to our work was also given by the recent crystallographic and DFT study 3 on the related complex Cp*La{CH(SiMe3)2}2. Although this molecule is electronically more saturated relative to the title compounds, the presence of agostic interactions with the alkyl ligand -Si-C bonds was unequivocally established. Close contacts between lanthanide centers and the CH(SiMe3)2 ligand have been observed in numerous instances. 4,5 However, theoretical studies of these systems have remained rather limited, possibly because of the difficulty in handling the f shell for lanthanide systems. Recently, however, Maron and Eisenstein have shown that the use of "large core" relativistic effective core potentials (RECPs) for the lanthanide atom, i.e., one which folds the 4f n shell into the core, greatly simplify the computational effort and yet provide sufficiently accurate results for lanthanide model systems. [6][7][8] The above mentioned DFT study of the model CpLa{CH(SiMe3)2}2 system employed the ADF program in which the core electrons are present in the form of a relativistic frozen core potential. The calculations provided a very good match between the optimized and experimental structures, in particular the metric parameters related to the metal-ligand agostic interactions.

Computational Details

All calculations were carried out using the B3LYP functional and employing relativistic effective core potentials (RECP) on the La and Sm atoms. Unless specifically noted in the text, all calculations used the "large core" RECP in which the 5s 2 5p 6 6s 2 5d 1 electrons were explicitly treated as "valence" electrons with the remaining electrons replaced by the RECP. 9,10 In the case of Sm the RECP also replaces the electrons in the partially filled 4f 5 shell which do not participate actively in the bonding. A contracted [5s 4p 3d] valence basis set is employed for Sm. The starting basis set for the CH(SiMe3)2 ligands is the 6-31G basis. The effects of d functions on C and Si (6-31G*), p functions on the hydrogen atoms, and polarization f functions on the Sm have also been explored.

For one set of calculations a "small core" RECP was used on Sm in which the 4f 5 shell was explicitly treated as valence electrons. 11 This ECP also includes the core-like 4s 2 4p 6 4d 10 shells of the same principal quantum number, along with the outer electrons already treated in the "large core" RECP. A contracted [5s 5p 4d 3f] Gaussian basis was used for Sm in this case.

The geometry optimizations were performed without symmetry constraints (C1 symmetry). The resultant structures correspond to stationary points with no net force on each nucleus, but second derivatives with respect to nuclear distortions were not computed. All calculations were carried out using Gaussian98. 12 In the later discussion on van der Waals radii, standard radii 13 were taken for C (1.70 Å)

and Si (2.10 Å ). For the lanthanides values were obtained from covalent radii in the alkyl complexes using the relation 14 RvdW = Rcov + 0.80 Å, the La-C (2.51 Å ) and Sm-C (2.33 Å) bond lengths in the structures of this paper, and Rcov = 0.77 Å for carbon. This yields van der Waals radii of 2.54 Å (La) and 2.36 Å (Sm).

Results and Discussion

In discussing the results of the calculations it is useful to consider the experimental The geometries were optimized with use of the B3LYP functional and various basis sets, the most thorough study being carried out for the Sm model system. The results pertaining to 1 and 2 are reported in Table 1 for Sm and in Table 2 for La. For each metal system, the relevant experimental parameters from the X-ray study are also reported in the corresponding table. The experimental crystallographic determinations are of sufficient quality to unambiguously establish the molecular geometry. The La structure has been determined more accurately and provides relatively reliable distance and angles to non-hydrogen atoms. The Sm structure, on the other hand, is of much lower quality (vide infra). The hydrogen atoms, however, were not located nor refined freely for either structure. The optimized geometries for the R = Me and H systems with R = CH3, shown in Figure 1 (a) and (b), are qualitatively the same for both metals while some important differences are observed for the simpler model compound 3 (with all H substituents), see Figure 1 (c). The geometric parameters for the latter model compounds are collected in Table 3.

<Table 1, Table 2 and Table 3; Figure 1> (a) Calculations on the full molecules Ln[CH(SiMe3)2]3

The geometry of the full molecule 1 has been optimized only with the simplest basis set, which utilizes a large core Sm or La RECP (see Computational Details) and 6-31G bases for the Si, C and H atoms. Although no polarization functions were used for any of the atoms, the calculated geometries reproduce quite well the experimental geometry. The most peculiar molecular feature is the pyramidal arrangement of the central LnC3 core. The computed C-Ln-C angles are identical to the experimental ones and are very close to the ideal value for a tetrahedral geometry. One Si(Me)2-CH3 arm of each alkyl ligand (facing the metal from the apex side of the LnC3 pyramid and containing the silicon atom labeled as Si in Figure 1), finds itself in relatively close proximity to the metal center, suggesting the existence of electronic (agostic)

interactions. The other Si(Me)2-CH3 arm (with the silicon atom labeled as Si' in Figure 1), on the other hand, is farther away from the metal. In spite of the sterically based rationalization given in the original report, 1 we feel that the molecule should prefer a trigonal planar structure on the basis of purely steric arguments.

All bond distances are slightly longer than the experimental ones, as usually found at this level of theory, while all bond angles are reproduced to within 3°. A notable exception is the CH-Si-CH3 angle which is calculated 4.5° greater than the experimental value for the Sm structure, while experimental and computed values for the same parameter are in much better agreement in the La structure. The poorer quality of the Sm structure leads us to believe that the experimental value of this parameter might have been determined with a lower precision than suggested by the reported esd's. In particular, we suspect that the position of the metal-bound C atom is erroneous for the Sm structure. In fact, other parameters related to this atom are also badly reproduced by the calculation, while the same parameters are correctly reproduced for the corresponding La structure. This is the case for the Ln-C distance and for the C-Si and C-Si' distances. It is particularly notable that the latter two distances are inexplicably quite different from each other (beyond the reported esd's) in the experimental Sm structure, while they are essentially identical in the La structure, as well as in all our computational results (including those on the simpler model compound at all levels of theory, vide infra). Because of the close distance, the presence of Ln•••H interactions may also be questioned, (see Scheme 1). Although the accurate position of the H atoms is not available from the X-ray studies, the computational results help the analysis of this interaction. Two indicators are strongly suggesting that the Ln•••(H-C) interaction is rather repulsive in nature. The first one is the nearly staggered orientation, with respect to the Ln•••Si axis, adopted by the CH3 fragments closest to the metal center. The dihedral Ln•••Si-C-H angles are in the proximity of ±60° for all calculations on both metal systems (see Table 1 andTable 2), while an attractive interaction would be expected to place a single H atom as close as possible to the metal center, The data in Table 1 and Table 2 and a comparison of parts (a) and (b) of Figure 1 illustrate the validity of the chosen model to represent the essential metal-ligand interactions. All optimized parameters obtained at the 6-31G level are close to those obtained for the real molecule at the same level of theory. The largest differences can be attributed to the steric effect of the bulkier SiMe3 groups in the real molecule, as clearly illustrated by the much smaller Ln-C-Si' and Si-C-Si' angles in the model compound.

The presence of agostic interactions between the

The effect of an increase in the quality of the basis set was explored for the model system and only for the Sm system. As shown in Table 1, No significant changes for the Sm-ligand interactions ( bonding and agostic) result from the addition of the f functions, thus these interactions can adequately be described by using only the metal s, p and d orbitals.

A final calculation was also carried out with the "small core" RECP and basis for the Sm atom in which the unpaired electrons in the 4f 5 shell were explicitly treated along with all other electrons outside the [Ar]3d 10 core (see Computational Details), but without polarization functions on the H, C and Si atoms. With respect to the large core calculation without polarization functions, the major noticeable difference is a slight shortening of all bond distances related to the Sm atom, whereas all other distances display more marginal changes. This result is analogous to that obtained by Maron and Eisenstein on the simple Ln(NH2)3 model systems. 6 Finally we note that a second local minimum for the structure is also obtained for the model complex using the large core RECP. Starting from a different initial geometry, the molecule optimized to a new conformer (Figure 2) that is characterized by agostic interactions with only two of the alkyl substituents, while the third one adopts a staggered orientation of the As shown in Figure 1 (c) and Table 3, the replacement of the residual CH3 groups with H atoms in the model system discussed in the previous section maintains the pyramidal structure of the central LnC3 core. The degree of pyramidalization is quite comparable with those of the two systems discussed above and with those of the experimental structures. This new modification, however, has the more profound effect of replacing the agostic -Si-C interaction with a -Si-H interaction. In this respect, the present model is no longer chemically meaningful. 3).

(d) Discussion of the electronic structure.

In the preceding sections the discussion has focused on the structural results of the calculations. In this section some aspects of the electronic structure of the complexes will be discussed briefly with particular attention to Sm[CH(SiMeR2)(SiR3)]3 molecules where R = Me and H. For a neutral Sm atom the ground state has the atomic configuration …[4f 5 ] 5s 2 5p 6 6s 2 5d 1 . In the present calculations the electrons in the partially occupied 4f shell denoted by brackets are replaced by the RECP, and the filled outer core 5s and 5p electrons are explicitly treated in the calculations. These shells will be omitted in the following discussion. From a standard Mulliken population analysis of the molecular DFT wavefunction in Sm[CH(SiMeR2)(SiR3)]3 using the 6-31G basis, one obtains the following atomic populations for Sm: R = Me 6s 0.42 6p 0.53 5d 1.50 R = H 6s 0.41 6p 0.31 5d 1.32 leading to overall atomic charges of +0.54 (R = Me) and +0.84 (R = H) for Sm and hence overall charges of -0.18 and -0.28 on each of the three alkyl ligands, respectively. In addition to observing that this analysis portrays the Sm as much less ionic than the formal Sm(III) valence state would indicate, the other notable feature is the significant amount of 5d admixture in the molecular orbitals. Analogous results are obtained for La complexes where a slight (0.1 electron) increase in d population is obtained.

<Figure 3>

This 5d participation in the bonding is particularly evident in the highest occupied molecular orbitals. In the B3LYP calculations the highest 3 occupied orbitals represent the Sm-C bonding orbitals. The lower, more strongly bound orbitals are describing the Si-C-H bonds in the ligands. The set of 3 highest occupied orbitals is comprised of a degenerate set of MOs (e symmetry) and another MO (a symmetry) at slightly lower energy. Contour plots of these orbitals are shown in Figure 3 for the case Sm[CH(SiMeH2)(SiH3)]3. For the e set a decomposition of the electron density in each MO gives 17 percent Sm(5d) character and 83 % ligand character. For the a orbital the relative populations are 10 % Sm(6s) and 88 % ligand character. The two highest orbitals provide 0.68 e -of the total 1.32 e -of 5d population on Sm.

In the figure the e MOs, which are predominantly s-p hybrid orbitals on the alpha carbon atoms, interact with the 5d orbitals on Sm. For the a MO the symmetric combination of similar orbitals on the alkyl ligands interacts with the 6s orbital on Sm, where the nodal structure in the region of the 5s is also evident.

The lanthanide trihalides 15 present a somewhat analogous situation to the tris-alkyl complexes in that there is a very soft mode corresponding to motion between planar and pyramidal structures. The results from calculations [16][17][18][19] show that many of the LnF3 molecules are pyramidal. The calculations also show a greater tendency for heavier trihalides to be planar, but there is not unanimity in the predictions. For example, DFT calculations 19 After submission of this manuscript, we have learned about an independent study by Perrin, Maron, Eisenstein and Lappert 20 on the same system. This study describes the bonding, including to the -Si-C bond, as essentially purely electrostatic on the basis of a natural bond orbital (NBO) analysis. 21 The same study, however, highlights the importance of the lanthanide d orbitals in enforcing the pyramidal configuration at the metal center through a 2 nd order Jahn-Teller effect.

Conclusions

The

  Ln center and the Si(CH3)2-CH3 arm is strongly suggested by several parameters. The most notable ones are the Ln-C-Si angles. These are much smaller (<La-C-Si = 104˚ (calcd), 102˚ (exptl)) than the Ln-C-Si' angles (<La-C-Si' = 121˚ (calcd), 121˚ (exptl)) and smaller than the ideal tetrahedral value. The Ln-C-Si' angles are comparatively larger than the tetrahedral value. A second parameter is the Si-CH3(•••Ln) distance (namely, the distance to the methyl group which is closest to the lanthanide metal center). This distance is significantly longer (0.033 -0.035 Å calc., 0.05 -0.09 Å exptl.) than all other Si-CH3 distances. These two features are clear indications of an electronic interaction between the metal center and the -SiC (i.e., the SiCbond. It is notable that the angular distortion is strongest at the C atom and weakest at the Si atom, as might be expected for an agostic interaction. On the other hand, the bond distance trends analyzed above clearly indicate that the Sm atom interacts with the -SiC bond (which is lengthened) and not with the -CSi bond. Additional indications for the presence of -SiC agostic interactions are the relatively short Ln•••Si and Ln•••C contacts that are significantly shorter than the sums of the van der Waals radii. For example, using an estimated radius of 2.5 Å for La (see Computational Details), the sums of van der Waals radii for La•••Si (4.6 Å) and La•••C (4.2 Å) interactions are significantly larger than the observed La•••Si 3.4 Å and La•••C 3.1 Å interatomic distances. With a van der Waals radius of 2.33 Å for Sm, the same analysis shows the observed Sm•••Si and Sm•••C contacts of 3.3 and 3.0 Å would also be less than the sums of van der Waals radii.

  with a dihedral angle close to zero. The Ln-C-Si-C dihedral angle, on the other hand, is close to zero, in agreement with the presence of the -Si-C interaction. The second indicator is the tilting of the Me group away from the metal, as shown by the greater than tetrahedral value for the Si-C-H angles relating to the two closest H atoms, and smaller value for the third angle. An analogous situation is evident from the DFT results for the compound (C5H5)La[CH(SiMe3)2]2.3 It is interesting to observe that the C-H distances for the two C-H bonds closest to the Ln center are marginally longer than the third C-H bond. While the difference is quite small, this phenomenon is consistently found for all Me groups within the same molecule and also for the model system at all levels of theory (vide infra). This trend is consistent with the carbon rehybridization caused by the distortion. In fact, the C-H bonds with a greater Si-C-H angle have a greater C p character, leading to bond lengthening, while those with a smaller Si-C-H angle have a greater C s character, leading to bond shortening. Therefore, there are no indications whatsoever for the presence of -C-H agostic interactions in these molecules, while the metal center establishes -Si-C interactions with all three alkyl ligands. Once again, an analogous situation exists for compound (C5H5)La[CH(SiMe3)2]2. 3(b) Calculations on the Ln[CH(SiH2Me)(Si'H3)]3 model system.

  additional calculations were carried out by adding polarization functions on the CH3 hydrogen atoms (p), on all silicon atoms (d), on the C atoms (d, plus the two above polarization functions on H and Si), and on the Sm atom (f). The energetic gain and the structural modifications are substantial only when d functions are introduced on the Si and C atoms, see Table 1, while they are very minor when p functions are added to the H atoms and especially when f functions are added to the Sm atom. The added polarization functions do not qualitatively affect the nature of the agostic interactions, although the Sm•••Si and Sm•••C contacts do get marginally shorter, especially upon adding d functions on the Si atoms. The C-Si distances also shorten significantly upon addition of the Si d function.

SiH2Me group relative to

  the Sm-C bond. The loss of one agostic interaction with a Si-C bond is compensated by the strengthening of the other two Sm•••(Si-C) interactions. In fact, the two Sm•••Si and Sm•••C distances become 3.409 and 3.183 Å on average, respectively, while those to the non-interacting alkyl group are lengthened to 3.57 and 4.03 Å. The total energy of this optimized geometry is slightly higher (0.4 kcal/mol) relative to the triply agostic conformer in the 6-31G basis with p polarization functions on H. This small difference suggests a facile pathway for the exchange of the interacting Si-C bonds in the parent compound by rotation about the C-Si bond and is consistent with the observation of single resonances in the 1 H and 13 C NMR spectra of compound La[CH(SiMe3)2]3 for all SiMe3 groups. 1 Obviously, a second process involving the inversion of the LaC3 pyramid or rotation about the Sm-C bond must also be invoked to account for the experimental result. <Figure 2> (c) Calculations on the Ln[CH(SiH3)2]3 model system.

  The -Si-H agostic interactions in this simpler model are stronger than their -Si-C counterparts discussed above. This is clearly indicated by the shorter Ln•••Si contacts (by over 0.2 Å) and by the smaller Ln-C-Si angles (by 10°). Additional indications of the presence of strong agostic interactions are the close Ln•••H contacts and the significant lengthening of the interacting Si-H bond relative to the other non-interacting bonds (see Table

  on LnCl3 molecules give pyramidal structures (but with larger bond angles than trifluorides), while other SCF and MP2 calculations show trichlorides to be planar. Analysis of the orbitals from the DFT calculations on LnX3 species for X=F, Cl, Br and I showed as much as 15-20 percent d orbital character for the cases X= Cl and I. From the experimental and theoretical studies one could conclude that the most ionic LnF3 molecules show the greatest tendency for nonplanar structures.Additional information regarding possible agostic interactions can be obtained from bond overlaps from the Mulliken population analysis. In the real molecule (R=Me) the Sm-C bond overlap is 0.075 for the agostic methyl groups compared to 0.238 for the Sm-C bonds and 0.01-0.02 for all other Sm-C interactions. The Sm-H bond overlaps for the two closest H atoms in each agostic methyl group are 0.026 and 0.011, several times the magnitude of other Sm-H interactions.

  present study has shown that the unusual structural features of compounds Ln[CH(SiMe3)2]3 are associated with -Si-C agostic interactions established by the metal center with all three alkyl ligands rather than -C-H agostic interactions. The interaction between the lanthanide metal and the -C-H bonds are in fact repulsive. The essential features of these structures can therefore be correctly reproduced by using the Ln[CH(SiH2Me)(Si'H3)]3 model system. The replacement of the last Me group in each alkyl ligand with a H atom, while maintaining the basic experimental geometry with a triple -agostic interaction, no longer provides a faithful chemical representation of the real molecule. The use of several ameliorated basis sets with the addition of polarization functions on the H, C, Si and Sm atoms does not allow a significantly better description of the molecular features. Finally, the excellent match of the optimized data with the experimentally better determined La structure and the poorer match with the lower quality Sm structure shows an example of the utility of computational chemistry for the assessment of experimental data.

  Figure 1. Views of the optimized geometries for compounds Ln[CH(SiR2R')(Si'R3)]3 (Ln =

Figure 2 .

 2 Figure 2. A top view of the Sm[CH(SiH2CH3)(Si'H3)]3 conformer with only two Sm-alkyl

Figure 3 .

 3 Figure 3. Contour plots of the three highest occupied orbitals of Sm[CH(SiH2CH3)(SiH3)]3

  Figure 1

Table 1 .

 1 Comparison of optimized and experimental structural parameters for

	Sm[CH(SiR2Me)(Si'R3)]3. a
		R = Me				R = H
		6-31G	exp	6-31G	6-31G	6-31G	6-31G	6-31G	6-31G
					(H**)	(Si*)	(H**C*Si*)	Sm(fpol )	Sm SC
	Sm-C	2.459	2.33(2)	2.449	2.450	2.457	2.460	2.439	2.420
	Sm•••Si	3.49	3.326(7)	3.46	3.45	3.43	3.43	3.45	3.40
	Sm•••C	3.24	3.04(2)	3.32	3.29	3.28	3.27	3.31	3.24
	Sm•••H	2.96	-	3.01	2.99	2.95	2.92	3.02	2.95
		3.15	-	3.20	3.23	3.29	3.25	3.24	3.18
	C-Si	1.906	1.80(3)	1.900	1.898	1.861	1.860	1.899	1.896
	C-Si'	1.906	1.94(3)	1.900	1.899	1.864	1.864	1.899	1.894
	Si-CH3(•••Sm)	1.956	1.95(3)	1.950	1.954	1.921	1.923	1.951	1.951
	Si-CH3(other)	1.923	1.86(5)	-	-	-	-	-	-
	C-H(•••Sm)	1.103	-	1.103	1.100	1.105	1.103	1.102	1.102
		1.099	-	1.099	1.097	1.100	1.098	1.099	1.099
	C-H(other)	1.097	-	1.096	1.094	1.097	1.096	1.096	1.096
	C-Sm-C	110.2	110.3(10)	108.4	108.4	108.2	108.4	108.7	107.6
	Sm-C-Si	103.9	106.6(12)	104.6	104.4	104.2	104.4	104.6	103.3
	Si-C-H(other)	108.0	-	108.6	108.7	108.7	109.6	108.5	108.4
	Sm-C-Si-C	14.0	14.4(17)	15.2	13.9	18.1	17.1	15.4	15.6
	Sm•••Si-C-H	-57.3	-	-54.8	-55.2	-51.9	-52.0	-55.7	-55.1
		68.8	-	70.7	70.2	73.5	73.0	69.8	70.4
	Rel. energy			0	-0.01962 -0.21630 -0.24904 -0.00451
	(hartrees)

a Values are given as averages of geometrically equivalent structural parameters.

Table 2 .

 2 Comparison of optimized and experimental structural parameters for La[CH(SiR2Me)(Si'R3)]3. a

		R = Me	R = H
		6-31G	exp	6-31G
				(H**C*Si*)
	La-C	2.554	2.516(10)	2.552
	La•••Si	3.54	3.411(3)	3.52
	La•••C	3.31	3.12	3.34
	La•••H	3.02	-	3.01
		3.20	-	3.29
	C-Si	1.903	1.841(10)	1.860
	C-Si'	1.903	1.837(11)	1.861
	Si-CH3(•••La)	1.958	1.923(10)	1.925
	Si-CH3(other)	1.923	1.875(15)	-
	C-H(•••La)	1.103	-	1.102
		1.099	-	1.099
	C-H(other)	1.097	-	1.095
	C-La-C	111.5	109.9(4)	108.9
	La-C-Si	104.1	101.9(4)	105.0
	La-C-Si'	120.6	121.0(6)	116.6
	Si-C-Si'	118.3	121.6(6)	114.8
	CH-Si-CH3	110.3	109.7(4)	111.8
	CH-Si-R	113.6	113.4(6)	111.5
		114.7	115.2(5)	114.6
	Si-C-H(•••La)	114.2	-	114.1
		113.0	-	112.6
	Si-C-H(other)	107.8	-	109.1
	La-C-Si-C	13.4	13.9(7)	13.6
	La•••Si-C-H	57.0	-	53.9
		-68.9	-	-71.3

a Values are given as averages of geometrically equivalent structural parameters.

Acknowledgements. We acknowledge support from Laboratory Directed Research and Development at Los Alamos. Los Alamos National Laboratory is operated by the University of California for the U.S. Department of Energy under Contract W-7405-ENG-36. R.P. is grateful to CINES (Montpellier, France) for granting free computer time and the G. T. Seaborg Institute for Transactinium Science for sponsoring a stay at Los Alamos National Laboratory. We thank Odile Eisenstein for sending us a preprint of her parallel study. 20


Masaru Nagaso
email: nagaso@lma.cnrs-mrs.fr

Joseph Moysan

Christian Lhuillier

Dimitri Komatitsch

Four dimensional wave propagation simulation in liquid metal with realistic thermal fluctuation

Keywords: Ultrasound, Fluctuating liquid media, Non-Destructive Testing (NDT), numerical modeling, spectral finite elements (SEM), Large-Eddy Simulation (LES)

In this paper, we performed a four-dimensional numerical simulation to calculate wave propagation in a thermal fluctuating liquid metal in order to obtain detailed knowledge on a wave propagation in coolant material of a Sodiumcooled Fast reactor (SFR). The wave and the medium are described in three spatial and one temporal dimensions. We made use of a massive data set to

describe the fluctuating situation of the medium. This data set was provided by Computational Fluid Dynamics (CFD) with a Large-Eddy Simulation (LES) model, which calculated the temperature field with a higher spatial resolution than Reynolds-Averaged Navier-Stokes turbulence models (RANS). This data set was furthermore obtained from other studies on a numerical and physical experiment called PLAJEST that created mixing jets of liquid metal in order to simulate the status of running SFRs. Because of the limitation of computational resources, previous acoustic studies applied to such a medium could only use the spatial-temporally averaged fluctuating heterogeneity of a medium calculated by RANS turbulence model. This limitation may overlook wave fluctuation because of the difference of the resolution between computational fluid dynamics and acoustic wave length. Our numerical effort allowed us to study the most realistic acoustic wave propagation in liquid metal than in any former studies. A new

Introduction

The fourth generation of Sodium-cooled Fast Reactors (SFR) is actively developed within a framework of an international cooperative research project [START_REF] Gif | Framework agreement for international collaboration on research and development of Generation-IV nuclear energy system[END_REF][START_REF] Gif | Technology roadmap update for Generation IV Nuclear Energy Systems[END_REF]. In France, a research and development project of this fourth generation SFR, so called ASTRID project, is now being carried out [START_REF] Hora | Determination of Lamb wave dispersion curves by means of Fourier transform[END_REF].

As SFRs use liquid sodium for their coolant material, which is completely opaque, the method to investigate the running state of SFRs is quite limited. Also because of the high reactivity of liquid sodium to air, it is important to reduce the number of mechanical parts inside of SFR design in order to reduce the risk of sodium leakage and the difficulty of maintenance. For those requirements, measurement methods using acoustic and elastic waves (i.e. non-destructive testing) are regarded as techniques which have advantages for monitoring SFRs [START_REF] Baqué | Development of tools, instrumentation and codes for improving periodic examination and repair of SFRs[END_REF][START_REF] Baqué | ASTRID in service inspection and repair: Review of r&d program and associated results[END_REF].

Before now, an acoustic method has been applied for the monitoring of the position of subassembly heads. The so-called SONAR monitoring system was developed by the Alternative Energies and Atomic Energy Commission (CEA) [START_REF] Berton | Continuous monitoring of the position of two subassemblies heads of Phenix at 350 MWth power and 550 deg C temperature[END_REF]. Using a standard emitter-receiver setup, the SONAR system measures the distance between a target object and an acoustic probe thanks to the conversion of the round-trip time of flight. Several studies developed theoretical or experimental approaches to improve ultrasonic telemetry. Wedges of objects could be used to create reflection points of interest. Lü et al. compared several high-frequency approximations by simulating acoustic scattering by an immersed rigid wedge.

They confirmed that for a planar wedge, these approximations produce the same predictions [START_REF] Lu | Models comparison for the scattering of an acoustic wave on immersed targets[END_REF][START_REF] Lü | Numerical comparison of acoustic wedge models, with application to ultrasonic telemetry[END_REF]. They concluded that Physical Theory of Diffraction (PTD) [START_REF] Ufimtsev | Fundamentals of the physical theory of diffraction[END_REF] is more efficient as it requires integration only over the illuminated part of the scattering surface and they recommended it for implementation in modelling software. From a more practical issue, the modelling of such signals and its comparison with experiments is a key point to help designers to define better geometries for monitoring situations [START_REF] Ploix | Investigation of ultrasonic backward energy from various edges as a function of their 2d/3d geometry and of the incidence angle, for application to ultrasonic thermometry at the outlet of a tube[END_REF].

Acoustic method is also regarded as a powerful tool for temperature measurement of sodium coolant, because of its higher measurement frequency than thermo-couples [START_REF] Mcknight | Remote temperature measurement[END_REF]. Thus, acoustic measurement techniques have a great advantage for application to SFR monitoring but we still need to advance in the knowledge on wave propagation in SFRs to increase the measurement accuracy and reliability, and to expand the application area of acoustic measurements. The need for advanced studies on wave propagation in SFRs comes from the complexity of the propagation medium status as acoustic characteristics of liquid sodium vary dynamically and in a complex manner during an SFR's operation. The aim of this study is to develop a new methodology for this monitoring that could be effective in any other liquid metals.

In the sodium coolant, the three fields, flow velocity field, temperature gradient field and gas bubble density field, may affect a wave propagation. In this study, we only focus on temperature gradient, as the effect from the flow velocity field is weaker than the effect of the temperature gradient field, in the case of an SFR's operation situation [START_REF] Massacret | Simplified modeling of liquid sodium medium with temperature and velocity gradient using real thermal-hydraulic data: Application to ultrasonic thermometry in Sodium Fast Reactor[END_REF], and the temperature gradient may cause an undesirable inflection [START_REF] Massacret | étude d'une méthode ultrasonore d'estimation des températures locales du sodium liquide en sortie coeur RNR-Na[END_REF][START_REF] Massacret | Modelling of ultrasonic propagation in turbulent liquid sodium with temperature gradient[END_REF]. The effect from the density of gas bubbles has been studied in [START_REF] Cavaro | Microbubble cloud characterization by nonlinear frequency mixing[END_REF][START_REF] D'hondt | Acoustical characterisation and monitoring of microbubble clouds[END_REF]. In our study, we considered that the density of gas bubbles is constant. The result of a recent study on the thermal-hydraulic state of SFRs [START_REF] Haubensack | Calculs thermohydrauliques du collecteur chaud de la chaudière ASTRID au nominal en soutien l'instrumentation coeur CFV-BD-16/10 et ISIR[END_REF] shows that the maximum variation rate of flow rate is about 4.5 m s -1 and the maximum variation rate of temperature may be about 23.4 • C par 10 cm. For SONAR application, the propagation distance is around 270 mm. Ultrasonic telemetry in sodium for distance from 1 m to 5 m are sought or tested in various reactors [START_REF] Tarpara | Under-sodium ultrasonic viewing for fast breeder reactors: A review[END_REF], [START_REF] Chaitanya | Detection of subassembly bowing in a fast breeder reactor using glancing angle ultrasonic imaging[END_REF]. Thus this variation may be not negligible for the possible application of acoustic measurement in SFRs.

When considering the effect of temperature gradient, we need to take into account liquid sodium's temperature-dependent properties reported in [START_REF] Sobolev | Database of thermophysical properties of liquid metal coolants for GEN-IV[END_REF]. This study shows that the density difference of liquid sodium with temperature change between the temperature range from the normal melting point to the normal boiling point may be calculated with the linear relation:

ρ [kg m -3 ] = 1014 -0.235 • T [kelvin], (1)
at normal atmospheric pressure, where ρ is the density of sodium and T is the sodium temperature in Kelvin degrees. The sound speed of liquid sodium decreases monotonically with temperature, which is caused by the decreasing of the inter-atomic interactions. In the range of the normal melting -boiling point (371-1155°K), the sound speed in pure liquid sodium may be described based on a linear relation:

c p [m s -1 ] = 2723 -0.531 • T [kelvin]. (2)
where c p is the celerity of ultrasonic waves in meters per second.

We have carried out several numerical studies on the wave propagation in order to advance our simulation method for liquid sodium. For an SFR's study, numerical calculation have great importance as the real size experiment of SFR with actual liquid sodium is not possible, as there is no operating SFR for an experimental study. Thus a realistic modeling of wave propagation is essential for designing an acoustic measurement system for SFR application.

In the former studies which we have conducted for researching a sodium temperature heterogeneity effect on a wave propagation, we have improved the modeling method of this heterogeneity. In the early stage, we applied stochastic methods which randomly generate a fluctuating temperature field using Gaussian Random process in [START_REF] Fiorina | Application of the Gaussian beam summation method to the study of the ultrasonic wave propagation in a turbulent medium[END_REF][START_REF] Iooss | Statistical moments of travel times at second order in isotropic and anisotropic random media[END_REF][START_REF] Iooss | Numerical simulation of transit-time ultrasonic flowmeters: uncertainties due to flow profile and fluid turbulence[END_REF][START_REF] Lü | Stochastic simulation of the high-frequency wave propagation in a random medium[END_REF][START_REF] Nagaso | Ultrasonic thermometry simulation in a random fluctuating medium: Evidence of the acoustic signature of a one-percent temperature difference[END_REF]. This method has a strong advantage of not consuming a large amount of computer memory and time when generating a fluctuation with high resolution. However as remarked in [START_REF] Iooss | Numerical simulation of transit-time ultrasonic flowmeters: uncertainties due to flow profile and fluid turbulence[END_REF], generated fields may not well model the geometrical heterogeneity, e.g. a mean flow field occurring due to an interaction with internal structures of SFR.

N. Massacret [START_REF] Massacret | étude d'une méthode ultrasonore d'estimation des températures locales du sodium liquide en sortie coeur RNR-Na[END_REF] applied a Computational Fluid Dynamics (CFD) calculation result with a Reynolds-Averaged Navier-Stokes (RANS) turbulent model. A RANS model requires relatively lower quantity of computational resources than others, e.g. Direct Numerical Simulation (DNS) and Large-Eddy Simulation (LES), by filtering the small turbulences. Thus the spatial resolution is not enough when compared with the acoustic wavelength, and not suitable for verifying the acoustic effect from smaller fluctuation pattern.

In this study, we carried out three-dimensional numerical simulations with application to a more realistic fluctuating propagation medium. In the SFRs, the thermal-hydraulic situation is generated by the sodium jets with a high temperature and surrounding sodium with a lower temperature, and mixing phenomena between them. As the modeling target of our study, we selected an experimental and numerical study called PLAJEST, since it targets the same object in the same condition, i.e. the upper-core region of an SFR in operation.

PLAJEST is a collaboration between the Japanese Atomic Energy Agency (JAEA), the U.S. Department of Energy (DOE), and by the Alternative Energies and Atomic Energy Commission (CEA). 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 thus may cause cracks in adjoining structures of SFRs. The configuration of this experiment is shown in Figure 1 (a). After the experiment made by JAEA [START_REF] Kimura | Study on convective mixing phenomena in parallel triple-jet along wall -Comparison of temperature fluctuation characteristics between sodium and water[END_REF][START_REF] Kimura | Study on thermal striping phenomena in triple-parallel jet -Transfer characteristics of temperature fluctuation in sodium and water based on conjugated numerical simulation[END_REF], a thermal-hydraulic analysis was carried out by the CEA in its Service de Thermo-hydraulique et de Mécanique des Fluides (STMF) ( [START_REF] Angeli | Large-eddy simulation of thermal striping in WAJECO and PLAJEST experiments with TrioCFD[END_REF], Figure 1 (b)). In that analysis, the results indicate a good fit on time-averaged normalized temperature and its fluctuation between experimental and numerical values.

We consider these CFD thermal-hydraulic computation results are well validated by the experiment and that allows us to process full-wave propagation numerical simulations in realistic conditions. CEA/STMF used a Large-Eddy Simulation turbulence model for this simulation, thus we may utilize these temperature fields with higher spatio-temporal resolution. In the first part of this study we analyse the temperature field issued from the mixing jets and we

show that it could be characterized by two representative points of the merging of jets. VisIt, an open source visualization tool for massive scientific data [START_REF] Childs | High Performance Visualization: Enabling Extreme-Scale Scientific Insight[END_REF].

In the second section, we firstly introduce the numerical method applied for wave propagation modeling and the technical details used for preparing the finite element mesh and the pre-processings for the CFD result data set.

Then in the third section, we explain the method for analyzing the fluctuation of the thermal field and the acoustic field calculated with our model. In the fourth section, we show the details of the calculation setups used for our wave propagation calculations. We also show an example of the calculated 3D wave field propagating at one single moment of the fluctuating temperature field. In the fifth section, we detail the analysis method based on second derivative properties of several wave parameters. We demonstrate that ultrasonic measurements in a through transmission configuration can find the two characteristic points of the temperature field. It opens the way to prospect industrial design to monitor such thermal-hydraulic situations using ultrasonic transducers.

Numerical method

Spectral element method

Our numerical wavefield modeling method is a time domain spectral-element method [START_REF] Patera | A spectral element method for fluid dynamics: Laminar flow in a channel expansion[END_REF], which is a type of finite-element method using higher order basis functions (generally higher than 4th order) than an ordinal finite-element method using 1st or 2nd order basis functions.

We have chosen this method because this may model the heterogeneity of physical parameters, i.e. temperature and density values in our case, with higher resolution and lower number of elements, thus we may efficiently model the spacial domain of wave propagation. This method directly increases the limit of domain size and the frequency of modeled waves may becomes higher than the ordinal finite-element methods.

Following our previous studies [START_REF] Nagaso | Ultrasonic thermometry simulation in a random fluctuating medium: Evidence of the acoustic signature of a one-percent temperature difference[END_REF][START_REF] Nagaso | Wave propagation simulation in the upper core of sodium-cooled fast reactors using a spectralelement method for heterogeneous media[END_REF], we use the software SPECFEM3D [START_REF] Komatitsch | Méthodes spectrales et éléments spectraux pour l'équation de l'élastodynamique 2D et 3D en milieu hétérogène (Spectral and spectralelement methods for the 2D and 3D elastodynamics equations in heterogeneous media)[END_REF] which is the one of the most well utilized softwares for full-wave wave propagation simulation because of its very matured and advanced development status especially for the use on HPC with the great numbers of CPUs [START_REF] Komatitsch | A 14.6 billion degrees of freedom, 5 teraflops, 2.5 terabyte earthquake simulation on the Earth Simulator[END_REF] and GPUs [START_REF] Komatitsch | Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA[END_REF][START_REF] Komatitsch | High-order finiteelement seismic wave propagation modeling with MPI on a large GPU cluster[END_REF] as well.

Mesh generation and interpolation of the temperature field

The CEA/STMF department used a numerical code for CFD calculations For their calculation, a tetrahedral unstructured staggered mesh was used.

Temperature field values are defined at the center of each TrioCFD's tetrahedral mesh element, and flow velocity values are defined on the vertex nodes. We thus had to transfer these values to our hexahedral mesh for SPECFEM3D. To do so, we used interpolation onto each node of the SPECFEM3D hexahedral mesh using the simulation data management tool called MEDCoupling. MEDCoupling is part of the pre-/post-processing platform SALOME (http://www.salome-platform.

org) and is also available independently as a library. Figure 2 shows the temperature field data transfer and mesh generation steps as a pre-process for SPECFEM3D simulation. Our hexahedral mesh for SPECFEM3D is built using the meshing software CUBIT developed by Sandia National Laboratories (USA).

It is possible to select an arbitrary volume to be extracted from the entire geometry and meshing is completed automatically, including the assignment of material characteristics and absorbing surface flags (the fluctuation of the temperature field may be defined later). After finishing preparation of mesh data, we carried out the temperature field transfer, i.e. interpolation of temperature values defined at the barycenter of each tetrahedral finite element to corner nodes of our hexahedral spectral elements. The flow velocity data is not used for our simulation because we apply the frozen fluid hypothesis, as in [START_REF] Massacret | Simplified modeling of liquid sodium medium with temperature and velocity gradient using real thermal-hydraulic data: Application to ultrasonic thermometry in Sodium Fast Reactor[END_REF].

The conversion of the temperature field from TrioCFD tetrahedral mesh to SPECFEM3D nodes is done by using MEDCoupling.

Determination of the element size to use in our simulations is done based on two conditions, which are the CFL condition (Equation 3) and the number of elements per one wave length:

C p ∆t ∆x gll ≤ α , (3)
where ∆t is the time step and ∆x gll is the minimum interval between two Gauss-Lobatto-Legendre (GLL) grid points. GLL points are the additional collocation points inserted between the corner nodes of a general finite element. The position of those points are defined by numerically solving the equation:

(1 -x 2 )P n (x) = 0, (4)
where P n denotes the derivative of the n-th order Legendre polynomial (i.e.

d dx P n (x)).
Then the recursion relation for these polynomials is

P n+1 (x) = 2n + 1 n + 1 xP n (x) - n n + 1 P n-1 (x), P 0 (x) = 1, P 1 (x) = x. (5)
Detailed explanations about GLL polynomials can be found in [START_REF] Davis | Methods of numerical integration (Computer science and applied mathematics[END_REF][START_REF] Canuto | Spectral methods[END_REF].

We selected the averaged Courant number α = 0.4 and the wave celerity C p = 2416.268 m s -1 (in sodium with the lowest temperature value 274.5 • C in the CFD simulation) for the calculation of the mesh size and time step duration.

In equation 3, ∆x gll is not the mesh size itself, it is the interval between GLL grid points inside the spectral elements. This led us to use a mesh size ∆x = 8.05 × 10 -4 m and a time step of 2.3 × 10 -8 s. We simulate a total of 5000 steps in order to have a sufficient total physical duration for the waves to travel through the entire simulation domain. The mesh used in our simulations thus 3 250 000 spectral elements and a total number of GLL grid nodes of 215 320 764.

The figure 3 represents one of the heterogeneous temperature fields at one time step.

Medium fluctuation and acoustic fluctuation analysis 185

The phenomenon of mixing flows has been actively studied by e.g. [START_REF] Durve | Numerical investigation of mixing in parallel jets[END_REF], [START_REF] Zang | On the wake-like vortical arrangement and behaviour associated with twin jets in close proximity[END_REF] and [START_REF] Ghahremanian | Near-field mixing of jets issuing from an array of round nozzles[END_REF]. These studies are not directly related to the PLAJEST experiment nor to liquid metal flows, but numerous thermal-hydraulic studies show that they have common thermal-hydraulic behaviors [START_REF] Massacret | Modelling of ultrasonic propagation in turbulent liquid sodium with temperature gradient[END_REF]. We therefore use the methodology proposed in these studies to analyze thermal-hydraulics. They 190 identified the mixing state of flows and categorized them into three types based Durve and al. [START_REF] Durve | Numerical investigation of mixing in parallel jets[END_REF] also carried out a comparative study on several models for predicting the mean temperature field and temperature fluctuation field caused by mixing phenomena of the three jets [START_REF] Durve | Thermal striping in triple jet flow[END_REF].

Following different authors we choose to analyze the temperature in the medium using an index called the Temperature Fluctuation Intensity (TFI). We use the definition of TFI as

T F I(r) = 1 N N i=1 (T (i, r) -T (r)) 2 , (6)
where r = (x, y, z) is the spatial position vector, i is the time step number, N is the total number of time steps, T (i, r) is the temperature value at time step i and position r, and

T (r) = 1 N N i=1 T (i, r) is the average temperature at r.
We process the PLAJEST Computational Fluid Dynamics (CFD) data to calculate the TFI index in order to compare the global behavior of the jets with these previous studies. this TFI field, it is found that from the shape of the TFI profile depending on the altitude, it seems to be possible to define three zones in a similar way as [START_REF] Durve | Thermal striping in triple jet flow[END_REF]. First, the two high TFI zones arise around altitude z = 0.05 m and these two zones are completely separated. Around altitude z = 0.08 m or lower altitude, the beginning of merging of the high TFI zones becomes clear (the lowest TFI value at point b between the two peeks a and c starts to increase).

Here there would be a specific point that we call a merging point of the TFI zone. The merging of these two zones is confirmed when the altitude becomes higher, and then around altitude z = 0.18 m and z = 0.19 m these two zones are completely merged (the peaks a and c of the 1D TFI curve become a single one).

We call this altitude a combining point of the TFI zone. These two altitudes are estimated of course on this numerical data.

In the following section, we carry out a more detailed analysis to verify the variations of the TFI profiles and the definition of these two points. The Figure 5 shows the insonified volume (central altitude z = 0.1 m) and PLAJEST's geometry. The point B is the origin (point (0,0,0)) of the cartesian coordinate system we used. The numbers in meters show the distance from the source to each y-z receiver plane. The near field limit, called Io, is calculated to be 67.17 mm using the empirical equation 2, with an average temperature of 333.167 • C (606.317 • K), considering a transducer with a 1 inch diameter and 1 MHz frequency. In this virtual setup, only the first two observation surfaces are in the near field. Thus, in the following analysis we mainly have acoustics virtual measurements in the far field. In the PLAJEST coordinates the limit of the near field is x = -0.061 m. axis indicates time. The visual information reveals that the wave fronts having passed through heterogeneous liquid sodium are deflected and deformed but that the amount of wave deformation is not so large. Waveforms are different between waves passing from the end of the near field region to a far longer distance in the far field region (almost 4 × Io, where Io is the near field limit).

The wavefront in the near field has its shape that is more complex than in the far field, as expected by usual theory of wave field distribution. For this 3D visualization, we used a VTK file (Visualization Toolkit: an open-source, freely available software system for 3D computer graphics, image processing, and visualization. http://www.vtk.org) that we displayed with the visualization software VisIt.

Analysis of time-varying temperature fields of PLAJEST

Selection of the CFD time steps

The CFD calculation carried out by the CEA STMF has approximately 10 s in total, from 200.000 s to 210.197 with a 0.001 s interval. The initial 200 s of the calculation was dedicated to ensuring stabilization of the flow and are thus excluded in our analysis. Because of the limitation of allocated computation time on the supercomputer that we use for this work, it was not possible to run the temperature field interpolation and wave propagation calculation processes for all of these CFD time steps. Instead, we had to extract several time steps of the temperature field with a wider time interval from the CFD results. In order to select the time step interval for our acoustic simulation, we used the power spectrum density curve obtained in [START_REF] Angeli | Large-eddy simulation of thermal striping in WAJECO and PLAJEST experiments with TrioCFD[END_REF] which indicates the temperature history at x = -0.015 m (between the left and center jets), y = 0.09 m (middle point on the y axis) and z = 0.1 m. The peak is found around 3 Hz. Thus, to be sure to include the frequency of temperature fluctuation up to 5 Hz, we extracted the temperature fields with a 0.1 s interval (i.e. 10 Hz). In a first approximation the peak was estimated to be around 2 Hz, and in that case we could expect to have 5 points per period to keep the peak at 2 Hz. This is the limit of Shannon's sampling criterion.

Computing TFI data from a large spatial-temporal 3D data set and comparison with acoustic simulation results

As introduced in Equation 6, we resort to an index called Temperature Fluctuation Intensity (TFI), as used in [START_REF] Angeli | Large-eddy simulation of thermal striping in WAJECO and PLAJEST experiments with TrioCFD[END_REF]. Because of the very large numbers of total time steps and also the huge number of mesh nodes included in the CFD calculation results, the standard way to calculate the TFI value based on this equation is not very efficient. To calculate that TFI, we thus selected and implemented another, more advanced algorithm: the online algorithm, which we will briefly describe in this section. Figure 7 shows the calculated TFI field in 3D.We find that there are two areas where the TFI value becomes high between the sodium jets at altitude z = 0.08 m to 0.16 m. It should also be noted that the TFI field is not symmetric with respect to the x center.

The online algorithm calculates some field value from serial data sequentially and based on a single step [START_REF] Knuth | The Art of Computer Programming[END_REF]. This algorithm can be required for serial data for which each step needs a large amount of computer memory and/or when the number of serial data is so large that it is very expensive to perform an entire loop of calculations more than twice. We thus applied this algorithm to compute the TFI, i.e. the standard deviation of the temperature value at a given point, because computing a standard deviation implies several loop over the whole time steps, first to calculate the mean temperature, and second to calculate the difference between a temporal value and the mean value. The PLAJEST CFD data comprise 10,000 time steps of 3D volume data with 2,039,769 mesh elements. One entire loop calculation for this data takes about 10 hours. By applying the online algorithm for the calculation of the TFI, we only need to do this long loop calculation once.

The definition of standard deviation of temperature at a given position r is

σ 2 (r) = 1 N N i=1 (T (r, i) -T (r)) 2 , ( 7
)
where N is the total number of time steps and T (r) = 1 N N i=1 T (r, i) is the average temperature at position r. In the online algorithm, the averaged value

The standard deviation calculation is then finalized as

σ r = M r,N N . ( 10
)
Figure 8 shows the TFI values on the x axis of y = 0.09 m at several altitudes (i.e. z positions), and Figure 4 represents these TFI curves drawn in the case of the 2D TFI field. The maximum TFI value is found at around altitude z = 0.13 m. The change of shape of the curves depending on the distance from the exit of the jets seems to match with the result of [START_REF] Durve | Thermal striping in triple jet flow[END_REF], i.e. the curves at low altitude have two TFI peaks, and then these peaks gradually merge when the z altitude increases. We will further analyze this effect in the next sections of this article.

Calculating the "Cumulated TFI" (CTFI) value

The TFI calculated above is the index that evaluates the intensity of the fluctuation at one spatial point, while acoustic wave propagation will be affected not only by one position but by the whole state along its propagation path.

Thus, in order to find the appropriate thermal-hydraulic index for comparison in the case of a propagating wave, we define a very new index, which we call the cumulated TFI (CTFI). The CTFI is the value that indicates the amount of TFI that the acoustic wave experiences along the central axis. We define the CTFI at the position x p in the propagation direction by

I c (x p , z alt , R) = xp xs R 0 2π 0 I(r, θ, z alt )w(r)dθdrdx (11)
where I c is the CTFI value, I(x, r, θ) is the TFI value at x, r, θ (a graphical representation of these parameters is on Figure 10 Using a larger w(r) allows us to take into account the whole ultrasonic beam. A more complex function will be needed to take into account the beam divergence.

In this work, we only use the first definition of the CTFI, i.e. the CTFI in 1D, as a first analysis. The other possible choices may be examined in future work.

The left image on Figure 9 represents the CTFI values versus altitude z.

Each line represents the x position, i.e. propagation distance (the acoustic source is placed at x = -0.128 m). The farther the x position becomes, the larger the magnitude of the CTFI becomes as well. We find that the peak of the CTFI positions is around 0.13 m to 0.15 m. For the x position = 0.00 m just after the first (left side of) the high TFI zone and the middle of the central jet, the CTFI peak is slightly shifted to a higher z altitude. This is caused by the slight difference of the shape and position of the high TFI zones, as one can see in Figure 5.

The right image of Figure 9 shows the second derivatives of the CTFI curves.

The second derivatives has been calculated by calculating twice the central finite differences i.e.: f

(x i ) = f (xi+1)-f (xi-1)

2h

, where h = x i+1 -x i = x i -x i-1 . The forward/backward finite difference was used at the each end of the x. The inflection points are found around z = 0.08 m and around z = 0.18 m for the position x = 0.105 (The yellow line on Figure 9). These inflection points are the same altitudes as what we defined as the merging point and combining point of the standard TFI. Thus, from this result, we find that it is possible to define the merging and combining points of the TFI as the inflection points of the second derivatives of the CTFI curve. We propose thus a numerical criteria to found this specific thermal-hydraulic points.

In the following part of this work, we will study the acoustic fluctuation state based on these merging and combining points. Our objective is to propose several acoustical measurement to find the best sensitive to the behaviour of the cumulative data and so the best sensitive to find these two particular points. showing the mean and standard deviation values. We confirm that the standard deviation is the largest at altitude z = 0.13 m (number 3 of Figure 11). This is in good agreement with the peak z position of the CTFI in Figure 9. At lower altitude z = 0.04 m and 0.09 m, the distribution of the impact points exhibits directivity, i.e. the standard deviation for the y-axis direction is larger than for the z-axis direction. This result illustrates the fact that the 3D temperature fluctuation pattern before maturing of the mixing state has directivity. It is a very interesting result that indicates we can no more support classical modelling of the fluctuation of the acoustic celerity field using an isotropic Gaussian random process for the regions where the flow is still strong. From this result, we can confirm the observation by [START_REF] Iooss | Numerical simulation of transit-time ultrasonic flowmeters: uncertainties due to flow profile and fluid turbulence[END_REF] of the non-applicability of an isotropic Gaussian random field for this region in a quantitative way.

Standard deviation analysis

The left image of Figure 12 Figure 14 shows the standard deviations on the scatter plots and those moving average curves (left), and the second derivatives (right) of the times of flight. As the other derivative curves presented previously, i.e. for maximum amplitude and r, these curves also exhibit inflection points at an altitude between z = 0.11 m and z = 0.12 m for the first point and z = 0.19 m for the second one. So, for this data set, the time of flight seems not to be the better acoustic indicator to find the characteristic points of the thermal-hydraulic fluid but these two measurements are still in rather good agreement with the previous results. As the time of flight (TOF) is directly linked to variation of the sound velocity and is often the best ultrasonic indicator, we have to be cautious with this last result. We calculated the TOF at the impact point which corresponds to the maximum amplitude in one plane. Another solution would have been to calculate the minimum TOF in this plane. Of course it will need further studies as, for a practical point of view, we should consider the resulting TOF after the integration of all 'sounds paths' on the surface of the transducer.

Conclusions

Developing a very strong numerical scheme, with huge computation power, we studied the effect of a realistic heterogeneous temperature field on wave propagation in four dimensions (i.e., three spatial dimensions and time). We first simulated wave propagation at a fixed instant of a Plajest's CFD simulation, changing the altitude of the volume scanned by transducer, to analyse variations caused by the flows of liquid metal on wave propagation. To quantify the effect of the 3D temperature fied fluctuations on an acoustic path we defined a new measurement index, called CTFI (Cumulated Temperature Fluctuation Intensity). We exhibited specific altitude ranges, corresponding to the merging of jets, at which a strong effect may occur on acoustic wave propagation. We thus defined two altitudes at which the state of temperature fluctuation changes: one is the merging point, where the two high TFI zones between each of three jets start to merge, and the other is the combining point of the TFI, where the merging of the TFI zones is completed. We demonstrate that these altitude can be detected as inflection points of the second derivative of the CTFI curves.

We then carried out the same acoustic simulation but for multiple CFD We also found that the second derivatives of the standard deviation of acoustic indicators have inflection points at altitudes close to those at which we defined the merging and combining points. The results for the maximum amplitude and the amount of deviation r are better than for time of flight. The r parameter will require an array of transducers to be measured with high resolution. At this moment the recommended parameter would the the maximum value. Of course a combination of the two parameters will increase the confidence in the measure.

In the future it will be possible to increase the validation of the use of acoustic measurements to monitor complex thermal-hydraulic fields by developing spectral analysis to characterize the temporal fluctuation of the fluid. It would require to refine the time-step between two acoustic simulations.

Acknowledgements

This study was granted access to the French HPC resources of TGCC/CINES under allocations 2016107626, 2016047165 and 2016047748 made by GENCI. We

Figure 1 :

 1 Figure 1: (a) Geometry of the PLAJEST CFD simulation. (b) Snapshot of the CFD result at time 200.0 s, x-y cross-sectional plane at y = 0.09 m. Visualization for this image is done with

called

  Figure 1 b. Their simulation results are in good agreement with the experiment results obtained by JAEA in terms of normalized time-averaged temperature, normalized time-averaged temperature fluctuations, spectral power density and standard deviation of temperature values.

Figure 2 :

 2 Figure 2: Explanation of data processing for mesh generation and preparation of the heterogeneous medium to use for our acoustic calculations.

Figure 3 :

 3 Figure 3: One of the heterogeneous temperature fields used for our 3D wave propagation calculations. This field is taken from time step number 10, and the central altitude of the calculation domain is 0.1 m from the sodium outlet. One mesh contains 3,250,000 spectral elements and the total number of GLL grid points is 215,320,764.

Figure 4

 4 Figure 4 is the 2D cross-section at y = 0.09 m of the calculated 3D TFI field from the CFD results of PLAJEST. In this figure, a set of 1D TFI curves are superimposed to the 2D cross-section TFI image, which are calculated at several altitudes. There are three jets in the configuration of PLAJEST (the center position of the jets are x = -0.070 m, 0.0 m and 0.070 m). Between each jet, two zones with high TFI value arise by the interaction of these jet flows. From

Figure 4 : 4 . Ultrasonic simulation setups 4 . 1 .

 4441 Figure 4: TFI field on a 2D cross-section at y = 0.09 m and 1D profiles. Each 1D profile is factored on vertical direction so that 0.1 on z altitude axis becomes the maximum TFI value for visualizing purpose. Points a, b, and c are the reference points used for 1D sequential analysis

Figure 5 :

 5 Figure 5: Relation of the positions in an insonified volume (central altitude z = 0.1 m) and the geometry of PLAJEST. The points A, B and C are the positions of the reference points used for temporal analysis.

Figure 6 :

 6 Figure 6: Visualization of the 3D wave front using contouring and clipping, at position x = A) 0.105 m and B) -0.035 m calculated at altitude z = 0.14 m, time = 200.010 s of PLAJEST with a heterogeneous medium temperature.

5. 2 .

 2 Treatment of the massive amount of calculations and management of the resultsIn this study, all our simulations were performed on two of the largest supercomputers in Europe: CURIE (CEA TGCC) and OCCIGEN (CINES), both part of GENCI (Grand Équipement National de Calcul Intensif). The computation domain was divided into 256 parts, and parallelized calculations were carried out. The average duration for an acoustic simulation is about 26 minutes, excluding mesh generation and the interpolation processes of the temperature fields, which are done once and for all. The duration of the interpolation of a temperature field from a TrioCFD result to SPECFEM3D is approximately 20 hours using a single CPU for one acoustic simulation altitude of one time step.We carried out 70 time steps of interpolation of the 3D temperature field, for 22 difference altitudes, resulting in a total of 1540 acoustic simulations to perform.As a result, the total time needed for a complete simulation of wave propagation over 7 s of the variable thermal-hydraulic regime is approximately equal to 687 hours, i.e. 29 days. Because of the huge number of calculations, the result data cannot be conserved as a 3D volume data because of the limitation of allocated storage on supercomputers. Instead of storing all results in 3D, we first selected the 2D planes in which we will analyze the acoustic signals. We then defined time windows for each y-z plane to cut the received signals depending on the arrival time of the wave front. These time-windowed signals were then gathered as a single HDF5 binary file. The HDF5 (Hierarchical Data Format) format is a standard and widely-used binary file format that has been developed to manage extremely large and complex data collections. By using it, one can access the results faster than with other standard file formats such as e.g. ASCII, json, csv, pickle (Python-friendly binary data format) etc.

Figure 7 :

 7 Figure 7: Visualized 3D TFI field.

Figure 8 :

 8 Figure 8: TFI values on the x axis and at y = 0.09 m at several altitudes.

  ), x s is the x coordinate of acoustic source, z alt is the altitude (along the z axis) of the center of the acoustic source plane, R is the distance from the central axis, and w(r) is a weight function to make TFI values near the central axis have more effect and TFIs far from the central axis less effect. One can define several versions of the CTFI, for instance: 1. CTFI in 1D (integrated on the central x axis), with w(r) = 1 and R = one mesh element size, 2. CTFI in 3D A (integrated in a domain where the acoustic beam passes), with w(r) = 1 and R = the radius of the acoustic source, i.e. 1.27 cm, 3. CTFI in 3D B (integrated in a cylindrical volume where the acoustic beam passes), with a weighting function w(r) = 0.54 + 0.46cosπ r R , hamming window which is the same apodization function used for simulating the acoustic plane source in this simulation, and R = the radius of the acoustic source, i.e. 1.27 cm in our simulation.

Figure 9 :

 9 Figure 9: CTFI values on the z axis at y = 0.09 m and second derivatives. Each line indicates the distance from the acoustic source.

6 .

 6 figure[START_REF] Ploix | Investigation of ultrasonic backward energy from various edges as a function of their 2d/3d geometry and of the incidence angle, for application to ultrasonic thermometry at the outlet of a tube[END_REF]. We call this point as the impact point. We define the following three acoustic indicators:

  425

Figure 11

 11 Figure11shows the impact points, i.e. the position where the pressure value becomes maximum in the y-z receiving plane. In order to show temporal changes of the positions of the impact points more clearly, the positions of the

Figure 10 :

 10 Figure 10: Definition of the impact point, of the amount of deviation r, and of the deviation angle θ.

  impact points are linked by arrows in temporal order, and digits are added to indicate the order in which the position of the impact point changes for the seven first ones. Histograms for the y and z axis directions are also placed,

  shows the standard deviation of maximum amplitudes of simulated acoustic signals at the center of the y-z planes for each z altitude as the scatter plots and the smoothed lines calculated after a moving average triangular window with a window length of 8. The standard deviation for the max amplitude and the other parameters were calculated from all results of 70 time steps of the acoustic propagation simulations. The right image shows the second derivatives of the average curves calculated by the same way with CTFI of Figure 9. The curves of standard deviations on the left have similar peak positions as the CTFI curves of Figure 9 and also exhibit the same position of their inflection points (right) when computing their second derivatives. Thus, the inflection points of standard deviation curves of maximum amplitudes occur at the altitudes of the merging point (z = 0.08 m) and combining point (z = 0.18 m) of the TFI values. It is clearly a very interesting result as it proves that it is possible to find representative points of the thermal-hydraulic field using acoustic measurements.

Figure 13

 13 Figure 13 is the analysis of the deviation length r. The left image shows the scatter plots which are the standard deviation values of r, and those smoothed curves calculated with the same method with the maximum amplitude. The right image shows the second derivative of these smoothed curves. The peak positions are also similar with the maximum amplitude curves, i.e. the peaks of the x = 0.035 m, 0.08 m and 0.105 m curves occur around altitude z = 0.13 m to 0.14 m. Considering the inflection points of standard deviation curves of deviation length r we obtain the same results as they occur at the altitudes of (z = 0.08 m) for the first one and (z = 0.18 m) for the second one. This acoustic indicator is also well adapted to find the altitudes of merging point and combining point of the thermal-hydraulic fluid. Finally, let us analyze the simulation results based on the fluctuation of times of flight, i.e. the receiving time of the signal. We define the receiving time (the

Figure 11 :

 11 Figure 11: Movement of acoustic impact points in the y-z receiving plane at x = 0.105 m, and histograms for the y and z directions. Movements of only the seven initial time steps are indicated with red arrows and with digits in magenta.

Figure 12 :

 12 Figure 12: Left: Standard deviations of maximum amplitude at the center of each y-z receiving plane (scatter) and those smoothed curves with a moving average window. Right: Second derivatives of smoothed standard deviation curves.

Figure 13 :

 13 Figure 13: Left: Standard deviation of deviation length r (scatter) and smoothed curves with a moving average window. Right: Second derivatives of the smoothed standard deviation lines.

Figure 14 :

 14 Figure 14: The standard deviation (left) and second derivative (right) of time of flight.

Figure 15 summarizesFigure 15 :

 1515 Figure 15 summarizes the results of the analysis and compares the standard deviation curves and second derivative curves of CTFI and the three types of acoustic indicators only at x = 0.105 m. This figure shows that the peaks of

  time steps in order to investigate the relation between acoustic measurements and the thermal-hydraulic state. We calculated the standard deviations of several acoustic indicators (maximum amplitude, amount of deviation r, and time of flight) to study the fluctuation of the acoustic signals and to correlate with the thermal field fluctuations. As a first result of our standard deviation analysis, we found that at the low altitude where the mixing of flows is not enough matured, the temperature heterogeneity may cause directive fluctuation of acoustic propagation. This means that the application of an isotropic Gaussian random field to simulate temperature field is not an efficient solution as the fluctuating values generated by an isotropic Gaussian random field always have a normal distribution and so produces no directivity.

thank also Antoine Gershenfeld and Pierre Emmanuel Angéli, from CEA / DES / IRESNE / DM2S / STMF, for sharing their thermal-hydraulic data of Plajest and their advises to analyse these data. This work was developed within the framework of the MISTRAL joint research laboratory between Aix-Marseille University, CNRS, Centrale Marseille and CEA.