M Hadi

S Pailhès

R Debord

A Benamrouche

Emmanuel Drouard

T Gehin

C Botella

J.-L Leclercq

P Noe

F Fillot

V M Giordano

Transient thermal conductivity in PECVD SiN x at high temperature: the thermal signature of an on-going irreversible modification

Keywords: thermal conductivity, silicon nitride, PECVD, thin film, thermoreflectance, hydrogen desorption

teaching and research institutions in France or abroad, or from public or private research centers.

Transient thermal conductivity in PECVD SiN x at high temperature: the thermal signature of an on-going irreversible modification

Introduction

Amorphous silicon nitride (a-SiN x , with 0 < x < 1.33) films have attracted en ever growing interest in recent years thanks to their appealing physical, chemical, structural, optical and electrical properties, which make them suitable for applications in a very large panoply of technologies. Largely used into the di-5 electric stacks of integrated circuits, as passivation or capping layers [START_REF] Benoit | Hydrogen desorption and diffusion in pecvd silicon nitride. application to passivation of cmos active pixel sensors[END_REF], and as gate insulators in thin film transistors [START_REF] Duffy | Preparation, properties and applications of chemically vapor deposited silicon nitride films[END_REF][START_REF] Sinha | Reactive plasma deposited si-n films for mos-lsi passivation[END_REF][START_REF] Hasegawa | Relationship between the stress and bonding properties of amorphous sin x: H films[END_REF], this material has become today an important alternative to silicon for integrated nonlinear optical applications, such as waveguides [START_REF] Wang | Nonlinear silicon nitride waveguides based on a pecvd deposition platform[END_REF], and is at the focus of research for photonic architectures in energy harvesting applications [START_REF] Braña | Enhancing efficiency of c-si solar cell by coating nano structured silicon rich silicon nitride films[END_REF]. The possibility of tuning its properties, 10 and specifically its band gap, by tuning the nitrogen content makes it a very attractive material with manageable vibrational, electrical and optical properties [START_REF] Ibaraki | Properties of amorphous semiconducting a-s i: H a-si n x: H multilayer films and of a-si n x: H alloys[END_REF][START_REF] Lowe | The electronic properties of plasmadeposited films of hydrogenated amorphous sin x (0 ≤ x ≤ 1.2)[END_REF][START_REF] Hasegawa | Si-h and n-h vibrational properties in glow-discharge amorphous sin x : H films (0 ≤ x ≤ 1.55 )[END_REF][START_REF] Hasegawa | Bonding configuration and defects in amorphous sin x : H films[END_REF].

Still, the very same variability of these properties represents a technical chal-lenge, as it strongly depends on the growth technique used and the resulting film morphology, texture, and chemical composition, which are not always perfectly mastered. Many are the techniques used for growing amorphous SiN x thin films, which lead to quite different materials, with a large variability in properties such as density, elastic moduli, intrinsic stress and thermal conductivity. Among such techniques, the one of choice for the preparation of SiN x for electronic and optic devices is the Plasma Enhanced Chemical Vapor Deposition (PECVD), which has the advantage of growing the material at low temperature. Here as well, however, the choice of the growth parameters, such as the deposition temperature, the radio-frequency, and the plasma pressure, can dramatically affect the properties of the deposited film [START_REF] Stoffel | Lpcvd against pecvd for micromechanical applications[END_REF][START_REF] Jehanathan | Thermal stability of plasma enhanced chemical vapor deposited silicon nitride thin films[END_REF]. Such large variability comes not only from the morphology or the exact Si:N ratio, but as well from an often uncontrolled content of hydrogen, which is bounded to both Si and N. Indeed, it has been reported that hydrogen is almost always present and its content and relative amount of the two types of bonds significantly affect the structural, electrical and mechanical properties of the thin film [START_REF] Bae | Characteristics of low-temperature silicon nitride (sinx: H) using electron cyclotron resonance plasma[END_REF][START_REF] Soh | Local bonding environment of plasma deposited nitrogen-rich silicon nitride thin films[END_REF][START_REF] Braun | Hydrogen effects on the thermal conductivity of delocalized vibrational modes in amorphous silicon nitride (a-SiNx:h)[END_REF]. A panoramic of the different physical properties of SiN x , depending on the growth technique, can be found in two recent reviews from Kaloyeros et al. [START_REF] Kaloyeros | Silicon nitride and silicon nitride-rich thin film technologies: Trends in deposition techniques and related applications[END_REF][START_REF] Kaloyeros | Review-silicon nitride and silicon nitride-rich thin film technologies: State-of-the-art processing technologies, properties, and applications[END_REF]. Here it is possible to appreciate the advantage of the PECVD method over the others, concerning thermal applications of SiN x . Indeed, PECVD a-SiN x has been reported to present a very low thermal conductivity, about 0.7 W/m K at room temperature [START_REF] Lee | Heat transport in thin dielectric films[END_REF], which is of relevance, as such films contribute to the thermal balance of the electronic devices, and thus play a role in solving the thermal management problem in microelectronics. This value has been reported for thicknesses larger than 100 nm, and represents the smallest thermal conductivity among all possible growing techniques for such thicknesses, the values being almost the double for the Atmospheric Pressure Chemical Vapor Deposition (APCVD) and the Chemical Vapor Deposition at 1173 K [START_REF] Kaloyeros | Review-silicon nitride and silicon nitride-rich thin film technologies: State-of-the-art processing technologies, properties, and applications[END_REF][START_REF] Lee | Heat transport in thin dielectric films[END_REF], and going up to 13 W/mK in Low Pressure Chemical Deposited SiN x [START_REF] Zhang | Thermal conductivity and diffusivity of freestanding silicon nitride thin films[END_REF][START_REF] Sikora | Highly sensitive thermal conductivity measurements of suspended membranes (sin and diamond) using a 3ω-völklein method[END_REF][START_REF] Sikora | Erratum:"highly sensitive thermal conductivity measurements of suspended membranes (sin and diamond) using a 3ω-völklein method[END_REF][START_REF] Jain | Measurement of the thermal conductivity and heat capacity of freestanding shape memory thin films using the 3ω method[END_REF]. Such a difference has been partially explained by a larger presence of voids, or micro-pores, in the microstructure of PECVD SiN x with respect to the other growing techniques.

Still, it is expected that hydrogen content and exact stoichiometry play a role as well.

Several reports exist in the literature on the room temperature thermal conductivity of PECVD SiN x , with values largely scattered between ∼ 0.45 and 4.5 W/mK, a consequence of the different deposition parameters used in the different works, often not mentioned apart from the temperature [START_REF] Kim | Thermal property evaluation of a silicon nitride thin-film using the dual-wavelength pump-probe technique[END_REF][START_REF] Stojanovic | Thin-film thermal conductivity measurement using microelectrothermal test structures and finite-element-model-based data analysis[END_REF][START_REF] Eriksson | Thermal characterization of surface-micromachined silicon nitride membranes for thermal infrared detectors[END_REF].

Recently, Braun et al. have pointed to the hydrogen content as the critical parameter behind such variability. They have reported an exhaustive study of the effect of hydrogen content on the room temperature thermal conductivity of PECVD SiN x deposited at 673 K. By annealing the films at different temperatures above 923 K, in order to gradually release hydrogen, the authors show that thermal conductivity is essentially determined by the hydrogen content, which is responsible for the increase of the average inter-atomic bond lengths and thus a decrease of the density. More specifically, after releasing all hydrogen at the highest annealing temperature, the thermal conductivity increases by a factor as large as 2.5 independently on the actual Si-N composition. Concerning the temperature dependence of the thermal conductivity, the most complete study to our knowledge is the one of Lee et al. [START_REF] Lee | Heat transport in thin dielectric films[END_REF], who investigated it in PECVD SiN x for a deposition temperature of 573 K, gas pressure of 1.8 Torr and RF power of 16W. They have reported a thermal conductivity increasing from ∼ 0.45 to ∼ 0.7 W/mK for thicknesses between 20 and 100 nm, and then stabilizing for thicker samples. This behavior was ascribed to an effect of the interface thermal resistance, whose importance, stronger for thinner films, would become negligible for films thicker than 100 nm. The reported values, indeed, were obtained neglecting all interface contributions. As such, the increase with the thickness would only be apparent, an artefact of the analysis. Similarly, the observed increase of the thermal conductivity with temperature between 78 and 400 K could be overestimated, as it included the interface thermal resistance, which is expected to decrease with temperature.

While the work from Braun et al. has helped in understanding the large variability of results in the literature, and has pointed to the hydrogen content as the key parameter determining the room temperature value for the thermal conductivity, it cannot fully account for the large scattering of reported values over an order of magnitude. It appears then that other systematic studies still need to be performed in order to get an exhaustive understanding of the thermal properties of PECVD SiN x and their dependence on the growth conditions.

The technological interest of such study is evident for engineering optimized stackings in microelectronics and for controlling the thermal balance in many applications. Indeed, the zoology of different results for different, not necessarily known, growth conditions hinders any understanding of the dependence of this quantity on the growth parameters. Moreover, the investigation of its temperature dependence, as a function of the growing parameters, is critical for applications where high temperatures are typically needed.

In this work, we present an exhaustive investigation of the temperature dependence of the thermal conductivity of PECVD SiN x between 300 and 773 K, as a function of the deposition temperature, chosen between 298 and 573 K, keeping all other growth parameters fixed. As we will show, a linear dependence of the room temperature thermal conductivity with the deposition temperature is found, which can be ascribed to an increase of the elastic moduli at a high deposition temperature. Surprisingly, for all deposition temperatures, the thermal conductivity exhibits a steep decrease with temperature above 473 K, featuring a minimum around 673 K, before increasing back to values comparable to those at room temperature. As we will show, this behavior, observed only the first time the sample is heated, is associated to an irreversible modification of the thin film, and may be related to a partial desorption of the hydrogen trapped in a-SiN x during the deposition. The article is organized as follows: in section 2 we give experimental details on sample preparation, measurement techniques and analysis method, in section 3 we present our experimental results, and in section 4 a discussion on the microscopic origin of the observed behavior is provided. Conclusions and perspectives are given in section 5. A first set of samples was deposited at 573 K, with three SiN x thicknesses: 200, 292 and 501 nm. This last thickness was also deposited on SiO 2 substrate, in order to establish the independence of the thermal transport measurements on the kind of substrate. A second set was deposited at substrate temperatures of 298 K, 323 K, 373 K, 473 K, 523 K, and 573 K and thickness of 300 nm. Layers of ∼ 14 and ∼ 30 nm were also deposited in the same conditions for density measurements by X ray reflectivity (XRR). Samples were then characterized by FTIR spectroscopy to identify the presence of oxygen and hydrogen, as well as the evolution of chemical bonds with deposition temperature. Experimental details and spectra for XRR and FTIR characterization are reported in the Supplementary Material.

Experimental Methods

Amorphous

Thermal conductivity was measured using the thermoreflectance technique [START_REF] Yuan | Nanosecond transient thermoreflectance method for characterizing anisotropic thermal conductivity[END_REF][START_REF] Cahill | Analysis of heat flow in layered structures for time-domain thermoreflectance[END_REF] with the commercial apparatus TF-LFA from Linseis. The pump laser had a pulse width of 8 ns, a wavelength of 1040 nm and an energy of 2.6 mJ, while the probe laser was continuous, and had a wavelength of 523 nm and 3.8 µW power. The signal was collected in time by a Si-PIN-Photodiode with an aquisition rate of 2.5 GHz. Samples were coated with a 5 nm thick Cr layer, as a diffusion barrier, and a 100 nm thick Au layer, which was impulsively heated by the pump laser. The time relaxation of its temperature was deduced from the time evolution of its temperature-dependent reflectivity, as measured using the probe laser. High temperature measurements were performed under Argon atmosphere, in a furnace chamber with optical aperture, able to go up to 773 K. The thermoreflectance signal was collected at temperature steps of 50 K, after the time needed for temperature stabilization. For each sample and temperature, we collected 8 thermoreflectance signals, changing the position on the sample, in order to get rid of eventual inhomogeneities and increase our statistics.

The signal has been fitted using the Cahill model adapted to time domain measurements [START_REF] Cahill | Analysis of heat flow in layered structures for time-domain thermoreflectance[END_REF], which takes into account the temperature-dependent contribution to temperature relaxation from all the layers in the stacking, here Au, Cr, SiN x and Si. The only unknown parameter in the fitting procedure is the thermal conductivity of SiN x . It is worth remembering that the total thermal resistance of the stacking corresponds to a series of thermal resistances:

R tot = R Au + R Au-Cr + R Cr-SiNx + R SiNx + R SiNx-Si + R Si (1)
Where R A-B is the thermal resistance at the interface between A and B. In our program, we fit an effective thermal conductivity [START_REF] Lambropoulos | Thermal conductivity of dielectric thin films[END_REF] for SiN x which includes all the interface thermal contributions:

k meas = d SiNx R SiNx,ef f (2)
with

R SiNx,ef f = R SiNx + R IF = d SiNx k SiNx + R IF (3)
and d the SiN x layer thickness. By fitting our measured values of

d SiNx
kmeas as a function of the SiN x thickness, we can get from the slope the real SiN x thermal conductivity and from the intercept the interfacial thermal resistance. This latter resulted to be between 2 and 6 ×10 -8 m 2 K/W, in good agreement with what reported by Lee et al. [START_REF] Lee | Heat transport in thin dielectric films[END_REF]. Within our error bar, we can consider its temperature dependence negligible. More details on the fitting model and procedure, and the slope method for the extraction of the thermal conductivity are reported in the Supplementary Material.

Results

Sample Characterization

In order to study the dependence of the thermal transport properties on deposition temperature, we have prepared a set of samples deposited on silicon at different temperatures T dep between 298 and 573 K, for which, from literature, it is expected that both film microstructure and local chemical structure will be different.

Indeed, at low deposition temperatures, the incoming free radicals will have less mobility and will be essentially stuck where they arrive on the surface, with no optimization of their position, leading to a fractal-like network with many embedded micro-voids [START_REF] Liao | Oxidation of silicon nitride prepared by plasma-enhanced chemical vapor deposition at low temperature[END_REF][START_REF] Walmsley | Effects of deposition temperature on the mechanical and physical properties of silicon nitride thin films[END_REF]. Increasing deposition temperature, the larger mobility allows a more packed arrangement, leading to higher densities, and indeed linearly increasing densities going from 1.8 to 2.9 g/cm 3 were reported for deposition temperatures between 373 and 573 K [START_REF] Walmsley | Effects of deposition temperature on the mechanical and physical properties of silicon nitride thin films[END_REF].

In Fig. S1 we report our experimental density, measured by XRR as a function of deposition temperature. Only a very weak increase is observed, by 4%, suggesting that in our growth conditions the microporosity is very weakly affected by the deposition temperature.

Our FTIR characterization revealed the presence of hydrogen bounded to both Si and N for deposition temperatures larger than 323 K, and oxygen incorporated in the film at all deposition temperatures, also bound to Si and N.

Increasing deposition temperature leads to a reduction of the oxygen content and an increase of the Si-N bonds, in agreement with previous reports [START_REF] Soh | Chemical structure of low-temperature plasma-deposited silicon nitride thin films[END_REF]. Instead, we could not quantify any change in the hydrogen content (See Supplementary Material).

Thermal conductivity study

We first measured the thermal conductivity of the as-grown films upon annealing up to 773 K ("first annealing"). After cooling to room temperature, the measure was repeated following the same thermal protocol ("second annealing"). In Fig. 1 we report some representative thermoreflectance spectra after subtraction of the electronic noise, together with their fit. The spectra have been collected upon heating on a 292 nm thick SiN x sample deposited on silicon at a deposition temperature of 573 K. In Fig. 3(a) we report the temperature dependence of the effective thermal conductivity for the same sample, as representative for all measured thicknesses, the overall behavior being similar among all the films. As shown in the Supplementary Material, the thermal conductivity extracted using the slope method only slightly differs from the effective one for a thickness of 292 nm, reported here, which can thus be used as a reference for our discussion. First of all, our room temperature values, k 300K = 0.63 ± 0.03 W.m -1 K -1 (on the as-grown) and k 300K = 0.70 ± 0.03 W.m -1 K -1 (after the first annealing ramp), are consistent with the ∼ 0.7 W.m -1 K -1 reported by Lee [START_REF] Lee | Heat transport in thin dielectric films[END_REF] on a 252 nm PECVD SiN x film deposited at the same temperature, with ρ = 2.2 g/cm 3 . On the other hand, Lee reports an increasing thermal conductivity with temperature, which slightly flattens above 400 K, matching the usual behaviour observed in amorphous materials [START_REF] Wingert | Thermal transport in amorphous materials: a review[END_REF]. In our sample, we can at least say that above 300 K, we don't observe an increase in thermal conductivity, but a surprising decrease with temperature. On the first heating ramp, κ(T ) decreases linearly upon heating from 300 K to 450 K and then drops abruptly down to a minimum of ∼ 0.45 W/mK at 673 K, before increasing again. This behavior is very specific to the first heating ramp indicating that it results from an irreversible modification of our sample. Upon a second annealing, κ(T ) exhibits a linear decrease up to ∼ 450 K following almost the same slope as for the first heating ramp and then saturates at higher temperature. This latter behavior of κ(T ) is then reproducible if the sample is annealed again.

The presence of this minimum is quite intriguing. A first possibility is that it is an artefact due to a non-correct account for thermal interface resistance.

To rule out this hypothesis, we have deposited the same sample on a SiO Once assessed that the peculiar temperature dependence is real and an intrinsic property of our PECVD SiN x , it is interesting to see how it depends on the deposition temperature. We first report in Fig. 3 the measured room temperature thermal conductivity (κ 300K ) as a function of the deposition temperature (T dep ). κ 300K increases with deposition temperature with a linear rate of change of 0.27 ±0.02 W/m.K 2 . Despite this, the global temperature dependence during the first annealing is found to be the same in all samples, independently on the deposition temperature, as can be observed in Fig. 4 (a), where we report κ(T ) in all samples during the first annealing (data have been vertically shifted for optical convenience). The same curves are reported in Fig. 4 (b) after normalization to the respective room temperature value.

It may be seen that in all samples κ(T ) undergoes a first linear decrease in the 300 K -450 K temperature range followed by a drop, a minimum and an increase at higher temperatures.The slope of the first linear decrease is found to increase with the deposition temperature by a factor as large as ∼ 50% at the highest deposition temperature, going from -0.

Discussion

First of all, we shall comment on the dependence of the room temperature thermal conductivity on the deposition temperature. The observed increase can be due to a smaller porosity of SiN x when deposited at higher temperatures, as reported in the literature [START_REF] Liao | Oxidation of silicon nitride prepared by plasma-enhanced chemical vapor deposition at low temperature[END_REF][START_REF] Walmsley | Effects of deposition temperature on the mechanical and physical properties of silicon nitride thin films[END_REF], so that we can model our material as effectively composed of SiN x and micropores [START_REF] Liao | Oxidation of silicon nitride prepared by plasma-enhanced chemical vapor deposition at low temperature[END_REF]. The effective medium theory allows to estimate the expected thermal conductivity of a porous medium from the zeros of the function where k 1 , k 2 , k e and x are the thermal conductivity of non porous SiN x , pores and the effective thermal conductivity of the porous medium respectively and

(1 -x) k 1 -k e k 1 + 2k e + x k 2 -k e k 2 + 2k e = 0 (
x is the pores volume fraction [START_REF] Gong | A novel effective medium theory for modelling the thermal conductivity of porous materials[END_REF]. As we don't know the thermal conductivity of an equivalent non-porous SiN x , we use k 1 = 0.7 W/mK for estimating how much porosity should change for giving the observed thermal conductivities.

For pores, we take k 2 = 0. Using the measured thermal conductivity values at room temperature, we obtain pores volume fractions ranging from 3% up to 24% from the highest to the lowest deposition temperature. Such a large change in porosity is however not compatible with the measured weak density change, which is at most 4%. This could be explained by the limitations of the model, which does not take into account the additional thermal resistance introduced by the SiN x -pore interfaces. Still such a strong difference in the porosity estimated from density and from thermal conductivity points to another origin for the large change in this latter.

As we have seen, the chemical structure evolves with the deposition temper-ature: oxygen content decreases and Si-N bonds increase. SiN x being a dielectric, the main heat carriers are the atomic vibrations, whose quasi-particle is the phonon. Not only the interfaces can scatter phonons, but also the atomic defects in the Si-N network, such as oxygen and hydrogen bonds. Decreasing their content means decreasing scattering from defects, leading thus to larger thermal conductivities. Moreover, the Si-N bond content increasing with deposition temperature, the elastic constants are expected to increase as well, due to the higher rigidity of this bond, and indeed an increase of the Young modulus by a factor of 2 has been reported in literature for deposition temperatures between 373 and 573 K [START_REF] Walmsley | Effects of deposition temperature on the mechanical and physical properties of silicon nitride thin films[END_REF]. From the kinetic theory for a phonon gas, the thermal conductivity reads k T = 1 3 C v vl, with C v the specific heat, v the sound velocity and l the phonon mean free path. Considering that the longitudinal and transverse velocities are given by v L/T = M L/T /ρ with M L/T the longitudinal and shear elastic modulus and ρ the mass density, we expect the thermal conductivity to increase with the square root of the elastic moduli. A rough estimation using the Young moduli reported in the literature gives a change of the thermal conductivity by a factor of ∼50% between 373 and 573 K, which well corresponds to the observed increase.

Second, we shall comment on the unusual decrease with temperature of the thermal conductivity, independent on the deposition conditions. As mentioned, this is opposite to what is usually observed in amorphous materials, where thermal conductivity at high temperature increases due to the increase of the specific heat, while the phonon mean free path is almost constant at its minimum value.

It is in fact a very weak decrease, as the slope of the thermal conductivity is in the 0.4-0.6 10 -3 W/mK 2 range for an absolute value of about 0.7 W/mK. A similar behavior was previously reported on CVD amorphous Si 3 N 4 , and explained as the result of a relevant interfacial thermal resistance, three orders of magnitude larger than in our case, and of its temperature dependence. The intrinsic thermal conductivity was shown to weakly increase with temperature [START_REF] Griffin | The effective transverse thermal conductivity of amorphous si 3 n 4 thin films[END_REF]. In our case, the interfacial thermal resistance is quite small and its temperature dependence negligible within our experimental uncertainty. Indeed, as mentioned and as shown in the Supplementary Material, the intrinsic thermal conductivity is not very different from the effective one as measured on the 292 nm thick film. Still, the slope of the initial decrease is reduced in the intrinsic thermal conductivity, leading to an almost constant value with temperature prior to the decrease leading to the minimum.

We come finally to the presence of a minimum in the thermal conductivity during the first annealing, at ∼ 673 K, which disappears in the second annealing ramp in all samples, independent on deposition temperature and substrate. A possible origin of this behavior is that our amorphous SiN x starts crystallizing at high temperature. However, the crystallization temperature in PECVD SiN x has been reported to be much larger, ∼ 1148 K [START_REF] Jehanathan | Crystallization and compositional changes in amorphous PECVD sin x thin films[END_REF]. Moreover, even if some incipient crystallization were present, this should lead to an increase of the thermal conductivity, being the one of the crystalline phase about 26 W/mK [START_REF] Kaloyeros | Review-silicon nitride and silicon nitride-rich thin film technologies: State-of-the-art processing technologies, properties, and applications[END_REF].

A second possibility is the activation of a temperature induced structural relaxation. The higher temperature induces higher atomic mobility, which can lead to slight atomic rearrangements, responsible for the evolution of elastic and thermal properties. This is expected to have a role for temperatures larger than the deposition temperature, thus we would expect the steep decrease of the minimum to appear at lower temperatures for lower deposition temperatures. This is however not the case. Moreover, in case of structural relaxation, we would expect a monotonic and asymptotic behavior, from the unrelaxed to the more relaxed state, while here we face a steep decrease followed by a similarly steep increase back to an unperturbed value.

In fact, the observed phenomenology resembles to a temporary thermally induced perturbation of our system: thermal transport is progressively hindered by some temperature-activated mechanism and then comes back to the unperturbed behavior when the mechanism is over. Being this phenomenon independent on deposition temperature, we cannot expect it to be related to the presence of oxygen nor the amount of Si-N bonds, which depend both on this latter. On the other hand, we also know that there is hydrogen in our sample, both bound to Si and N, and under a non-bonded molecular form, trapped within the porous structure [START_REF] Jafari | Composition limited hydrogen effusion rate of a-sinx: H passivation stack[END_REF], and we could not see important differences among samples deposited at different temperatures. In order to see whether hydrogen could be responsible for the minimum in the thermal conductivity, we have measured the FTIR spectrum of our SiN x deposited at 573 K before and after the annealing, but we could find indication of only a very slight decrease of the Si-H and N-H bonds content (see Supplementary Material). Still, FTIR is sensitive only to bounded hydrogen, so that we cannot rule out an effect of molecular H 2 . Detecting small quantities of molecular hydrogen in thin films is not trivial. The typical H-content in PECVD SiN is about 1.7 10 22 atoms/cm 3 [START_REF] Xie | Stability of hydrogen in silicon nitride films deposited by low-pressure and plasma enhanced chemical vapor deposition techniques[END_REF], of which only a few percent will be under the form of molecular hydrogen [START_REF] Jafari | Composition limited hydrogen effusion rate of a-sinx: H passivation stack[END_REF], i.e. about 10 20 atoms/cm 3 . This amount has remained elusive to our tentative Raman measurements. More advanced techniques are necessary to detect it in thin films, such as effusivity measurements [START_REF] Jafari | Composition limited hydrogen effusion rate of a-sinx: H passivation stack[END_REF], time of flight elastic recoil detection analysis [START_REF] Siketić | Quantitative analysis of hydrogen in thin films using time-of-flight elastic recoil detection analysis[END_REF] or a self-consistent analysis of the light element NRA data combined with an RBS analysis for heavy elements [START_REF] Lanford | Nuclear reaction analysis for h, li, be, b, c, n, o and f with an rbs check[END_REF].

An exhaustive study of the hydrogen effusivity as a function of temperature in PECVD SiN x deposited at 673 K with different stochiometries has been recently reported by Jafari et al. [START_REF] Jafari | Composition limited hydrogen effusion rate of a-sinx: H passivation stack[END_REF]. It is found that, starting from 673 K, the thermal energy is enough to break first the N-H bond and then, above ∼973 K, the more energetic Si-H bond, so that all hydrogen diffuses out of the sample between 673 and 1273 K. Interestingly, for N-rich compositions, a small peak appears in the hydrogen effusivity, with a maximum between 650 K and 670 K depending on the N content, accounting for about the 2-5% of the total hydrogen effusivity signal. This peak has been ascribed to the premature release of hydrogen from surface desorption of H 2 molecules. The correspondence of the temperature range, as well as the similar peak shape points to the diffusion and release of molecular hydrogen as the cause of the observed minimum in the thermal conductivity. The SiN x matrix with freely moving H 2 can be assimilated to the case of a superionic material, where some atoms are characterized by a liquid-like mobility. In these latter, the active ions movement may enhance phonon scattering, due to disorder and anharmonicity, as reported in the high temperature super-ionic phase of CuCrSe 2 [START_REF] Niedziela | Selective breakdown of phonon quasiparticles across superionic transition in (cucrse) 2[END_REF]. We may expect a similar phenomenology here. In our case the large H 2 mobility is limited to the temperature -or time at a given temperature -range within which all the molecular hydrogen diffuses out of the sample, thus the perturbation on thermal transport disappears once all the molecular hydrogen has left, being an irreversible process. The peak shape can be understood as due to a scattering probability proportional to the amount of mobile H 2 . The fact that we find a similar shape, and more specifically depth, of the minimum throughout all samples may indicate that the quantity of mobile H 2 is quite insensitive on the deposition temperature, or that there is a saturation effect.

Conclusions

We have here reported an exhaustive investigation of the thermal conductivity of PECVD SiN x thin films as a function of the deposition temperature, all other growing parameters being kept constants. A strong variation of the room temperature value, of a factor of almost 2, is found as a function of the deposition temperature, likely related to the increase of the Si-N bonds content, and thus a larger rigidity of the amorphous network at higher deposition temperature. Our data can explain the panoply of results on thermal conductivity reported in the literature, measured on samples grown in different conditions, and can help in aiming the best deposition temperature for the foreseen application, on the basis of the needed thermal properties. Moreover, we find that, independently on the deposition temperature, the effective thermal conductivity weakly decreases with temperature, apparently at odds with other reports [START_REF] Lee | Heat transport in thin dielectric films[END_REF]. This decrease is in fact reduced even more once the interfacial thermal resistance is properly taken into account, leading to an almost temperature independent thermal conductivity below 500 K. Indeed, we have observed the onset of an irreversible process at ∼500 K, leading to a minimum of the thermal conductivity at about 673 K, before it increases back to the unperturbed value. This peculiar behavior does not depend on the deposition temperature, indicating that poros-ity, Si-N bonds and oxygen bonds are not responsible for it. Instead, we can explain it as due to a premature release of molecular hydrogen, whose diffusivity is temperature-activated at about the same temperatures [START_REF] Jafari | Composition limited hydrogen effusion rate of a-sinx: H passivation stack[END_REF]. The shape of the minimum and its irreversibility can be understood in terms of an enhanced phonon scattering from mobile H 2 , which is proportional to the content of this latter. Once all the mobile H 2 molecules have left the sample, the process is over and the thermal conductivity comes back to the unperturbed value. Our results assess the stability of the temperature window in which the phenomenon takes place, and the important effect on the thermal conductivity, which decreases by 30%, together with its relatively long time constant (∼20 minutes ).

As reported in Braun et al. [START_REF] Braun | Hydrogen effects on the thermal conductivity of delocalized vibrational modes in amorphous silicon nitride (a-SiNx:h)[END_REF], for higher temperatures, such as Si-H and N-H bonds break down and all hydrogen diffuses out of the sample, an irreversible absolute modification of the thermal conductivity takes place once these point defects have been annihilated. We could thus expect, increasing temperature, a transient behavior even more pronounced than in the present study, as the amount of mobile hydrogen is larger. Our findings indicate that in order to use SiN x in electronic or optic applications, when application temperature is larger than 500 K, a preliminary annealing up to 723 K is fundamental for assuring its thermal stability. For application temperatures larger than 723 K, much higher annealing temperatures will be needed to get rid of all the bound hydrogen. This means as well that, in high temperature applications, it is not possible to tailor the thermal conductivity by tuning the hydrogen content, as suggested in Ref. [START_REF] Braun | Hydrogen effects on the thermal conductivity of delocalized vibrational modes in amorphous silicon nitride (a-SiNx:h)[END_REF], because this would introduce transient mechanisms leading to a time evolution of thermal properties. Finally, we stress here the high sensitivity demonstrated by PECVD SiN x thermal conductivity with respect to mobile hydrogen, which can open the way to its use to monitor hydrogen diffusion, and possibly to the application of PECVD SiN x as membranes for hydrogen stocking and release.

  SiN x materials have been deposited by Plasma Enhanced Chemical Vapor Deposition (PECVD) on Silicon substrate using a 200 mm OXFORD Plasmalab 80+ PECVD system. SiH 4 4% in He, NH 3 , Ar and N 2 were used as reactive gases plasma mixture. Plasma gases were adjusted with SiH 4 /He = 80 sccm, Ar = 200 sccm, NH 3 = 10 sccm and N 2 = 600 sccm. 13.56MHz RF power was held at 70 W under a pressure of 1.2 Torrs. Temperature of deposition was varied from 298 K to 573 K. Pre-conditioning of the deposition chamber was operated during 5 minutes before each introduction of sample.

Figure 1 :

 1 Figure 1: Thermoreflectance signal measured upon heating on a 292 nm thick SiNx sample deposited on silicon at a deposition temperature of 573 K. The signal, in blue, is reported after subtraction of the electronic noise and normalization to 1. The red curve is the fit result, performed using the Cahill method as described in the Supplementary Material. The temperature of the measurement is 298 (a), 373 (b), 473 (c), 573 (d), 673 (e) and 773 K (f).

Figure 2 :

 2 Figure 2: Temperature dependence of the thermal conductivity in a SiNx thin film. (a): Thermal conductivity (κ) of a 292 nm thick SiNx sample deposited on Si, reported as a function of temperature (T ) during the very first heating ramp on the as-grown sample (blue circles) and during a second heating ramp run after the very first one (red squares). (b) Thermal conductivity on 501 nm thick SiNx deposited on Si (blue circles) and SiO 2 (black filled squares), reported as a function of temperature during the very first heating ramp. All error bars correspond to the standard deviation of the distribution of values obtained from 8 different measurements at each temperature.

2 substrate,Figure 3 :

 23 Figure 3: Room temperature thermal conductivity (κ 300K ) of SiNx on Si grown by PECVD as a function of the deposition temperature (T dep ). The error bars correspond to the standard deviation of the distribution of values obtained from 8 different measurements at each temperature. The solid black line is a linear fit.

Figure 4 :

 4 Figure 4: Temperature dependence of the thermal conductivity of SiNx thin films for different deposition temperatures. (a) Thermal conductivity (κ(T ) as a function of temperature (T ) during the very first heating ramp of samples deposited on Si at different temperatures. The deposition temperature is indicated next to each curve. Dashed lines follow the linear decrease in the first temperature range, solid lines mark the shape of the minimum (b) The same curves in (a) are here reported after normalization to the room temperature value. The color/symbol legend is the same as in (a). The error bars correspond to the standard deviation of the distribution of values obtained from 8 different measurements at each temperature.

Figure 5 :

 5 Figure 5: Time evolution of the thermal conductivity of a SiNx thin film at high temperature. The measured thermal conductivity of SiNx on Si grown by PECVD at 523K heated to 673K as a function of time. The red line is the fit of the data to an exponential function. The error bars correspond to the standard deviation of the distribution of values obtained from 8 different measurements at each temperature.

Aknowledgements

The authors are grateful to D. Cahill for the fruitful discussions. Funding: this work was supported by the Idex Lyon Scientific Breakthrough program [IPPON project] and the Agence National de la Recherche [MAPS-ANR-20-CE05-0046].


Cédric Leblond

Vincent Mélot

Jean-François Sigrist

Christian Lainé

Bruno Auvity

Hassan Peerhossaini

Transient fluid forces on a rigid circular cylinder subjected to small-amplitude motions

Keywords: circular cylinder, shock loading, transient dynamics, analytical method, fluid forces, history effect, confinement effect

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the case of shock loading that results from underwater explosion, naval propulsion devices can be subjected to highly accelerated, high frequency motions. These components can be in contact with a fluid, which is the case for heat exchangers, nuclear propulsion reactors, or squeeze-film dampers. In order to improve design margins and ensure safety, long life, and satisfactory operating performance of the shock loaded components, the precise knowledge of the fluid effects on the solid components is of major importance. Since the body motion modifies the fluid flow, and the fluid flow influences the body motion, complex fluid/structure interaction problems have to be accounted for. In order to choose an appropriate fluid model, design engineers must identify the meaningful physical phenomena occurring in their particular geometry. This can be done by estimating the fluid forces on simplified representative geometries, for which analytical results are available.

The purpose of the present paper is to formulate and group such models for a rigid circular cylinder subjected to relatively small amplitude transient motions along a radial line, in an annular or infinite fluid domain. Numerous analytical and numerical methods are now available for the harmonic dynamics of this geometry ͓1,2͔. It also continues to be an active domain of interest particularly to studies of flow-induced vibrations. However, transient dynamics have comparatively received little attention. Although the solution of a time-dependent problem can be obtained from the solution of the corresponding harmonic problem using Fourier synthesis methods or series of resonance modes ͓3,4͔, timedomain methods can provide more physical insight. Such methods are currently receiving considerable attention ͓5-8͔, in particular, due to increasing computer capabilities and constantly improving computational methods.

This paper deals with several aspects of fluid effects in one of the simplest two-dimensional fluid problems, illustrated in Fig. 1. The motion of the inner cylinder is imposed and the outer cylinder, when present, is fixed. Consequently, this is not a fluidstructure interaction problem and emphasis can be especially placed on the forces acting on the moving cylinder. Highly confined flows, such as those encountered in squeeze-film dampers or rotating machinery, are not discussed here. Related fluid effects can be found in other articles ͓9-14͔. In this study, the imposed motion is assumed sufficiently small, that no boundary-layer separation occurs ͓15͔. Hence, there is neither wake formation nor its force related modifications ͓16,17͔. Two-phase flows are also not considered. Some elements concerning the associated damping effect can be found in a review article ͓18͔.

The problem of interest is governed by seven independent dimensional parameters: R 1 , R 2 , , , c, U, and . R 1 and R 2 are the geometrical parameters denoting the inner and outer cylinder radii, respectively; , , and c are the fluid parameters representing the density, the kinematic viscosity, and the speed of sound in the fluid at rest. Lastly, U and denote the characteristic velocity and pulsation frequency of the imposed inner motion. For instance, U may be the maximum velocity and the pulsation with the maximum energy in the power spectrum of the transient motion. Since there are three fundamentals units, length, mass and time, the problem is driven by four independent dimensionless numbers. While there are of course infinite ways to build such numbers The present paper treats the transient fluid forces experienced by a rigid circular cylinder moving along a radial line in a fluid initially at rest. The body is subjected to a rapid displacement of relatively small amplitude in relation to its radius. Both infinite and cylindrically confined fluid domains are considered. Furthermore, non-negligible ampli-tude motions of the inner cylinder, and viscous and compressible fluid effects are ad-dressed, successively. Different analytical methods and models are used to tackle each of these issues.

For motions of non-negligible amplitude of the inner cylinder, a potential flow is assumed and the model, formulated as a twodimensional boundary perturbation problem, is solved using a regular expansion up to second order. Subsequently, viscous and compressible effects are handled by assuming infinitesimal amplitude motions. The viscous fluid forces are formulated by solving a singular perturbation problem of the first order. Compressible fluid forces are then determined from the wave equation. A nonlinear formulation is obtained for the non-negligible amplitude motion. The viscous and com-pressible fluid forces, formulated in terms of convolution products, are linked to fluid history effects induced by wave propagation phenomena in the fluid domain. These mod-els are expressed with dimensionless parameters and illustrated for a specific motion imposed on the inner cylinder. The different analytical models permit coverage of a broad range of motions. Hence, for a given geometry and imposed displacement, the appropriate fluid model can be identified and the resulting fluid forces rapidly estimated. The limits of these formulations are also discussed.

from the seven dimensional numbers, it will be seen below that the four following ones are particularly meaningful:

␣ = R 2 R 1 , = U R 1 , ␤ = R 1 2 , ⍀ = R 1 c ͑1͒
2 Motion of Non-Negligible Amplitude in a Potential Flow 2.1 General Equations. Here, a rapidly moving circular cylinder in a two-dimensional annular incompressible and inviscid fluid region is considered. The fluid forces on the inner cylinder for an infinitesimal displacement are found in Ref. ͓19͔, in which fluid advection and geometrical deformation are neglected. An extension of this model to larger amplitude motions is proposed by formulating the problem as a boundary-perturbation problem ͓20͔ solved up to the second order. Only the main steps of the derivation are described here; a more detailed account can be found in Ref. ͓21͔. Another way of accounting for largedisplacement effects can be found in Refs. ͓13,14͔. In this section, the fluid forces are obtained by solving the dimensional equations and are then expressed in a dimensionless form at the end.

Since the flow is assumed inviscid and initially irrotational, it remains irrotational at later times ͓22͔ and can be described by the Laplace equation for the velocity potential ⌽ as follows:

ٌ 2 ⌽ =0 ͑2͒
The normal component of the fluid velocity at the boundaries must equal the velocity of the cylinders, so that the boundary conditions can be written as ٌ⌽ • e r = 0 on the fixed outer cylinder ͑3͒ ٌ⌽ • n = e ˙͑t͒e x • n on the moving inner cylinder ͑4͒

where e ˙͑t͒ is the inner-cylinder velocity, e x the unit vector along the x-axis, and n the unit outward normal, as illustrated in Fig. 1.

In order to make this problem analytically tractable, the position of the inner cylinder and the unit normal in Eq. ͑4͒ must be explicitly expressed.

2.2 Boundary-Perturbation Formulation. By considering the triangle OOЈA in Fig. 1, the position of the inner circular cylinder C͑t͒ can be characterized by the polar equation

r c ͑,t͒ = e͑t͒cos + R 1 ͱ1-e 2 ͑t͒ R 1 2 sin 2 ͑5͒
where e͑t͒ is the inner-cylinder displacement. This formula is only valid for ͉e͑t͒ / R 1 ͉ Ͻ 1. Expanding the square root in terms of series gives the relation

r c ͑,t͒ = R 1ͭ 1 + cos e͑t͒ R 1 + ͚ n=1 ϱ ͑-1͒ n ͑sin ͒ 2n 1 n! ͟ k=0 n-1 ͩ 1 2 -kͪ ϫ ͫ e͑t͒ R 1 ͬ 2n ͮ ͑6͒
By introducing the perturbation parameter ͑t͒,

͑t͒ = e͑t͒ R 1 ͑7͒
approximate positions of the inner cylinder can be evaluated at a given order in powers of ͑t͒. Hence, by neglecting third-order terms and higher in Eq. ͑6͒, the second-order approximation for the inner-cylinder position takes the form r c ͑,t͒ϷR 1 ͓1+͑t͒cos -1 2 2 ͑t͒sin 2 ͔͑ 8͒ It is also useful to write in explicit terms the unit outward normal n on the moving inner cylinder so as to express the boundary condition ͑4͒. Thus, the parametric curve ⌿ of C͑t͒ is introduced as follows:

C͑t͒: ۋ ⌿͑͒ = O + r c ͑,t͒e r ͑,t͒͑ 9͒
where O denotes the center of the outer cylinder. The unit tangent T to C͑t͒ at the position is given by: Here, ␣ is the dimensionless geometrical parameter called the confinement number. , often called the Keulegan-Carpenter number, quantifies the displacement of the inner cylinder in relation to its radius and takes into account fluid advection and geometrical deformation. ␤, the Stokes number, is the ratio between unsteady and viscous effects, and ⍀ is used to take into account the fluid compressibility in unsteady flows.

Section 2 of this paper considers motions of non-negligible amplitude in an incompressible and inviscid flow. The problem is formulated as a boundary-perturbation method for the velocity potential and a regular expansion is used to solve it up to the second order. Hence, fluid advection and geometrical deformations due to the inner motion are approximately taken into account. The resulting fluid force expression depends nonlinearly on the imposed instantaneous displacement, velocity, and acceleration, and can be seen as an extension of the Fritz model ͓19͔, which is limited to very small amplitude motions. Its dimensionless formulation highlights the Keulegan-Carpenter number as the meaningful parameter in quantifying large-displacement effects. In Sec. 3, the fluid viscosity is considered assuming a small amplitude imposed motion that is rapid enough that the boundary layer thickness can be taken as much less than the radial clearance. The corresponding fluid forces are obtained using a singular perturbation method solved at the first order. Its expression includes a convolution product whose amplitude is quantified by the Stokes number. Section 4 takes fluid compressibility into account, assuming a small amplitude motion and inviscid flow. The forces are derived by solving the wave equation for the velocity potential and take the form of a convolution product whose kernel depends on the compressibility number. The viscous and compressible fluid force expressions extend the harmonic models of Chen ͓2͔ for transient motions. In Sec. 5, the fluid forces obtained using these models are illustrated and compared for a specific imposed inner-cylinder motion.

T͑͒ = ⌿Ј͑͒ ʈ⌿Ј͑͒ʈ ͑10͒ with ⌿Ј͑͒ = r c Ј͑͒e r + r c ͑͒e ͑11͒
The prime denotes a derivative according to . The unit normal n orthogonal to T can then be evaluated. Its truncation at the second order in powers of ͑t͒ is

n͑,t͒Ϸ R 1 ʈ⌿Ј͑͒ʈ ͫ cos + ͑t͒cos 2 -3 8 2 ͑t͒͑cos -cos 3͒ ͬ e x + R 1 ʈ⌿Ј͑͒ʈ ͫ sin + ͑t͒sin 2 -3 8 2 ͑t͒͑sin -sin 3͒ ͬ e y

͑12͒

The boundary condition on the moving inner cylinder ͑4͒ can now be seen as the extreme boundary condition of the following family:

ٌ⌽͓r n ͑,t͒,͔ • n͓r n ͑,t͒,͔ = e ˙͑t͒e x • n͓r n ͑,t͔͒͑ 13͒
where r n ͑ , t͒ is the nth-order approximation of the power series expansion of the inner-cylinder position r c ͑ , t͒. Performing a Taylor expansion of Eq. ͑13͒ about R 1 permits turning the original problem into an equivalent sequence of problems with solutions ⌽ 0 , ⌽ 1 , ⌽ 2 , ⌽ 3 , . . . that can be found recursively. Each ⌽ n must satisfy the following Laplace equation and outer boundary condition:

ץ 2 ⌽ n ץr 2 + 1 r ץ⌽ n ץr + 1 r 2 ץ 2 ⌽ n ץ 2 =0 ͑14͒ ͯ ץ⌽ n ץr ͯ r=R 2 =0 ͑15͒
Due to Eqs. ͑8͒, ͑12͒, and ͑13͒ and several Taylor expansions, the boundary conditions on the inner cylinder up to the second order are given by

ͯ ץ⌽ 0 ץr ͯ r=R 1 = e ˙͑t͒cos ͑16͒ ͯ ץ⌽ 1 ץr ͯ r=R 1 = 2␣ 2 ␣ 2 -1 e ˙͑t͒cos 2 ͑17͒ ͯ ץ⌽ 2 ץr ͯ r=R 1 = 2␣ 2 ͑␣ 2 -1͒ 2 ͑␣ 2 +1͒ e ˙͑t͒cos + 3␣ 2 ͑␣ 4 +1͒ ͑␣ 2 -1͒ 2 ͑␣ 2 +1͒ e ˙͑t͒cos 3 ͑18͒
where Eqs. ͑17͒ and ͑18͒ are anticipated knowledge of the solutions ⌽ 0 and ⌽ 1 , respectively. Hence, each ⌽ n is the solution of the Laplace equation on a fixed annular fluid region with Neumann boundary conditions. The solution of the general problem is then achieved using

⌽ = ͚ n=0 ϱ n ͑t͒⌽ n ͑19͒ be found in Ref. ͓21͔.
Once the velocity potential is determined, the pressure in the fluid domain can be obtained from the Bernoulli relation ͓22͔

p =- ץ⌽ ץt - 1 2 ٌ͑⌽͒ 2 + C ͑20͒
where is the fluid density and C a constant available in the whole fluid domain. Since the flow is assumed inviscid, the integrated fluid forces per unit length on the moving inner circular cylinder can be written as

F͑t͒ =-͵ 0 2 p͑r c ͑,t͒͒I ញ • n͑,t͒ʈ⌿Ј͑͒ʈd ͑21͒
with I ញ the identity matrix. By performing Taylor expansions about r = R 1 in the pressure term of the above equation, the fluid forces up to the second order can be found. Only the e x -direction component of F is nonzero and its dimensional expression is given by

͑t͒ =-R 1 2 ␣ 2 +1 ␣ 2 -1 e ¨͑t͒ + 2␣ 2 ͑␣ 2 -1͒ 3 ͑␣ 2 +1͒ e͑t͒e ˙2͑t͒ - 4␣ 4 ͑␣ 2 -1͒ 3 ͑␣ 2 +1͒
e 2 ͑t͒e ¨͑t͒͑ 22͒

The first term on the right-hand side, which is the main-order term, is exactly the Fritz model ͓19͔. It is proportional to the inner-cylinder acceleration and consists only in an added mass effect. The next two terms are displacement-dependent nonlinear corrections. The first one, proportional to the square of the innercylinder velocity, is induced by the fluid advection, and the second one can be seen as an added mass induced by relatively large amplitude motions. For an infinite fluid domain ͑␣ → ϱ͒, the force expression reduces to -R 1 2 e ¨͑t͒, the well-known fluid displaced term. It is fruitful to express the fluid forces with dimensionless quantities. By normalizing t, e͑t͒, e ˙͑t͒, e ¨͑t͒, and ͑t͒ according to Table 1, the resulting second-order nondimensional fluid forces can be written as

͑t͒ =- ␣ 2 +1 ␣ 2 -1 e ¨͑t͒ + 2 2␣ 2 ͑␣ 2 -1͒ 3 ͑␣ 2 +1͒ ͓e͑t͒e ˙2͑t͒ -2␣ 2 e 2 ͑t͒e ¨͑t͔͒

͑23͒

The Keulegan-Carpenter number appears explicitly and allows quantification of the large-displacement nonlinear effects.

Viscous Model

3.1 Problem Formulation. Small amplitude motions of the inner-circular cylinder in a viscous incompressible fluid are now addressed. For an infinite fluid domain, the harmonic problem was solved by Stokes ͓23͔. An approximate solution is proposed here for transient motion in a confined flow that is achieved using a matched asymptotic expansion at first order, also known as the composite boundary-layer solution ͓15͔. The incompressible

U Displacement U / Acceleration U Potential R 1 U Pressure R 1 U Force R 1 2 U
If the perturbation parameter is sufficiently small, the series will converge rapidly and a few terms will suffice for a good approximation of the solution.

Fluid Forces Evaluation.

The results are obtained up to the second order, so the solution is expected to be valid in cases where ͓e͑t͒/R 1 ͔ 3 Ӷ1. Here,⌽ 0 ,⌽ 1 , and⌽ 2 are found by solving Eqs. ͑ 14͒-͑18͒ by separation of variables. Their expressions can Navier-Stokes equations are the starting point. All variables are normalized as indicated in Table 1. The nondimensional governing equations take the form

ץu ץt + ͑u • ٌ͒u =-ٌp + ␤ -1 ٌ 2 u ͑24͒ ٌ • u =0 ͑25͒
where u denotes the fluid velocity and p the pressure. The Keulegan-Carpenter number and the Stokes number ␤ appear explicitly. Since a small amplitude motion in relation to the innercylinder radius is assumed here, the Keulegan-Carpenter number is small and the advection term can be neglected. Moreover, the no-slip boundary conditions can then be applied on fixed boundaries as follows:

͉u͉ r=1 = e ˙͑t͒e x , ͉u͉ r=␣ =0 ͑26͒
The inner-cylinder motions are also assumed to be rapid enough that the viscous boundary layer will be much smaller than the radial clearance. Hence, the fluid domain can be divided into two parts: one far from the boundaries, where viscous effects can be neglected, and the other closer to the boundaries, where the viscous dissipation is a main-order term ͓24͔.

Outer Resolution.

Far from the boundaries, the following outer expansion in ␤ -1/2 is introduced into the nondimensional Navier-Stokes equations:

ͩ u͑r,,t͒ p͑r,,t͒ ͪ=ͩ u 0 ͑r,,t͒ p 0 ͑r,,t͒ ͪ+␤ -1/2 ͩ u 1 ͑r,,t͒ p 1 ͑r,,t͒ ͪ+O͑␤ -1 ͒͑ 27͒

Keeping only the main-order terms gives

ץ t u 0 =-ٌp 0 ͑28͒ ٌ • u 0 =0 ͑29͒
Since Eq. ͑28͒ has lost the second-order spatial derivative, Eq. ͑26͒ cannot be fully satisfied and only the normal projection of the boundary conditions is kept. The fluid flow in this outer region can therefore be characterized by a velocity potential ⌽ 0 governed by the following system:

ץ 2 ⌽ 0 ץr 2 + 1 r ץ⌽ 0 ץr + 1 r 2 ץ 2 ⌽ 0 ץ 2 =0 ͑30͒ ͯ ץ⌽ 0 ץr ͯ r=1 = e ˙͑t͒cos , ͯ ץ⌽ 0 ץr ͯ r=␣ =0 ͑31͒
The solution of this system is readily obtainable by separation of variables; it is given by

⌽ 0 ͑r,,t͒ =- 1 ␣ 2 -1 ͩr+ ␣ 2 r ͪe ˙͑t͒cos ͑32͒
The nondimensional main-order velocity components far from the boundaries can then be deduced as follows:

u 0 r = 1 ␣ 2 -1 ͩ ␣ 2 r 2 -1ͪe ˙͑t͒cos ͑33͒ u 0 = 1 ␣ 2 -1 ͩ ␣ 2 r 2 +1ͪe ˙͑t͒sin ͑34͒ r u ¯͑,,t͒ u ¯͑,,t͒ p ¯͑,,t͒ = 0 u ¯0 ͑,,t͒ p ¯0͑,,t͒ + ␤ -1/2 u ¯1 ͑,,t͒ u ¯1 ͑,,t͒ p ¯1͑,,t͒ + O͑␤ -1 ͒ ͑35͒
where u ¯, u ¯, and p ¯are, respectively, the radial velocity, the azimuthal velocity, and the pressure in the boundary-layer. The introduction of the boundary-layer variable together with the inner expansion, Eq. ͑35͒, account for the different orders of magnitude of the radial and azimuthal boundary layer velocities: the radial velocity is normalized with the viscous diffusion velocity given by ͱ and the azimuthal velocity with the characteristic cylinder velocity. Introducing Eq. ͑35͒ into the Navier-Stokes equations and retaining the main-order terms yield the following system of equations:

ץp ¯0 ץ =0 ͑36͒ ץu ¯0 ץt =- ץp ¯0 ץ + ץ 2 u ¯0 ץ 2 ͑37͒ ץu ¯1 ץ + ץu ¯0 ץ =0 ͑38͒
As expected, the unsteady and diffusion terms in Eq. ͑37͒ are now of the same order. Equation ͑36͒ states that the pressure at the main order is constant across the boundary layer. Therefore, the pressure term in Eq. ͑37͒ can be linked to the outer variables by considering the azimuthal component of Eq. ͑28͒ as follows:

-

ץp ¯0 ץ = lim r→1 ͩ-1 r ץp 0 ץ ͪ=ͯ ץu 0 ץt ͯ r=1 = ␣ 2 +1 ␣ 2 -1 e ¨͑t͒sin ͑39͒
Only the azimuthal boundary-layer velocity u ¯0 is required for evaluation of the first-order fluid forces, since the inner radial velocity contributes to the higher-order terms. Thus, u ¯0 must satisfy the no-slip boundary condition on the moving cylinder and must match the outer azimuthal velocity. It is fully characterized by the system

ץu ¯0 ץt = ץ 2 u ¯0 ץ 2 + ␣ 2 +1 ␣ 2 -1 e ¨͑t͒sin ͑40͒ ͉u ¯0 ͉ =0 =-e ˙͑t͒sin ͑41͒ lim →ϱ u ¯0 = ␣ 2 +1 ␣ 2 -1 e ˙͑t͒sin ͑42͒
Taking the Laplace transform of the above system, the inner azimuthal solution can be found in the Laplace domain,

u ˆ0 ͑,,s͒ = ͩ-2␣ 2 ␣ 2 -1 e -ͱ s + ␣ 2 +1 ␣ 2 -1 ͪ e ¨͑s͒ s sin ͑43͒
where s is the Laplace variable and the circonflex denotes the Laplace transform of the corresponding temporal function.

Evaluation of Fluid Forces.

The uniformly valid firstorder dimensionless azimuthal velocity in the fluid domain is obtained by summing the inner and outer solutions ͑34͒ and ͑43͒ and removing the common limit ͑42͒, which is counted twice ͓24͔,

u ˆ͑r,,s͒ = ͫ 1 ␣ 2 -1 ͩ1+ ␣ 2 r 2 ͪ-2␣ 2 ␣ 2 -1 e -ͱ ␤s͑r-1͒ ͬ e ¨͑s͒ s sin

͑44͒

The resultant dimensionless fluid forces on the inner cylinder can be evaluated in the Laplace domain by the formula ͓22͔

where u 0 and u 0 ⌽ denote, respectively, the radial and azimuthal velocities.

Inner Resolution.

Due to the no-slip boundary condition, a viscous boundary-layer is formed close to the moving cylinder.

In this fluid part, the boundary layer variable =␤ 1/2 ͑r-1͒ is injected into the Navier-Stokes equations. Moreover, the following inner expansion is performed for the nondimensional dependent variables

F ˆ͑s͒ =- 1 ͵ 0 2 ͉p ˆ͉r=1 e r d + ␤ -1 1 ͵ 0 2 ٌ͉͑u ˆ+ ٌ͑u ˆ͒T ͉͒ r=1 • e r d

͑45͒

The first term on the right-hand side can be expressed in terms of the fluid velocities, Eqs. ͑33͒ and ͑44͒, using integration by parts. Keeping only the main-and first-order terms in the power expansion of ␤ -1/2 gives the integrated fluid forces in the Laplace do- main. The nondimensional fluid forces in the temporal domain can then be found by using the convolution theorem and some inverse Laplace transform properties. Only the e x -component of F is nonzero and is given by

͑t͒ =- ␣ 2 +1 ␣ -1 e ¨͑t͒ -␤ -1/2 ␣ 2 ␣ 2 -1 4 ͱ ͵ 0 t e ¨͑͒ ͱ t - d ͑46͒
The first term is the small amplitude potential added mass effect found in the previous section. The second one, inversely proportional to the square root of the Stokes number, is the first-order viscous rectification. The total history of the imposed motion is involved so that its effect can be attributed neither to added mass nor viscous damping alone. The convolution product kernel, known explicitly, indicates that the imposed motion at time has decreasing influence on the fluid forces at time t as t -increases. This type of convolution product ͑which accounts for viscous ef-fects͒ is well known and can be found in Ref. ͓25͔ for an infinite plate. For an infinite fluid domain, the resulting dimensionless fluid forces are then given by

͑t͒ =-e ¨͑t͒ -␤ -1/2 4 ͱ ͵ 0 t e ¨͑͒ ͱ t - d ͑47͒
which consists of the displaced fluid term and a first-order viscous correction. This expression can also be found, with higher-order viscous corrections, from existing harmonic models ͓26,27͔ based on a Fourier transform method ͓28͔. Nevertheless, for relatively high Stokes number flows, the first-order term suffices to account for the viscous effects ͓29͔.

Compressible Model

4.1 Problem Formulation. Here, small amplitude motions of the inner circular cylinder in an inviscid compressible fluid are considered. The solution of the associated harmonic problem can be found in Ref. ͓2͔, for instance. The fluid dynamics is assumed to be governed by the wave equation for the velocity potential

ٌ 2 ⌽ -⍀ 2 ץ 2 ⌽ ץt 2 =0 ͑48͒
where the variables are normalized according to Table 1 and the compressibility number ⍀ appears explicitly. The boundary conditions are of Neumann type and can be written as

ͯ ץ⌽ ץr ͯ r=1 = e ˙͑t͒cos , ͯ ץ⌽ ץr ͯ r=␣ =0 ͑49͒
4.2 Evaluation of Fluid Forces. The system of Eqs. ͑48͒ and ͑49͒ can be solved by using the Laplace transform for the temporal variable and performing a separation of variables. The resulting nondimensional velocity potential in the Laplace domain is then given by ⌽ ˆ͑r,,s͒ = e ¨͑s͒ ks

Y 1 Ј͑␣k͒J 1 ͑kr͒ -J 1 Ј͑␣k͒Y 1 ͑kr͒ Y 1 Ј͑␣k͒J 1 Ј͑k͒ -Y 1 Ј͑k͒J 1 Ј͑␣k͒ cos ͑50͒
convolution theorem yield the nondimensional e x -component of the fluid forces in the temporal domain as follows:

͑t͒ =-͵ 0 t C͑t -͒e ¨͑͒d ͑51͒
where the kernel C͑t͒ must be evaluated. This requires a Laplace inversion according to a contour integration based on the Bromwich integral as follows:

C͑t͒ = 1 2i ͵ -iϱ +iϱ C ˆ͑s͒e st ds ͑52͒ with C ˆ͑s͒ =- 1 k Y 1 Ј͑␣k͒J 1 ͑k͒ -J 1 Ј͑␣k͒Y 1 ͑k͒ Y 1 Ј͑␣k͒J 1 Ј͑k͒ -Y 1 Ј͑k͒J 1 Ј͑␣k͒ ͑53͒
This inversion can be obtained numerically. Figure 2 shows the temporal evolution of the kernel C͑t͒ for ⍀ = 1 and a few different values of ␣. These functions have points of finite discontinuity that account for the reflections of hydrodynamic waves in the fluid domain. The intervals between the discontinuities are equal to the time for a wave to travel from the inner to the outer cylinder and back. For instance, for ⍀ = 1 and ␣ = 2, the pressure wave induced by the impulse motion of the inner cylinder at t = 0 is reflected back on the outer cylinder at t = 1 and influences the integrated fluid forces on the inner cylinder at t = 2. Furthermore, the more confined the fluid domain is, the more discontinuities the kernel contains. At the beginning of the process, that is, before the wave returns to the inner cylinder, the kernel is seen to be independent of the confinement ratio ␣. In evaluating the fluid forces at a time t in a confined domain, the motion of the inner cylinder at an early time is as important as that at time t, even if t -is large. On the other hand, for an infinite fluid domain, the motion at large tdoes not significantly influence the fluid forces at t, since the waves are not reflected back.

5 Fluid Forces on a Rigid Circular Cylinder Subjected to a Sine-Wave Shock 5.1 Sine-Wave Shock. The fluid forces obtained in Secs. 2-4 are illustrated and compared for a specific motion imposed on the inner cylinder that corresponds to a unique sinusoidal period of acceleration and is then stopped. Hence, the dimensional motion is fully determined by where k=i⍀s, the circonflex denotes the Laplace transform. Here, J 1 and Y 1 are the Bessel functions of order 1, respectively, of the first and second kinds. The prime denotes the derivative according to the function argument. Making use of Eq. ͑50͒ and applying the

e ¨͑t͒ = Ά U 2 sin t, ∀t ͓0,2/͔ 0, ∀t Ͼ 2/ · ͑54͒
where is the pulsation of the shock. In what follows, it is also convenient to express this motion in a nondimensional form. Using the normalization factor of Table 1 yields the simple representation

e ¨͑t͒ = Ά 1 2 sin t, ∀t ͓0,2͔ 0, ∀t Ͼ 2 · ͑55͒
fluid domain, the greater the discrepancies with the Fritz model, becoming as much as 20% for ␣ ഛ 2. The different limits of the curves in Fig. 4 can be explained by recalling that the nonlinear model was derived assuming an inner-cylinder displacement smaller than its radius and smaller than the radial clearance.

Viscous Effects.

Here, the viscous fluid force model given by Eq. ͑46͒ is applied to the specific motion of Eq. ͑55͒. The temporal evolution normalized by the maximum potential forces is displayed in Fig. 5 Even if this function is partly composed of a harmonic function, it can be considered as a purely transient motion since it starts from rest at t=0 and stops at t=2. It consists of a simple imposed translation and is representative of the early residual motion undergone by equipment during shock loading.

Non-Negligible Amplitude Motion Effects.

Here, the dimensionless fluid forces given by Eq. ͑23͒ are illustrated on the specific imposed motion defined by Eq. ͑55͒. Figure 3 shows the temporal evolution of these forces with confinement number ␣ =2 for different Keulegan-Carpenter numbers . They are normalized with the maximum value given by the potential Fritz model ͓19͔, which assumes negligible amplitude motions. It can be seen that the nonlinear and Fritz models predict similar forces when the inner cylinder is still close to its initial position, but show different behaviors as the inner cylinder approaches the outer one. This results in a force increase that becomes more significant as the Keulegan-Carpenter number increases. This effect is checked with a finite-volume numerical code able to take into account moving boundaries ͓30͔. The nonlinear model is in good agreement with these numerical results for inner-cylinder displacements up to 60% of the radial clearance. Some discrepancies are observed for higher displacements. More specifically, the nonlinear model underestimates the maximum fluid forces, and hence more than three terms should be kept in Eq. ͑19͒ to reproduce the strong geometrical deformation. Nevertheless, the nonlinear model predicts the trend better than the Fritz model.

Figure 4 shows the maximum fluid forces obtained with the nonlinear model as function of the Keulegan-Carpenter number for different confinement numbers ␣, and also normalized by the maximum Fritz model forces. At a given , the more confined the -R 1 ͑or, in dimensionless form, ␤ 1/2 ͑␣ -1͒ ӷ 1͒.

The maximum fluid forces normalized by the maximum potential forces as function of the Stokes number are illustrated in Fig. 6 for different confinement numbers. Large discrepancies are observed for low Stokes numbers, an effect that is less pronounced for confined domains. This is because the potential added mass increases more rapidly with the confinement than the viscous term. This result is not expected to be valid for strongly confined domains in which viscous forces become predominant.

Compressible Effects.

The compressible fluid forces of Eq. ͑51͒ are now applied to the specific motion described through Eq. ͑55͒. The temporal evolution, normalized by the maximum potential forces, is shown in Fig. 7 for ␣ = 2 and different compressible numbers ⍀. It can be seen that compressible waves are generated in the fluid domain and result in fluid force oscillation around the potential curve due to their multiple reflections on the boundaries. Once the inner cylinder has been stopped, the compressible waves continue to induce these harmonic oscillations without damping. For higher compressible numbers than those in Fig. 7, several compressible waves can be generated and clearly identified on a single curve. In an infinite fluid domain, the loss of energy induced by the outgoing waves has a damping effect ͓2͔.

The maximum fluid forces normalized by the maximum potential value are displayed in Fig. 8 as function of the compressible number ⍀ and different confinement numbers ␣. For small compressible numbers, the fluid forces tend to their potential value for all confinements. For large compressible numbers, they tend to zero due to the increasing difference between the imposed motion frequency and the fundamental fluid cavity acoustic frequency. Pseudoresonances are observed around ⍀ Ϸ 1. This effect is all the more pronounced as the domain becomes more and more confined. The maximum fluid forces are then seen to be three times greater than their corresponding potential value. Several peaks on the curve ␣ = 4 can be observed, a reminder that a compressible fluid domain cannot be viewed as a single-degree-of-freedom spring-mass system since its harmonic dynamics contains an infinite number of resonant frequencies ͓2͔.

Illustration of the Use of the Models on a Given

System. The three models were described independently in Secs. 2-4 as a function of their corresponding dimensionless numbers. They can, however, be grouped on a single figure for a particular geometry, so that their range of validity can be discussed and the significant fluid effects in a particular case can be identified. For instance, by fixing the geometrical parameters R 1 and R 2 , the fluid viscosity , the speed of sound c, and the imposed velocity U, the maximum normalized forces can be displayed as a function of the dimensional imposed pulsation . Hence, increasing causes the Stokes number and the compressibilty number increase and the Keulegan-Carpenter number decrease. Such curves are given in Fig. 9 with R 1 = 0.2 m, U = 0.5 m s -1 , c = 1500 m s -1 , =10 -6 m 2 s -1 , and several R 2 . For imposed pulsations such that ഛ 1s -1 , viscous effects must be taken into account. For 1 s -1 ഛ ഛ 10 2 s -1 , the flow is dominated by potential effects and for 10 2 s -1 ഛ ഛ 10 5 s -1 , fluid compressibility must be considered.

Conclusion

Transient fluid loading on a circular cylinder moving along a radial line is addressed. The fluid domain can be infinite or cylindrically confined. A nonlinear fluid force model that takes into account the geometrical deformation and fluid advection is de- All of these models emphasize some fluid effects in simple asymptotic cases. However, all physical phenomena are of course not treated in this paper. For instance, large-displacement motions in a potential flow and viscous effects are handled independently, and this prevents the models from tackling boundary-layer separation issues. Furthermore, the viscous fluid forces are derived by assuming laminar flow. This assumption does not hold for high-Stokes-number flows, even at very low Keulegan-Carpenter numbers, due to the presence of centrifugal hydrodynamic instabilities ͓31,32͔ in the boundary layer. In harmonic flows, these threedimensional phenomena are known to increase the fluid damping ͓33,34͔ by a factor 2. Another restriction comes from the monophasic assumption, which makes the models incapable of accounting for cavitation phenomena. The present paper gives, however, general guidance that can be used by design engineers for the purpose of estimating the fluid force level and, if required, selecting the meaningful fluid model to use in further investigation.

Fig. 1

 1 Fig. 1 Geometrical configuration

Fig. 2

 2 Fig. 2 Evolution of the compressible fluid forces kernel with ⍀ = 1, for different confinement numbers ␣: "¯… 1.5, "---… 2, and "-… 10

Fig. 3

 3 Fig. 3 Evolution of the nondimensional fluid forces normalized with the maximum of the Fritz fluid force model †19 ‡ for ␣ = 2. Comparison between the Fritz model "¯… and the nonnegligible displacement model for different Keulegan-Carpenter numbers : "---… 0.19 and "-… 0.28 Fig. 4 Maximum nondimensional fluid forces, obtained with the non-negligible displacement model and normalized with the maximum of the Fritz fluid force model †19 ‡, in function of the Keulegan-Carpenter number and for different confinement numbers ␣: "¯… 4, "-•-… 2, "---… 1.5 and "-… 1.1

Fig. 6

 6 Fig. 6 Maximum nondimensional fluid forces, obtained with the viscous model and normalized with the maximum of the Fritz fluid force model †19 ‡, in function of the Stokes number and for different confinement numbers ␣: "¯… ؕ, "-•-… 2, "---… 1.5 and "-… 1.1

Fig. 8

 8 Fig. 8 Maximum nondimensional fluid forces, obtained with the compressible model and normalized with the maximum of the Fritz fluid force model †19 ‡, in function of the compressibility number and for different confinement numbers ␣: "¯… 4, "-•-… 2, "---… 1.5 and "-… 1.1

Fig. 9

 9 Fig. 9 Maximum fluid forces obtained with the different models and normalized with the maximum of the Fritz fluid force model †19 ‡, in function of the dimensional imposed motion pulsation, for different confinement numbers ␣: "¯… 4, "---… 2, and "-… 1.1

Transient Fluid Forces on a Rigid Circular Cylinder Subjected to Small Amplitude Motions Cédric Leblond, Vincent Melot, Jean-François Sigrist, Christian Lainé

	DCNS Propulsion, Service Technique et Scientifique, 44620 La Montagne, France
	Bruno Auvity, Hassan Peerhossaini
	Laboratoire de Thermocinétique, CNRS UMR 6607, Rue Christian Pauc, BP 50609,44306 Nantes, France

Table 1 Normalization factors for the variables

 1

	Variable	Normalization factor
	Length	R 1
	Time	1 /
	Velocity

Downloaded from http://asmedigitalcollection.asme.org/pressurevesseltech/article-pdf/130/3/031302/5530324/031302_1.pdf by McGill University user on 15 November 2022