R Viennois
email: romain.viennois@umontpellier.fr

P Hermet

D Machon

M M Koza

D Bourgogne

B Fraisse

A P Petrović

D Maurin

Stability and lattice dynamics of Ruddlesden-Popper tetragonal Sr2TiO4

We report a combined experimental and theoretical lattice dynamics study of the Ruddlesden-Popper layered compound Sr2TiO4. From inelastic neutron scattering experiments we derive the generalized phonon density of states of Sr2TiO4. We also report its heat capacity, thermal expansion and the thermodynamic Grüneisen parameter using the calculated bulk modulus and find a large value of about 2. Using Raman scattering experiments under pressure, we discuss a probable structural distortion of the tetragonal structure above 11 GPa which could be due to non-hydrostatic compression. The mode Grüneisen parameters of the four Raman-active modes are determined and shown to be in reasonable agreement with those obtained by DFPT calculations. The temperature behavior of the Raman-active modes was studied, allowing us to determine the implicit volume and explicit anharmonic contributions.

Above 400 K, the implicit volume contribution dominates the temperature-induced variation of the four Raman-active modes, whereas below this temperature, the explicit anharmonic contribution is the dominant contributor to the highest energy mode. Our results underline the importance of the anharmonicity in vibrational-related properties of Sr2TiO4.

Introduction

The Ruddlesden-Popper (RP) structure forms a family of layered body-centered tetragonal crystals (space group I4/mmm) with general formula An+1MnO3n+1, where A is an alkaline or rare-earth metal and M is a transition metal. [START_REF] Ruddlesden | New compounds of the K2NiF4 type[END_REF][START_REF] Ruddlesden | The compound Sr3Ti2O7 and its structure[END_REF] These structures are derived from cubic perovskites and can be viewed as a natural superlattice made by successively stacking n AMO3 perovskite-type layers separated by one AO Rocksalt-type layer.

Ruddlesden-Popper compounds with n = 1 display a wide variety of microscopic transport properties, making them very interesting for technological applications in the energy and photonics sectors. In particular, we highlight high-temperature superconductivity [START_REF] Bednorz | Possible High Tc Superconductivity in the Ba-La-Cu-O System[END_REF]4 ; thermoelectricity, photoconductivity and photocatalysis in titanates [START_REF] Lee | Ruddlesden-Popper phases as thermoelectric oxides: Nb-doped SrO(SrTiO3)n (n = 1, 2)[END_REF][START_REF] Lee | Preparation and thermoelectric properties of heavily Nb-doped SrO(SrTiO3)1 epitaxial films[END_REF][START_REF] Wang | Thermoelectric properties of electron doped SrO(SrTiO3)n (n = 1, 2) ceramics[END_REF][START_REF] Paranthaman | Photoelectrochemical studies on the mixed oxides, SrTiO3, Sr2TiO4 and Sr3Ti2O7[END_REF][START_REF] Yi | Photochromism and visible light induced H2 generation in Sr2TiO4: Cr complexes[END_REF][START_REF] Castelli | Stability and bandgaps of layered perovskites for one-and two-photon water splitting[END_REF][START_REF] Jia | Composite Sr2TiO4/SrTiO3(La,Cr) heterojunction based catalyst for hydrogen production under visible light radiation[END_REF][START_REF] Sun | Photocatalytic Hydrogen Production over Chromium Doped Layered Perovskite Sr2TiO4[END_REF] ; protonic and/or oxygen conduction in titanates [START_REF] Shimura | High temperature protonic conduction in Sr2TiO4-based ceramics with K2NiF4-type structure[END_REF][START_REF] Sirikanda | effect of Co doping on oxygen permeation in Sr3Ti2O7 with Ruddlesden-Popper structure[END_REF] and nickelates [START_REF] Zhao | Electrochemical evauation of La2NiO4+ as a cathode material for intermediate temperature solid oxide fuel cells[END_REF][START_REF] Ananyev | Oxygen isotope exchange in La2NiO4[END_REF] ; and a range of optical properties, especially photoluminescence in titanates. [START_REF] Cockroft | Site-selective spectroscopy of defect chemistry in SrTiO3, Sr2TiO4 and Sr3Ti2O7[END_REF][START_REF] Inaguma | Photoluminescence of Praesodymium-Doped Srn+1TinO3n+1[END_REF][START_REF] Mu | Luminescence properties of Eu 2+ and Ho 3+ in Sr2TiO4[END_REF] During the last decade, various complex correlated properties have been discovered in Sr2MO4 compounds with d-electrons around the Fermi level. [START_REF] Matsuno | Variation of the Electronic Structure in Systematically Synthesized Sr2MO4 (M = Ti, V, Cr, Mn, and Co)[END_REF][START_REF] Weng | Electronic structure and optical properties of layered perovskite Sr2MO4 (M = Ti, V, Cr and Mn): An Ab initio study[END_REF][START_REF] Nakamura | eg-level splitting in a layered perovskite manganite as revealed by charge modulation spectroscopy[END_REF][START_REF] Sakurai | Synthesis Conditions and Magnetic Properties of Sr2CrO4 with the K2NiF4-type structure[END_REF][START_REF] Kim | Phase-Sensitive Observation of a Spin-Orbital Mott State in Sr2IrO4[END_REF][START_REF] Haskel | Pressure Tuning of the Spin-Orbit Coupled Ground State in Sr2IrO4[END_REF][START_REF] Matsuno | Metallic Ferromagnet with Square-Lattice CoO2 Sheets[END_REF][START_REF] Wang | Synthesis, structures, and magnetic properties of novel Ruddlesden-Popper homologous series Srn+1ConO3n+1 (n = 1, 2, 3, 4, and )[END_REF][START_REF] Luke | Time-reversal symmetry-breaking superconductivity in Sr2RuO4[END_REF][START_REF] Xia | High resolution polar Kerr effect measurements of Sr2RuO4: evidence for broken time-reversal symmetry in the superconducting state[END_REF][START_REF] Pustogow | Constraints on the superconducting order parameter in Sr2RuO4 from oxygen-17 nuclear magnetic resonance[END_REF][START_REF] Zhou | Orbital-Ordering Transition in Sr2VO4[END_REF][START_REF] Teyssier | Optical study of the phonons and electronic excitations in tetragonal Sr2VO4[END_REF][START_REF] Sugiyama | Hidden magnetic order in Sr2VO4 clarified with µ + SR[END_REF][START_REF] Karmakar | Novel P-T Phase Diagram of the Multiorbital Mott Insulator Sr2VO4[END_REF][START_REF] Teyssier | Teller induced nematic orbital order in tetragonal Sr2VO4[END_REF] Several of these compounds are Mott insulators. [START_REF] Matsuno | Variation of the Electronic Structure in Systematically Synthesized Sr2MO4 (M = Ti, V, Cr, Mn, and Co)[END_REF][START_REF] Weng | Electronic structure and optical properties of layered perovskite Sr2MO4 (M = Ti, V, Cr and Mn): An Ab initio study[END_REF][START_REF] Nakamura | eg-level splitting in a layered perovskite manganite as revealed by charge modulation spectroscopy[END_REF][START_REF] Sakurai | Synthesis Conditions and Magnetic Properties of Sr2CrO4 with the K2NiF4-type structure[END_REF][START_REF] Kim | Phase-Sensitive Observation of a Spin-Orbital Mott State in Sr2IrO4[END_REF] More interestingly unconventional superconductivity with broken time-reversal symmetry was discovered in Sr2RuO4 28- 30 and complex 2D spin-orbital ordering was found in Sr2IrO4 [START_REF] Kim | Phase-Sensitive Observation of a Spin-Orbital Mott State in Sr2IrO4[END_REF][START_REF] Haskel | Pressure Tuning of the Spin-Orbit Coupled Ground State in Sr2IrO4[END_REF] and Sr2VO4. [START_REF] Zhou | Orbital-Ordering Transition in Sr2VO4[END_REF][START_REF] Teyssier | Optical study of the phonons and electronic excitations in tetragonal Sr2VO4[END_REF][START_REF] Sugiyama | Hidden magnetic order in Sr2VO4 clarified with µ + SR[END_REF][START_REF] Karmakar | Novel P-T Phase Diagram of the Multiorbital Mott Insulator Sr2VO4[END_REF][START_REF] Teyssier | Teller induced nematic orbital order in tetragonal Sr2VO4[END_REF] In this latter case, the exact nature of this coupling remains to be understood. [START_REF] Zhou | Orbital-Ordering Transition in Sr2VO4[END_REF][START_REF] Teyssier | Optical study of the phonons and electronic excitations in tetragonal Sr2VO4[END_REF][START_REF] Sugiyama | Hidden magnetic order in Sr2VO4 clarified with µ + SR[END_REF][START_REF] Karmakar | Novel P-T Phase Diagram of the Multiorbital Mott Insulator Sr2VO4[END_REF][START_REF] Teyssier | Teller induced nematic orbital order in tetragonal Sr2VO4[END_REF] Despite these interesting electronic properties, the vibrational properties of these compounds have been poorly investigated and are far from being understood despite their importance in the understanding of the properties of these A2MO4 compounds. The aim of the present work is to combine both theoretical and experimental approaches to study the lattice dynamics of the simplest of these compounds, Sr2TiO4, in which the titanium d-orbitals are empty.

In the past, relatively few works investigating the lattice dynamics of Sr2TiO4 and its related properties have been performed. These mainly focused on the study of phonons at the Brillouin zone center by means of infrared-spectroscopy [START_REF] Burns | Raman and infrared studies of Sr2TiO4: A material isomorphic to (La,Sr)2CuO4 superconductors[END_REF][START_REF] Kamba | Composition dependence of the lattice vibrations in Srn+1TinO3n+1 Ruddlesden-Popper homologous series[END_REF][START_REF] Goian | Influence of the central mode and soft phonon on the microwave dielectric loss near the strain-induced ferroelectric phase transitions in Srn+1TinO3n+1[END_REF] and Raman spectroscopy. [START_REF] Burns | Raman and infrared studies of Sr2TiO4: A material isomorphic to (La,Sr)2CuO4 superconductors[END_REF][START_REF] Kamba | Composition dependence of the lattice vibrations in Srn+1TinO3n+1 Ruddlesden-Popper homologous series[END_REF][START_REF] Ge | Sol-gel synthesis and dielectric properties of Ruddlesden-Popper phase Srn+1TinO3n+1 (n = 1, 2, 3, )[END_REF][START_REF] Burns | Raman measurement of materials with K2NiF4 structure[END_REF][START_REF] Venkastewaran | Pressure dependence of the Raman modes in Sr2TiO4[END_REF] The phonon dispersion curves and density of phonon states have only been recently reported using DFT calculations. [START_REF] Fennie | Structural and dielectric properties of Sr2TiO4 from first principles[END_REF][START_REF] Huang | Structure Dependent Phase Stability and Thermal Expansion of Ruddlesden-Popper Strontium Titanates[END_REF][START_REF] Parker | Strain-induced incommensurate distortions in eptixial Ruddlesden-Popper-type Ba2TiO4[END_REF] Several additional DFT works investigated the elastic constants of Sr2TiO4 using stress-strain methods, although discrepancies of ~10 % were found between those studies. [START_REF] Music | Elastic properties of Srn+1TinO3n+1 phases (n=1-3,)[END_REF][START_REF] Lu | First-principles study on the elastic properties of Sr-Ti-O ceramics[END_REF][START_REF] Lu | First-principles calculations of the thermal properties of SrTiO3 and SrO(SrTiO3)n (n=1,2)[END_REF]

II. Experimental details

Ceramic samples of Sr2TiO4 were obtained from a stoichiometric mixture of 99.999 % SrCO3 and 99.99 % TiO2. The mixture was initially ball-milled and mixed for several hours in an agate jar, then calcined at 1200°C for 12 h. In a second step, the obtained product was again ball-milled and mixed over several hours, before sintering for 12 h at 1100°C. The resulting product was a single-phase Sr2TiO4 specimen, as shown by X-ray diffraction (XRD) characterization (see Fig. 1). The room temperature XRD experiment was performed with a PW1820 Philips diffractometer using copper K1 and K2 wavelengths and a 0.02° step. The lattice parameter and atomic positions were derived from a Rietveld refinement of the XRD patterns using the Fullprof software. [START_REF] Rodriguez-Carvajal | Recent advances in magnetic-structure determination by neutron powder diffraction[END_REF] We derive the lattice parameters a = 3.88267(2) Å ; c = 12.5838(1) Å, in good agreement with published results (a = 3.88656 Å ; c = 12.5975 Å [START_REF] Kawamura | Structural Origin of the Anisotropic and Isotropic Thermal Expansion of K2NiF4-Type LaSrAlO4 and Sr2TiO4[END_REF] or a = 3.8812 Å ; c = 12.601 Å [START_REF] Sohn | Cooperative interaction of oxygen octahedra for dielectric properties in the perovskite-related layered compounds Srn+1TinO3n+1, Can+1TinO3n+1 and Srn+1(Ti0.5Sn0.5)nO3n+1 (n = 1, 2, 3 and )[END_REF] ).

High temperature XRD experiments were performed with a Seifert diffractometer using copper K1 and K2 wavelengths and a 0.04° step. The heating was accomplished using a Pt plate as heater and the temperature was measured with a Pt/Rh0.1Pt0.9 S type thermocouple.

Raman scattering experiments under high pressure were performed using a homemade Raman spectrometer equipped with a Peltier cooled CCD and a green laser ( = 514.5 nm). High pressure measurements were performed within a membrane diamond anvil cell (DAC) with low fluorescence diamonds. The samples were placed into a 125 µm chamber drilled in an indented stainless-steel gasket.

The transmitting medium was a methanol:ethanol mixture with a 1:4 ratio, maintaining hydrostatic conditions up to 10 GPa. The fluorescence of ruby was used as a pressure gauge. [START_REF] Piermiarini | Calibration of the pressure dependence of the R1 ruby fluorescence line to 195 kbars[END_REF] Raman scattering experiments below room temperature were performed in a Horiba Jobin-Yvon T64000 spectrometer equipped with a liquid N2 cooled CCD and an Ar + laser ( = 514.5 nm). The low temperature measurements were performed in an Oxford Helium microstat under the objective lens (X50) of an Olympus microscope. The temperature was measured with a Si diode.

Raman scattering experiments above room temperature were performed in a Horiba Jobin-Yvon Labram Aramis spectrometer equipped with a Peltier cooled CCD and a blue diode laser ( = 473 nm). High temperatures were attained using the oven of a Linkham TS1500 heating stage under the objective lens (X50) of an Olympus microscope. The temperature was measured with a Pt resistance thermometer at the bottom of the oven.

In all the Raman experiments, the spectra were obtained in a backscattering configuration using a microprobe device which allowed the incident light to be focused into a spot of diameter several microns on the sample.

The inelastic neutron scattering (INS) experiments were performed at room temperature with a powder sample using the high-resolution IN6 time-of-flight (TOF) spectrometer at the Institute Laue Langevin in Grenoble, France. An incident neutron wavelength of 4.14 Å was utilized, offering an elastic resolution of about 0.2 meV. Standard corrections for background, sample holder, detector efficiency and self-attenuation have been applied to the TOF signal. The corrected response has been transformed to the dynamic structure factor by established computational routines within the software package LAMP. The generalized density of states (GDOS) has been computed from the dynamic structure factor within the incoherent approximation. [START_REF] Koza | Breakdown of phonon glass paradigm in La-and Ce-Fe4Sb12 skutterudites[END_REF] Heat capacity experiments were performed in a Physical Property Measurement System (PPMS) cryostat from Quantum Design, using a short relaxation method.

III. Computational details

The first-principles calculations were performed by means of density functional theory (DFT) and the local density approximation (LDA) using the ABINIT package. [START_REF] Gonze | ABINIT: First-principles approach to material and nanosystem properties[END_REF] We used highly transferable normconserving Teter pseudopotentials. [START_REF] Teter | Additional condition for transferability in pseudopotentials[END_REF] Strontium 4s, 4p and 5s electrons, Ti 3s, 3p, 3d and 4s electrons and O 2s and 2p electrons were treated as valence states.

The plane wave energy cutoff was 60 Ha. Brillouin zone integration was performed using a 8x8x8

Monkhorst-Pack mesh for the calculation of the zone-center phonon modes, while a 12x12x12 mesh was used to compute the formation energy and the electronic band gap. The structural relaxation was performed in the primitive body-centered tetragonal unit cell, until the maximum residual forces on each atom were less than 6x10 -6 Ha/Bohr and the pressure was less than 2x10 -8 Ha/Bohr 3 . From our DFT calculations with the LDA exchange-correlation functional, we find lattice parameters a = 3.8217 Å and c = 12.3145 Å, which are about 2 % lower than our experimental results and prior results in the literature 49,51 . In contrast, our c/a ratios are almost the same (c/a = 3.2223 in DFT against c/a = 3.241 in our experiments and 3.241 [START_REF] Kawamura | Structural Origin of the Anisotropic and Isotropic Thermal Expansion of K2NiF4-Type LaSrAlO4 and Sr2TiO4[END_REF] , 3.247 [START_REF] Sohn | Cooperative interaction of oxygen octahedra for dielectric properties in the perovskite-related layered compounds Srn+1TinO3n+1, Can+1TinO3n+1 and Srn+1(Ti0.5Sn0.5)nO3n+1 (n = 1, 2, 3 and )[END_REF] in prior experiments). Our results agree well with prior LDA calculations performed by Fennie and Rabe. [START_REF] Fennie | Structural and dielectric properties of Sr2TiO4 from first principles[END_REF] The dynamical matrix, dielectric constants, Born effective charges and elastic constants have been calculated using a variational approach to density functional perturbation theory. [START_REF] Gonze | Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory[END_REF] Phonon dispersion curves are interpolated according to the scheme described by Gonze et al. [START_REF] Gonze | Interatomic force constants from densityfunctional perturbation theory: the case of -quartz[END_REF] In this method, the dipole-dipole interactions are subtracted from the dynamical matrix before the Fourier transformation, so that only the short-range part is handled in the real space. A 4x4x4 q-grid in the irreducible Brillouin zone is employed for the calculations of the vibrational dispersion curves, while a denser 90x90x60 q-grid is used for the calculation of the phonon density of states and the thermodynamic functions. The Raman susceptibility tensors are obtained within a non-linear formalism using the 2n+1 theorem. [START_REF] Veithen | First Principles Study of Electro-optic Effect in Ferroelectric Oxides[END_REF] The infrared transmittance, infrared reflectivity and Raman spectra are calculated as described in refs. [START_REF] Hermet | Far-and Mid-Infrared of Crystalline 2,2'-Bithiophene; Ab Initio Analysis and Comparison with Infrared Response[END_REF][START_REF] Hermet | Magnetic, and Phonon Properties of Nickel Hydroxide[END_REF][START_REF] Hermet | First-Principles Calculations of the Non-Linear Optical Susceptibilities and Raman Scattering Spectra of Lithium Niobate[END_REF] . The isothermal mode Grüneisen parameters at the Brillouin zone center have been calculated via a finite differences method using 𝛾 𝑖 𝑇 = 𝐵 𝜔 𝑖 ∆𝜔 𝑖 ∆𝑃 , by computing the frequencies at two additional different pressure P = -1 GPa and 1 GPa, reoptimizing the lattice parameters and the atomic positions, then using the bulk modulus B obtained by DFPT calculations.

IV. Results and discussion

1) Lattice dynamics a) Acoustic phonons

We report in Table 1 the DFT results of the elastic constants Cij. From these Cij, we have calculated the bulk modulus Bh, the shear modulus Gh, the Young's modulus Eh and Poisson's coefficient h within the Hill approximation. [START_REF] Wu | Crystal structure and elastic properties of superhard IrN2 and IrN3 from first principles[END_REF] Our results differ from two recent works where a Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional was used. [START_REF] Music | Elastic properties of Srn+1TinO3n+1 phases (n=1-3,)[END_REF][START_REF] Lu | First-principles study on the elastic properties of Sr-Ti-O ceramics[END_REF] This is certainly related to the use of different exchangecorrelation functionals. Indeed, using the LDA exchange-correlation functional, we find larger elastic constants due to the smaller lattice parameters in our calculations. We expect that the experimental value will lie between our results and those obtained with a PBE exchange-correlations functional. We find a ratio Bh/Gh = 1.778, which is larger than the value of 1.75 proposed for the boundary between brittle and ductile materials. [START_REF] Pugh | Relations between the elastic moduli and the plastic properties of polycrystalline pure metals[END_REF] This means that Sr2TiO4 is slightly ductile and lies close to this boundary. This is surprising because ductile materials are generally metallic, whereas brittle materials are semiconducting.

63,64

We find a powder-averaged Poisson's coefficient h = 0.263, a value similar to those found in ZnSb or GaSb and typical of materials whose bonding is more covalent than ionic. [START_REF] Jund | Physical properties of thermoelectric zinc antimonide using first-principles calculations[END_REF] The knowledge of the Cauchy pressures also gives an insight into the bonding character along the different crystal directions of a material. [START_REF] Jund | Physical properties of thermoelectric zinc antimonide using first-principles calculations[END_REF] The Cauchy pressures in the ab plane and along the c axis can be defined as Pxy = C23 -C44 and Pz = C12 -C66, respectively. In the case of Sr2TiO4, one finds Pxy = 18.3 GPa and Pz = -14.5

GPa. Thus, the rather large and positive value of the Cauchy pressure in the ab plane is more characteristic of ionic compounds, whereas the rather large and negative value of the Cauchy pressure in the c direction is characteristic of covalent materials. [START_REF] Jund | Physical properties of thermoelectric zinc antimonide using first-principles calculations[END_REF] This means that the bonds are more directional and have a more covalent character in the c direction than in the ab plane. Despite this bonding anisotropy, one can determine [START_REF] Saha | Elastic and thermal properties of the layered thermoelectrics BiOCuSe and LaOCuSe[END_REF] that the anisotropy of the bulk modulus

𝐵 𝑎 𝐵 𝑐 = (𝐶 11 -𝐶 33 )+(𝐶 12 -𝐶 13 ) (𝐶 33 -𝐶 13 )
is as low as 1.021, and therefore that Sr2TiO4 has a quasi-isotropic bulk modulus.

b) Optical phonons at the Brillouin zone center

There are 7 atoms in the primitive unit cell and therefore 21 phonon modes, including 3 acoustic and 18 optical phonon modes. The decomposition into irreducible representations of the optical modes is:

 = 2 A1g  2 Eg  3 A2u  4 Eu  B2u.
We report in Table 2 our results for the phonon energies at the  point. There are 4 Raman-active modes with A1g and Eg symmetries, and 7 infrared-active modes with A2u and Eu symmetries. Our calculations are in very good agreement with those of Fennie and Rabe for the fully-relaxed cell, also based on the LDA exchange correlation. [START_REF] Fennie | Structural and dielectric properties of Sr2TiO4 from first principles[END_REF] We report in Figure 2 the experimental Raman spectrum of Sr2TiO4 at 20 K, together with the DFPTcalculated spectrum. There is a good agreement between experiments and calculations for both the positions and the relative intensities of the Raman lines located at 125 cm -1 , 210 cm -1 , 280 cm -1 and 570 cm -1 . Table 2 shows the exact frequencies of the Raman modes in both experiments and calculations.

The broader lines observed at 400-450 cm -1 and 700 cm -1 have been previously reported and are probably due to some second order Raman modes or defect-induced excitations. [START_REF] Burns | Raman and infrared studies of Sr2TiO4: A material isomorphic to (La,Sr)2CuO4 superconductors[END_REF][START_REF] Kamba | Composition dependence of the lattice vibrations in Srn+1TinO3n+1 Ruddlesden-Popper homologous series[END_REF][START_REF] Burns | Raman measurement of materials with K2NiF4 structure[END_REF][START_REF] Venkastewaran | Pressure dependence of the Raman modes in Sr2TiO4[END_REF] In the latter case, given the high frequencies of the modes and the broadness of these peaks, they are more likely to originate from defects on the O sub-lattice. Concerning the infrared-active modes, there is a reasonable agreement between calculations and experiments [START_REF] Kamba | Composition dependence of the lattice vibrations in Srn+1TinO3n+1 Ruddlesden-Popper homologous series[END_REF][START_REF] Burns | Raman measurement of materials with K2NiF4 structure[END_REF] for the four bands below 300 cm -1 , but a rather poor agreement concerning the three modes above 300 cm -1 . However, for these high-frequency modes, even the two experimental reports yield inconsistent frequencies. [START_REF] Kamba | Composition dependence of the lattice vibrations in Srn+1TinO3n+1 Ruddlesden-Popper homologous series[END_REF][START_REF] Burns | Raman measurement of materials with K2NiF4 structure[END_REF] c) Lattice dynamics in the full Brillouin zone

In Figure 3, we report the phonon dispersion curves, as well as the total and atom-projected phonon densities of states (DOS) of Sr2TiO4. The shoulder at 95 cm -1 is due to the flattening of the transverse acoustic modes at the Brillouin zone boundaries. The peak at 120 cm -1 is due to the flattening of the longitudinal acoustic mode at the Brillouin zone boundaries as well as the lower energy optical mode with Eg symmetry at the  point, which is a Raman-active vibration in the ab plane. The next two peaks in the phonon DOS at 160 and 175 cm -1 are due to the next-lowest energy optical mode with Eu symmetry, which is an infrared-active vibration in the ab plane. The subsequent features in the phonon DOS up to 575 cm -1 are due to many different optical branches. The last two features in the phonon DOS originate from the higher energy optical mode which is an infrared-active vibration in the ab plane,

with the peak at about 850 cm -1 caused by the dispersion of the LO mode. The phonon DOS is dominated below 140 cm -1 by Sr vibrations, between 140 and 210 cm -1 by Sr and O, between 210 and 570 cm -1 by O and Ti, and almost exclusively by O above 570 cm -1 .

When comparing our results with those of Huang et al, [START_REF] Huang | Structure Dependent Phase Stability and Thermal Expansion of Ruddlesden-Popper Strontium Titanates[END_REF] we find significant differences, especially at low frequencies. Notably, they found an optical mode at about 2 THz, i. e. about 67 cm -1 , in contrast with our calculations and Fennie's work [START_REF] Fennie | Structural and dielectric properties of Sr2TiO4 from first principles[END_REF] in which the lowest energy optical mode is at about 124 cm - The GDOS and Z'(E) are shown in Figure 4. There is a qualitative agreement up to 550 cm -1 and especially below 300 cm -1 , with a good correspondence between the positions of the different peaks and shoulders in the spectral densities. However, above 550 cm -1 , the agreement is quite poor. This may be due to the presence of O sublattice defects in the samples, as suggested by the additional broad bands in the Raman spectra, which are likely associated with the presence of defects.

d) Heat capacity

The heat capacity of Sr2TiO4 is reported in Figure 5, where our experimental results below 300 K together with Jacobs' results above 300 K 66 are compared with our DFPT calculations. The agreement is rather good up to 650 K. Above this temperature the experimental heat capacity is larger than the calculated one and increases linearly with temperature, whereas the calculated heat capacity saturates to the Dulong-Petit limit of 174.6 J/mole.K. Indeed, in the DFPT calculations, we calculate the heat capacity Cv at constant volume, whereas the experimental data correspond to the heat capacity CP at constant pressure. The difference between the two terms CP -CV is equal to V 2 BVT. Unfortunately, although the thermal expansion has been measured (see ref. 70 and below), this is not the case for the bulk modulus. With B obtained by DFPT calculations, we compute CP -CV = 26.9 J/mole.K at 1600 K. This is the same order of magnitude as the difference between the experimental and calculated data (12.94 J/mole.K).

2) Raman scattering under pressure

The Raman spectra of Sr2TiO4 as a function of pressure are shown in Figure 6. In contrast to the study reported in ref. [START_REF] Venkastewaran | Pressure dependence of the Raman modes in Sr2TiO4[END_REF] , which was limited to 7 GPa, we can observe Raman spectra with good signal-to-noise ratio up to the maximum pressure of 31 GPa. There is a broadening of the Raman peaks with increasing pressure. Above roughly 11-13 GPa, we observe the appearance of a new peak at about 200 cm -1 , a shoulder at about 660 cm -1 , and the decay of the small peak from the lowest energy Raman mode Eg at about 127 cm -1 . Above 14.5 GPa, a third new weak peak appears at 157 cm -1 . We also note that the broad peak at 520-570 cm -1 increases in intensity at high pressure. One cannot exclude that this peak is due to a new, first-order Raman mode. These new peaks could be related to a structural distortion of the tetragonal structure above 11 GPa. In a recent theoretical study, a phase transition was found upon compression from the I4/mmm structure stable at ambient conditions to the Cmcm structure, although at much higher pressure (66.6 GPa). [START_REF] Gao | NaCl-CsCl structural transition in Sr2PdO3 and Sr2TiO4[END_REF] However, the observed modification shows little change in the Raman spectra and may instead be related to a distortion associated with non-hydrostatic compression.

The pressure variation of the different mode frequencies is linear at low pressures and can be described by a 2 nd order polynomial function, similarly to the case of silicon. [START_REF] Ulrich | Lifetimes of Phonons in Semiconductors under Pressure[END_REF] The pressure derivatives d/dP of the four Raman-active modes derived from polynomial fits are reported in Table 4. These values and their trend compare well with data obtained by Venkastewan, [START_REF] Venkastewaran | Pressure dependence of the Raman modes in Sr2TiO4[END_REF] except for the highest energy A1g mode, for which we find a 20 % larger value. From the knowledge of d/dP, it is possible to determine the mode Grüneisen parameter. This will be discussed below.

3) Temperature dependent properties a) Raman scattering

We report in Figure 7 several selected Raman spectra between ambient temperature and 1473 K. The thermal variation of the frequencies and the full-width at half-maximum (FWHM) of the four Ramanactive modes are reported in Figures 8 and9. For the three lowest energy modes, a monotonic decrease of the peak positions with increasing temperature is visible. The behavior of the highest energy mode is more complex, with an upshift of the peak position from the lowest temperature up to about 200 K where it reaches a maximum, followed by a downshift above 200 K. We will quantitatively discuss the different contributions below.

b) Thermal expansion and the thermodynamic Grüneisen parameter

The thermal variation of the lattice parameters a and c of Sr2TiO4 is shown in Figure 10. The error bars of the lattice parameters are within the symbol size. From calibration measurements, the error of the temperature measurements has been estimated to be about 20 K at 900 K. One can see that at higher temperatures (from 1000 to 1400 K), the increase of the lattice parameters with temperature is stronger than linear. The ratio c/a is 3.239 at room temperature and increases slightly to 3.241 between 1000 K and 1200 K, but stays constant above 1200 K. Kawamura et al have determined the thermal expansion coefficients (x = (da/dT)/a, z = (dc/dT)/c and V = 2 x + z) from the change of the lattice parameters between 300 K and 1300 K and found x = 12.25 MK -1 , z = 12.7 MK -1 and V = 37.2 MK -1 . [START_REF] Kawamura | Structural Origin of the Anisotropic and Isotropic Thermal Expansion of K2NiF4-Type LaSrAlO4 and Sr2TiO4[END_REF] Applying the same procedure we find slightly lower values: x = 12.09 MK -1 , z = 12.34 MK -1 and V = 36.52 MK -1 . These reinforce Kawamura et al.'s finding that the thermal expansion is weakly anisotropic in Sr2TiO4 with x < z. Parameters derived from polynomial fits to the thermal variation of the thermal expansion coefficients x, z and V are reported in Figure 11. An estimation of the error bar is given at 900 K. One can see that the thermal expansion coefficients increase with temperature and are rather large. Between room temperature and 1700 K, z increases by about 45 %, which is much higher than the approximate 25 % increase in x. Over the same temperature range, V increases from 32.2 MK -1 to 42.2 MK -1 , which is higher than in SrTiO3 where an average of 32.2 MK -1 is reported between 300 and 2000 K. [START_REF] De Ligny | High temperature heat capacity and thermal expansion of SrTiO3 and SrZrO3 perovskites[END_REF] To shed light on the anharmonicity of Sr2TiO4, we compute the thermodynamic Grüneisen parameter  = VB/VCP 70 from the experimentally-determined parameters V, molar volume V, and heat capacity CP as well as the DFPT-calculated bulk modulus B (see Table 1). The results are reported in the inset of Figure 12. At room temperature,  is about 2. Note that using the bulk modulus B of Music et al, [START_REF] Hermet | Magnetic, and Phonon Properties of Nickel Hydroxide[END_REF] one derives a value of 1.75.  is increasing slightly at the highest temperature. The thermodynamic Grüneisen parameter  is rather large, meaning that the anharmonicity is significant in Sr2TiO4. We have also determined the thermodynamic Grüneisen parameters in the x and z directions, x and z as follows:

using elastic constants Cij calculated by DFPT. As can be seen in the inset of Figure 12, the anisotropy of the Grüneisen parameters is weak in Sr2TiO4, meaning that the anharmonicity is similar in both the ab plane and the c direction.

4) Mode Grüneisen parameters and anharmonicity

In the harmonic approximation, the frequencies of the vibrational modes are constants and do not change with volume, pressure or temperature. In the quasi-harmonic approximation, the frequencies of the vibrational modes vary upon volume change and therefore also upon temperature and pressure variations.

The thermal variation of the mode frequencies i is divided into two contributions, yielding: [START_REF] Lucazeau | Effect of pressure and temperature on Raman spectra of solids: anharmonicity[END_REF][START_REF] Liarokiapis | Raman study of the phonon anharmonicity in LaF3[END_REF] i(T) = i(T = 0 K) + i impl + i expl (1)

Here i impl is the implicit contribution accounting for the volume dependence of the mode frequencies, and i expl is the explicit contribution determined by the thermal population of the vibrational levels which depends on the anharmonicity.

In the isotropic Grüneisen approximation, one can define the isothermal Grüneisen parameter 𝛾 𝑖 𝑇 and the isobaric Grüneisen parameter 𝛾 𝑖 𝑃 of the i th mode as: [START_REF] Lucazeau | Effect of pressure and temperature on Raman spectra of solids: anharmonicity[END_REF][START_REF] Liarokiapis | Raman study of the phonon anharmonicity in LaF3[END_REF] 𝛾 𝑖 𝑇 = 𝐵 (

The isobaric Grüneisen parameter can be expressed as: [START_REF] Liarokiapis | Raman study of the phonon anharmonicity in LaF3[END_REF] 𝛾 𝑖 𝑃 = 𝛾 𝑖 𝑇 + 𝛾 𝑖 𝑉

One can apply the isotropic approximation because the bulk modulus is almost isotropic, as is the thermal expansion. Using Raman scattering experiments under high pressure, one can determine the isothermal Grüneisen parameters 𝛾 𝑖 𝑇 ; while using Raman experiments as a function of temperature, one can determine the isobaric Grüneisen parameters 𝛾 𝑖 𝑃 . When the quasi-harmonic approximation can be strictly applied, 𝛾 𝑖 𝑉 = 0 and 𝛾 𝑖 𝑃 = 𝛾 𝑖 𝑇 . Therefore, in the general case, 𝛾 𝑖 𝑉 is a measure of the anharmonic contribution to the frequency shift, i.e. the explicit effect, whereas 𝛾 𝑖 𝑇 corresponds to the volume contribution, i.e. the implicit effect. One can thus define the implicit fraction i = 𝛾 𝑖 𝑇 /𝛾 𝑖 𝑃 and evaluate the magnitude of the implicit and explicit contributions to each mode.

Using the bulk modulus computed by DFPT, one can determine the isothermal Grüneisen parameters 𝛾 𝑖 𝑇 of the different Raman-active modes, which are given in Table 3. We have also performed DFPT calculations of the isothermal Grüneisen parameters 𝛾 𝑖 𝑇 of the different modes at the Brillouin zone center. The theoretical results are given in Table 4. In the case of the Raman-active modes (see Table 3), there are some significant disagreements for the lowest energy Eg and A1g modes (about 15-20 %), but a good agreement for the highest energy Eg and A1g modes (about 3-5 %). Moreover, the theoretical calculations provide a good reproduction of the experimental trend. The values of the two lowest Ramanactive energy modes (involving mainly Sr motions) are larger than the value of the thermodynamic Grüneisen parameter  and about two times larger than those of the two highest energy Raman-active modes (for which the contribution of the oxygen atoms in the SrO layers is predominant). This means that the two lowest energy Raman-active modes are more anharmonic than the two highest energy Raman-active modes, and thus that vibration modes involving Sr atoms are the most anharmonic ones, contributing much more strongly to the thermal expansion. In contrast, the optically silent mode B2u

(which involves only c-axis motion of the oxygen atoms in the TiO2 layer) has a very low Grüneisen parameter, and therefore very weak anharmonicity.

In Fig. 8, for the three lowest energy Raman-active modes one observes a monotonic decrease of the peak positions with increasing temperature. We note that the Raman shift of the highest Eg mode stays constant from the lowest temperature up to 300 K, which indicates that the explicit and implicit contributions have opposite signs, cancelling out each other in this temperature range. The peak position of the highest energy A1g mode upshifts from the lowest temperature up to about 200 K where it reaches a maximum, followed by a downshift above 200 K. The behavior of this mode can be explained by the opposite signs of the implicit and explicit contributions. As the isothermal Grüneisen parameter 𝛾 𝑖 𝑇 of the highest energy mode is positive, the implicit volume term i impl is negative, which means that the explicit anharmonic term i expl should be positive and even the largest contribution, at least up to 200 K.

Since the thermal variation of the different modes is linear from room temperature until at least 800 K, one can determine the isobaric Grüneisen parameter 𝛾 𝑖 𝑃 from the slope obtained by a linear fit of the thermal variation of the frequencies (with the knowledge of the thermal expansion). Now, knowing both 𝛾 𝑖 𝑃 and 𝛾 𝑖 𝑇 , one can deduce 𝛾 𝑖 𝑉 and i of each Raman-active mode. The results are given in Table 3. All modes 𝛾 𝑖 𝑉 have negative values, meaning that the explicit anharmonic terms i expl have positive values for the 300-800 K temperature range. In this temperature range, the ratio i is much larger than 0.5, which means that the implicit volume contribution is the dominant contribution in the temperature dependence of the frequencies of the Raman-active modes of Sr2TiO4. The large values of the ratio i of the Raman-active modes found in Sr2TiO4 are typical of ionic materials, e. g. alkaline-earth fluorides 74-76 , a conclusion different from that which we obtain from the Poisson coefficient. However, as discussed above, the explicit anharmonic contribution is dominant at least up to 200 K for the highest energy mode. This is due to the low Grüneisen parameter of this mode. Both the large explicit contribution and the low mode Grüneisen parameter of this A1g mode (involving c-axis O motion)

indicate that bonding in this direction is less ionic/more covalent. This observation is in good agreement with our previous analysis of the Cauchy pressure, indicating less ionic/more covalent bonding in the c direction.

In order to better account for the thermal variation of the implicit and explicit contributions, we assume that the Grüneisen parameters are independent of the temperature and correspond to the isothermal Grüneisen parameters 𝛾 𝑖 𝑇 . One can then determine both contributions (in the isotropic approximation) as follows: [START_REF] Lucazeau | Effect of pressure and temperature on Raman spectra of solids: anharmonicity[END_REF][START_REF] Liarokiapis | Raman study of the phonon anharmonicity in LaF3[END_REF][START_REF] Balkanski | Anharmonic effects in light scattering due to optical phonons in silicon[END_REF] 𝜔 𝑖 𝑖𝑚𝑝𝑙 = 𝜔 𝑖 (𝑇 = 0𝐾). 𝑒 -𝛾 𝑖 𝑇 ∫ 𝛼 𝑉 (𝑇)𝑑𝑇 𝑇 0

If we assume that the volume thermal expansion V is independent of the temperature, we obtain:

𝜔 𝑖 𝑖𝑚𝑝𝑙 = 𝜔 𝑖 (𝑇 = 0𝐾). 𝑒 -𝛾 𝑖 𝑇 𝛼 𝑉 𝑇 (5)
𝜔 𝑖 𝑒𝑥𝑝𝑙 = 𝐴 (1 + 2 𝑒 𝑥 -1 ) + 𝐵 [1 + 3 𝑒 𝑦 -1 + 3 (𝑒 𝑦 -1) 2 ] (6)
where 𝑥 = ħ𝜔 2𝑘 𝐵 𝑇 and 𝑦 = ħ𝜔 3𝑘 𝐵 𝑇 .

Using the constant value of the volume thermal expansion determined between 300 K and 1273 K, V = 36.52 MK -1 , one can determine the implicit quasi-harmonic contribution from eq. ( 5) and hence deduce the explicit anharmonic contribution. The results are given in Figure 12, together with the fitting of the explicit contribution with eq. ( 6). We note that for all the modes, A > 0 and B < 0.

The FWHM i of the Raman modes can be expressed as a function of the residual term (0), the 3phonon interaction term 3 and the 4-phonon interaction term 4: [START_REF] Lucazeau | Effect of pressure and temperature on Raman spectra of solids: anharmonicity[END_REF][START_REF] Liarokiapis | Raman study of the phonon anharmonicity in LaF3[END_REF][START_REF] Balkanski | Anharmonic effects in light scattering due to optical phonons in silicon[END_REF] i = (T = 0 K) + 3 + 4

where the term ( T = 0 K) is constant, and 88,89,93

Δ 3 = 𝐶 (1 + 2 𝑒 𝑥 -1
) and

Δ 4 = 𝐷 [1 + 3 𝑒 𝑦 -1 + 3 (𝑒 𝑦 -1) 2 ] (8)
We show the results of modeling the FWHM i using eqs. (7-8) in Figure 9. It is difficult to draw quantitative conclusions for the lowest Raman-active mode since the FWHM have large errors, especially at high temperature. However, our data indicate that 𝐶 and 𝐷 are positive for the three lowest energy modes, whereas 𝐶 > 0 and 𝐷 < 0 for the highest energy mode.

Conclusion

We have reported a combined experimental and theoretical study of the lattice dynamics in the Ruddlesden-Popper layered compound Sr2TiO4. We have measured the generalized phonon density of states of Sr2TiO4 using inelastic neutron scattering experiments as well as the pressure and temperature dependence of the Raman-active modes of Sr2TiO4. We have determined a rather large thermodynamic

Grüneisen parameter of about 2 from the heat capacity and thermal expansion experiments, using the calculated bulk modulus. From our Raman scattering experiments under pressure, we report a probable structural distortion of the tetragonal structure above 11 GPa which could be due to non-hydrostatic compression. There is a reasonable agreement between the experimental mode Grüneisen parameters of the Raman-active modes with those obtained by DFPT calculations. The mode Grüneisen parameters of the two lowest energy Raman-active modes involving mainly Sr atoms are two times larger than the two highest energy Raman-active modes. By determining the pressure and temperature dependence of the Raman-active modes, we were able to determine their implicit volume contributions and explicit anharmonic contributions. Above 400 K, the implicit volume contribution dominates the temperature behavior of all four Raman-active modes. This is typical of ionic materials, but contradicts the conclusion drawn from the Poisson coefficient. However, our analysis of the Cauchy pressures indicates that the bondings are more ionic in the ab plane than in the c direction. Therefore, one can consider Sr2TiO4 to exhibit iono-covalent bondings. Below 400 K, the explicit anharmonic contribution is dominant for the highest energy mode, which mainly involves motion of the oxygen atoms along the c direction. This means that these bondings are less ionic/more covalent, in agreement with our analysis of the Cauchy pressures. experiments. [START_REF] Koza | Breakdown of phonon glass paradigm in La-and Ce-Fe4Sb12 skutterudites[END_REF][START_REF] Gonze | ABINIT: First-principles approach to material and nanosystem properties[END_REF][START_REF] Teter | Additional condition for transferability in pseudopotentials[END_REF] Note that the LO modes are labelled in parentheses. In the case of Kamba's data [START_REF] Koza | Breakdown of phonon glass paradigm in La-and Ce-Fe4Sb12 skutterudites[END_REF] , the modes labelled with * may not be true modes.

1 .

 1 For comparing our calculations with the experimental Generalized Density of States (GDOS) obtained from Inelastic Neutron Scattering (INS) experiments, we calculate the weighted neutron density of states 𝑍'(𝐸) using : 𝑍 ′ (𝐸) = ∑ 𝜎 𝑖 𝑀 𝑖 𝑍 𝑖 (𝐸) 𝑖 , where 𝑍 𝑖 (𝐸) are the atom-projected density of states, σi are the scattering cross sections and Mi are the atomic masses of each atom.

Figure captions Figure 1 :

 captions1 Figure captions

Figure 2 :

 2 Figure 2: Experimental (20K) and calculated (0K) Raman spectra.

Figure 3 :

 3 Figure 3: Calculated phonon dispersion curves, total and atom-projected density of states of Sr2TiO4.

Figure 4 :

 4 Figure 4: Experimental generalized density of states of Sr2TiO4 from INS experiments and total phonon

Figure 5 :

 5 Figure 5: Experimental heat capacity CP of Sr2TiO4 from the present work as well as Jacobs and Rajitha 89 , and computed heat capacity CV of Sr2TiO4 obtained from DFPT calculations.

Figure 6 :

 6 Figure 6: Evolution of the Raman spectra with pressure (left). The asterisks correspond to the new

Figure 7 :

 7 Figure 7: Evolution of the Raman spectra with temperature.

Figure 8 :

 8 Figure 8: Variation of the frequencies 1 (Eg), 2 (A1g), 3 (Eg) and 4 (A1g) of the Raman-active modes

Figure 9 :

 9 Figure 9: Variation of the FWMH 1, 2, 3 and 4 of the Raman-active modes with temperature.

Figure 10 :

 10 Figure 10: Thermal variation of the lattice parameters a and c of Sr2TiO4.

Figure 11 :

 11 Figure 11: Thermal evolution of the thermal expansion coefficients x and z (top) and the volume

Figure 12 :FiguresFigure 1

 121 Figure 12: Variation of 1 (Eg), 2 (A1g), 3 (Eg) and 4 (A1g) of the Raman-active modes with

Table captions Table 1 :

 captions1 Calculated mechanical constants Cij, the bulk modulus Bh, shear modulus Gh, Young's modulus Eh (in GPa) and Poisson coefficient h within the Hill approximation in Sr2TiO4.

Table 2 :

 2 Calculated frequencies of the optical phonons (in cm -1 ) compared to Raman and Infrared (IR)

Table 3 :

 3 Frequencies and Grüneisen parameters of the optical phonons (in cm -1 ) from the Raman scattering experiments.

Table 4 :

 4 Frequencies and Grüneisen parameters of the optical phonons (in cm -1 ) from DFPT calculations.

	Tables
			Elastic constants
			C11	308.8
			C33	287.4
			C12	73.8
			C13	91.1
			C44	72.8
			C66	88.3
			Bh	157.5
			Gh	88.6
			Eh	223.8
			h	0.263
			Table 1
	Symmetry modes	A1g	Eg	A2u	Eu	B2u
	DFT (present)
	DFT 58
	Raman 20 K
	(present)
	Raman RT
	(present)
	Raman, IR RT 56
	Raman, IR RT 53
	IR 10 K, film 54

Table 4

 4


Georgios Koutsoumpas

Anastassios Charitopoulos

Christos Ι Papadopoulos

Michel Fillon

Computational Evaluation of Dynamic Coefficients of Thrust Bearings; Effect of Artificial Texturing on Thermohydrodynamic Performance

Keywords: THD, Thrust bearing, Stiffness coefficient, Damping Coefficient, Artificial Surface Texturing, CFD

Fluid film thrust bearings are commonly used in industry, providing durable and reliable operation at high values of load carrying capacity, accompanied by low friction losses. A major advantage of hydrodynamic fluid film bearings, over other types of bearings, is their enhanced dynamic behaviour, especially under transient or impact loads. Currently, a systematic approach to identify the dynamic coefficients of thrust bearing geometrical configurations utilising high complexity CFD simulation data has not yet been established. It is therefore imperative to develop a method, capable of evaluating the dynamic characteristics of complex bearing designs and allow the evaluation of bearing response under transient loads. In the present work, a computational approach is proposed to estimate the stiffness and damping coefficients of fluid-film thrust bearings. A CFD-based ThermoHydroDynamic (THD) numerical model of the bearing is developed and utilised for performing an initial steady-state simulation at given rotational speed and thrust load, as well as subsequent transient simulations at increasing or decreasing thrust loads. The former simulation is used to calculate the stiffness coefficient of the bearing at the specified conditions, while the latter are appropriately post-processed to estimate the damping coefficient of the bearing at different values of rotor acceleration. The procedure is repeated at different operating conditions, yielding a map of the dynamic coefficients of the bearing. Finally, a single degree of freedom model is generated, which utilises the calculated values of dynamic coefficients to evaluate transient bearing performance under any given thrust load history. The proposed methodology is applied to compare the dynamic response characteristics of a conventional sector-pad tapered-land thrust bearing and a textured tapered-land thrust bearing of the same principal dimensions.

Introduction

Today, thrust bearing applications demand high efficiency designs, in order to minimise friction losses. To that end, bearings of smaller principal dimensions are selected, resulting in more heavily loaded designs, which operate at lower values of minimum film thickness. Thus, the dynamic behaviour of such bearings, becomes a significant factor of the design process, in order to avoid contact, due to impact loads. The transient response of thrust bearings under impact loads depend on the stiffness and damping coefficients, which are a function of bearing operating parameters (rotational speed, thrust load, oil temperature). In order to evaluate those coefficients, a systematic method needs to be established.

Recently, due to the increased computing power available, computational methods have proven to be a valid alternative to experimental methods which are expensive and difficult to set-up. Due to the complex computations involved in the physics of the problem, the available literature is quite recent. The two main approaches to computational methods involve either the use of the Reynolds' equations or the full Navier-Stokes (N-S) equations. The Reynolds' equations are derived from the N-S equations by introducing certain assumptions to simplify the physics of the problem. Therefore, they are simpler to solve, but their results are less accurate.

A study conducted by Storteig and White [START_REF] Storteig | Dynamic characteristics of hydrodynamically lubricated fixed-pad thrust-bearings[END_REF] on tapered-land thrust bearings solved the Reynolds' equation using a 1-D finite element technique considering the hot oil carry-over effect to calculate its dynamic coefficients. The results showed a relatively good agreement with coefficients calculated with 2-D methods.

Zhu and Zhang [START_REF] Zhu | A Preliminary Nonlinear Analysis of the Axial Transient Response of the Sector-Shaped Hydrodynamic Thrust Bearing-Rotor System[END_REF] examined the transient axial response of a sector-shaped thrust bearing using a finite difference method to solve the Reynolds' equation. The results highlighted the significance of a nonlinear analysis and that the dynamic coefficients depend strongly on the initial film thickness. Nica [START_REF] Nica | A thermohydrodynamic method of bearing analysis[END_REF] in his work managed to couple the Reynolds equation with the energy equation, and heat conduction equations. He concluded that the effect of the temperature distribution is crucial on the evaluation of bearing performance. In reference [START_REF] Jang | Determination of the dynamic coefficients of the coupled journal and thrust bearings by the perturbation method[END_REF], Jang and Lee presented the two methods that are commonly used to calculate stiffness and damping coefficients. The first, considered as the physical perturbation method, consists in differentiating the forces with respect to finite displacements (for stiffness) and finite velocities (for damping). Regardless of the simplicity of the method, it requires a significant amount of computing time, and the accuracy depends on the perturbation amplitude. The latter is the mathematical perturbation method and uses the finite element method on the perturbed Reynold's equation. In [START_REF] Srikanth | Modelling of Large Tilting Pad Thrust Bearing Stiffness and Damping Coefficients[END_REF], Srikanth et al. used a finite difference method to solve the Reynolds' equation and calculate the pressure in the oil film of tilting pad thrust bearings. Variations of the film thickness and runner axial velocity allow the calculation of stiffness and damping coefficients. Papadopoulos et al. [START_REF] Papadopoulos | Characterisation of stiffness and damping in textured sector pad micro thrust bearings using computational fluid dynamics[END_REF] proposed a computational method to determine the stiffness and damping coefficients of thrust bearings using Computational Fluid Dynamics (CFD). The method acknowledged both translational and tilting motion of the bearing. A hydrodynamic model was set in a steady-state condition (constant rotational speed and minimum film thickness values). Subsequently, a translational or tilting micro-perturbation was forced, altering the axial and tilting velocities. The damping coefficient was calculated by the ratio of the difference of load capacity (for the translational motion) or the average moment (for the tilting motion) to the velocity alteration. Finally, the model was allowed into a new steadystate which allowed the calculation of the stiffness coefficient using the minimum film thickness variations instead of the velocity. The proposed method was applied for smooth and textured pad bearing with results that confirmed published data. In [START_REF] Qin | Prediction of Dynamic Characteristics of Foil Thrust Bearings using CFD[END_REF], the dynamic coefficients of foil thrust bearings were calculated with CFD tools by imposing a sinusoidal motion to the rotor. Snyder and Braun [START_REF] Snyder | Comparison of Perturbed Reynolds' Equation and CFD Models for the Prediction of Dynamic Coefficients of Sliding Bearings[END_REF] studied the difference between transient Reynolds' and N-S equations through a simple slider under harmonic loading. Similar procedures were implemented in [START_REF] Takenaka | Application of CFD analysis for static and dynamic characteristics of hydrodynamic journal bearing[END_REF][START_REF] Li | Stiffness and damping coefficients for journal bearing using the 3D transient flow calculation[END_REF], regarding the dynamic coefficient of journal bearings. A similar approach to that proposed by Papadopoulos et al. [START_REF] Papadopoulos | Characterisation of stiffness and damping in textured sector pad micro thrust bearings using computational fluid dynamics[END_REF] was used in [START_REF] Vieira | Analysis of the influence of force and moment dynamic coefficients on lubricated thrust bearings[END_REF] by Vieira and Cavalca, who solved the Reynolds' equation, instead of the N-S equations, to determine the dynamic coefficients for both axial and tilting motions of a tapered-land thrust bearing.

In the present work, a novel procedure is proposed to calculate the transient response of fluid film thrust bearings using a CFD-based approach. In particular, a Thermohydrodynamic (THD) CFD method is employed to calculate the dynamic coefficients of a tapered-land and a surface textured bearing of the same principal dimensions, operating at the same conditions. Their performance under transient loads is calculated using a single degree of freedom lumped-mass model, which takes into consideration the dynamic coefficients calculated by the CFD-based THD simulations. A CFD-based THD approach was preferred over the simplified Reynolds' equation approach, since the latter does not account for pressure and temperature variations in the cross-flow (film thickness) direction, neglects fluid inertia forces and cannot model the hot-oil carry-over phenomenon between consecutive bearing pads. A THD approach is necessary in order to include heat effects, which are not negligible due to temperature differences and viscous dissipation effects, by additionally solving the energy equation, allowing the local calculation of viscosity to achieve higher precision.

Problem setup 2.1 Calculation of stiffness and damping coefficients

In order to calculate the stiffness and damping coefficients of a fluid film thrust bearing, a Thermohydrodynamic (THD) CFD analysis has been conducted (see Section 2.3 for details). The analysis consists of a steady state simulation followed by a transient one. In the steady state, the bearing is set at a constant minimum film thickness h and the load Fstatic is determined. This permits the calculation of the stiffness coefficient (Hooke's law), eq. (1):

𝑘𝑘(ℎ; 𝜔𝜔) = 𝐹𝐹 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 ℎ (1)
As a result, the stiffness coefficient is based on constant minimum film thickness and does not include any dynamic aspects, which are included in the damping coefficient. The procedure is repeated for various minimum film thickness values in order to generate a stiffness coefficient to minimum film thickness function.

Subsequently, for each minimum film thickness considered in the stiffness coefficient calculation, a perturbation of the minimum film thickness from the previous steady state condition (initial state) is examined using two transient models. In the first model, the rotor is accelerated axially with a constant acceleration, until a specific displacement Δx (namely 0.1 μm) from the initial state is reached, to approach the pad, leading to an increase of the generated force. The second model considers the same acceleration profile forcing the rotor away from the pad, decreasing the generating force. Due to zero initial axial velocity, the axial motion can be described as follows:

ℎ = ℎ| 𝑠𝑠=0 + 1 2 𝑥𝑥̈⋅ 𝑡𝑡 2 (2)
The transient force values Ftrans, resulting from each transient simulation, are considered to be formed by a static part, which has been calculated by the steady-state simulation, and a dynamic part, as it appears in the following equation:

𝐹𝐹 𝑠𝑠𝑡𝑡𝑠𝑠𝑡𝑡𝑠𝑠 (ℎ + 𝛥𝛥𝑥𝑥, 𝑥𝑥, 𝑥𝑥; 𝜔𝜔) = 𝐹𝐹 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 (ℎ + 𝛥𝛥𝑥𝑥; 𝜔𝜔) + 𝐹𝐹 𝑑𝑑𝑑𝑑𝑡𝑡𝑠𝑠𝑑𝑑𝑠𝑠𝑠𝑠 (ℎ + 𝛥𝛥𝑥𝑥, 𝑥𝑥, 𝑥𝑥̈; 𝜔𝜔) (3)
The dynamic part of the total force is used to calculate the dynamic coefficient as it appears in Eq. ( 4) and Eq. ( 5) leading to a damping coefficient for each transient simulation:

𝑐𝑐 𝑑𝑑𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑡𝑡𝑠𝑠𝑠𝑠𝑡𝑡𝑑𝑑 (ℎ, 𝑥𝑥, 𝑥𝑥; 𝜔𝜔) = 𝐹𝐹 𝑑𝑑𝑑𝑑𝑡𝑡𝑠𝑠𝑑𝑑𝑠𝑠𝑠𝑠,𝑑𝑑𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑡𝑡𝑠𝑠𝑠𝑠𝑡𝑡𝑑𝑑 𝑥𝑥̇= 𝐹𝐹 𝑑𝑑𝑑𝑑𝑡𝑡𝑠𝑠𝑑𝑑𝑠𝑠𝑠𝑠,𝑑𝑑𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑡𝑡𝑠𝑠𝑠𝑠𝑡𝑡𝑑𝑑 𝑥𝑥̈⋅ 𝑡𝑡 (4) 𝑐𝑐 𝑠𝑠𝑎𝑎𝑎𝑎𝑡𝑡𝑎𝑎𝑠𝑠𝑠𝑠ℎ𝑠𝑠𝑡𝑡𝑑𝑑 (ℎ, 𝑥𝑥, 𝑥𝑥; 𝜔𝜔) = 𝐹𝐹 𝑑𝑑𝑑𝑑𝑡𝑡𝑠𝑠𝑑𝑑𝑠𝑠𝑠𝑠,𝑠𝑠𝑎𝑎𝑎𝑎𝑡𝑡𝑎𝑎𝑠𝑠𝑠𝑠ℎ𝑠𝑠𝑡𝑡𝑑𝑑 𝑥𝑥̇= 𝐹𝐹 𝑑𝑑𝑑𝑑𝑡𝑡𝑠𝑠𝑑𝑑𝑠𝑠𝑠𝑠,𝑠𝑠𝑎𝑎𝑎𝑎𝑡𝑡𝑎𝑎𝑠𝑠𝑠𝑠ℎ𝑠𝑠𝑡𝑡𝑑𝑑 𝑥𝑥̈⋅ 𝑡𝑡 (5)
As demonstrated from the results listed in Table 1, the resulting damping coefficients are independent of the perturbation displacement value. The damping coefficients for the taperedland bearing at a constant minimum film thickness and acceleration have been calculated for three different perturbations. The percentage difference from the 0.1 μm case is calculated and is within the computational error margin.

Single Degree of Freedom (SDOF) model

The transient thrust bearing response has been predicted through a single degree of freedom lumped-mass model, developed with MATLAB/Simulink. In this specific case, only axial motions have been considered. The previously derived bearing-specific coefficient maps generated through the afore mentioned methodology, and are inputted into the SDOF model to calculate the restoring and damping forces at each time-step. The bearing setup is always under a state of compression and the rest position distance hrest has been considered. Initially, the rotor is placed at the rest position, thus, there is no external or reaction force. Any load displaces the rotor from this position to a distance x, depending on the magnitude, compressing the rotor, as shown in Fig. 1. Therefore: The model's motion is described by the following equation:

ℎ(𝑡𝑡) = ℎ 𝑡𝑡𝑟𝑟𝑠𝑠𝑠𝑠 -𝑥𝑥(𝑡𝑡) (6)
𝑚𝑚 ⋅ 𝑥𝑥̈+ 𝑐𝑐 ⋅ 𝑥𝑥̇+ 𝑘𝑘 ⋅ ℎ = 𝐹𝐹(𝑡𝑡) ⇒ 𝑥𝑥̈(𝑡𝑡) = 1 𝑚𝑚 [𝐹𝐹(𝑡𝑡) -𝑐𝑐 ⋅ 𝑥𝑥̇-𝑘𝑘 ⋅ ℎ] (7)
In Eq. ( 7), m expresses the rotor's mass, while k and c are the stiffness and damping coefficients, respectively, being calculated by the CFD-based THD computational approach.

The SDOF equation has been used to calculate the rotor axial acceleration at each time-step. Zero initial conditions are set for the displacement and the velocity. Time integration calculates the velocity and the displacement at each time-step. Subsequently, stiffness and damping coefficients are extracted from the maps created, and the restoring and damping forces are defined. While the stiffness coefficient only requires the minimum film thickness to be defined, the damping coefficient is also a function of the velocity and acceleration as stated above. To deal with the issue, all the interpolations for the damping coefficient relationship with acceleration have been done considering that acceleration is equal to the RMS value of the load's acceleration 𝑎𝑎 𝑙𝑙𝑎𝑎𝑠𝑠𝑑𝑑 (𝑡𝑡) = 𝐹𝐹(𝑠𝑠)

𝑑𝑑

. Depending on the direction of motion, either the distancing or the approaching coefficient is used (𝑥𝑥̇> 0 → 𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑐𝑐ℎ𝑖𝑖𝑖𝑖𝑖𝑖, 𝑥𝑥̇< 0 → 𝑑𝑑𝑖𝑖𝑑𝑑𝑡𝑡𝑎𝑎𝑖𝑖𝑐𝑐𝑖𝑖𝑖𝑖𝑖𝑖). The angular velocity ω has been considered to remain constant for each simulation. The above described methodology is represented graphically in Fig. 2. If a constant force is applied, the resulting minimum film thickness from the SDOF model, should be equal to that of the steady-state CFD model. Thus, the stiffness coefficient does not include any dynamic parts.

Thermohydrodynamic (THD) Model setup

The steady-state THD model has been based on that developed by Charitopoulos et al. in [START_REF] Fouflias | Performance comparison between textured, pocket, and tapered-land sectorpad thrust bearings using computational fluid dynamics thermohydrodynamic analysis[END_REF][START_REF] Charitopoulos | Design optimization of an automotive turbocharger thrust bearing using a CFD-based THD computational approach[END_REF] and solves the Navier-Stokes equations coupled with the Energy equation to account for thermal effects by simulating a conjugate heat transfer problem between lubricant, pad and rotor. Cavitation, viscous dissipation and the temperature dependency of the lubricant's viscosity have also been considered.

The equations, solved with the CFD code ANSYS CFX involving an incompressible flow, neglecting gravitational and external body forces, are:

Mass conservation equation:

𝛻𝛻 • 𝑢𝑢 �⃗ = 0 (8)
Momentum equations:

𝜌𝜌 𝐷𝐷𝑢𝑢 �⃗ 𝐷𝐷𝑡𝑡 = -𝛻𝛻𝑎𝑎 + 𝜇𝜇𝛻𝛻 2 𝑢𝑢 �⃗ (9)
Fluid domain energy equation:

𝜌𝜌 • 𝑐𝑐 𝑎𝑎𝑝𝑝 • 𝑢𝑢 �⃗ • 𝛻𝛻𝛻𝛻 = 𝛻𝛻 • �𝜆𝜆 𝑝𝑝 • 𝛻𝛻𝛻𝛻� -𝜏𝜏: 𝛻𝛻 • 𝑢𝑢 �⃗ (10)
Solid domain energy equation:

𝛻𝛻(𝜆𝜆 𝑠𝑠 • 𝛻𝛻𝛻𝛻) = 0 (11)
where, 𝑢𝑢 �⃗ is the fluid velocity vector (m/s), p is the fluid pressure (Pa), T is the fluid/solid temperature (K), τ is the viscous stress tensor (Pa), ρ is the oil density (kg/m 3 ), μ is the oil dynamic viscosity (kg/(m•s)), cpf is the oil specific heat capacity (J/(kg•K)), λf is the oil thermal conductivity (W/(m•K)), and λs corresponds to the thermal conductivity of the pad and the rotor (W/(m•K)).

Bearing Geometrical Configurations

The present study has considered two thrust bearing designs, both having an eight-sector pad configuration, whose main geometric characteristics are given in Table 2.

The first design is a tapered-land pad as shown in Fig. 3a. Each pad has a linear wedge where pressure buildup is developed and is followed by a land parallel to the rotor plane. Each pair of consecutive pads is separated by a groove. The bearing selected in this study has been derived from a parametric analysis (considering nominal minimum film thickness of 20 μm) to determine the optimal taper inclination and land extent, for the specific dimension, by maximising the load to friction torque ratio.

The second design is a parallel surface-textured bearing with the same principal dimensions (Fig. 3b). The density of the textured area has been set to 80%, the configuration selected comprise of twenty-four (24) orthogonal dimples, four (4) in the radial direction and six [START_REF] Papadopoulos | Characterisation of stiffness and damping in textured sector pad micro thrust bearings using computational fluid dynamics[END_REF] in the circumferential. The textured parameters, namely textured length, width, and depth have been selected by a genetic algorithm optimisation [START_REF] Papadopoulos | Evolutionary Optimization of Micro-Thrust Bearings With Periodic Partial Trapezoidal Surface Texturing[END_REF]. The selected parameter values are presented in Table 2 and are the resulting optimal geometry generated by Chalkiopoulos et al. [START_REF] Chalkiopoulos | Effects of thermal and mechanical deformations on textured thrust bearings optimally designed by a THD calculation method[END_REF].

Mesh Characteristics

Regarding the tapered-land bearing, a structured mesh, whose fluid film and pad regions are depicted in Fig. 4a, totaling 730,000 elements has been generated, consisting of 121 elements across the circumferential direction, 65 across the radial direction and 90 elements across the height, 10 out of which inside the film.

As for the textured pad bearing, a mesh study has been conducted to select the mesh parameters. Regarding the pad domain, which appears together with the fluid domain in Fig. 4b, a dense hexahedral mesh has been utilised, composed of 910,000 elements. The fluid domain is discretised with hexahedral elements, utilising 15 layers of elements in the film thickness direction, while in the longitudinal and transverse directions, 150 and 60 elements have been used, respectively. The total number of elements in the fluid and rotor domains is 340,000 and 350,000, respectively.

During the transient simulation, as the rotor is forced to approach or move away from the pad, the fluid film gets thinner and thicker, respectively. As a result, as the element number in the fluid film thickness remains constant, the elements' dimension in the thickness direction is altered depending on the rotor's movement.

Boundary Conditions

Both bearing designs have been modeled with the same boundary conditions. Regarding the fluid domain, the internal boundary (named inlet) has been considered a pressure inlet. The oil is fed at a constant pressure, however, depending on the pressure distribution in the fluid domain, oil may be pushed out in locations where the pressure is higher than the feeding value.

On the other side, the fluid is allowed only to exit through the outlet. The outer side of the grooves is covered from the pad's body; thus, a wall is placed which does not allow oil outflow. Since only one pad has been considered, the sides have been modeled with a rotational periodicity interface to account for the adjacent pads. Finally, the two walls that are on the top and bottom side are the rotor and pad interfaces, respectively. In both cases, thermal effects have been taken into account by considering that the heat flux on either side of the interface is the same. In addition, while the pad is fixed in space, the rotor rotates, dragging the fluid's upper layer due to the no-slip condition. Consequently, the fluid has the same velocity field as the rotor. The hydrodynamic and thermal boundary conditions are listed in Table 3 and Table 4, respectively, while the boundaries can be seen in Fig. 5.

Fluid domains Oil

Inlet

Sides

Rotational periodicity Rotational periodicity

Oil interface

Heat flux and temperature continuity Heat flux and temperature continuity

Material thermophysical properties

The lubricant considered in the present study is an ISO-VG 46 oil, whose thermophysical values are listed in Table 5. All properties, apart from viscosity, have been considered temperature independent. The viscosity temperature dependence has been approximated using the ASTM chart equation, an empirical equation based on McCoull and Walther's equation:

�(𝜈𝜈 + 𝑏𝑏) � = 𝑎𝑎 -𝑐𝑐(𝛻𝛻) (12)
where, ν is the kinematic viscosity (cSt), T is the temperature (K), a = 9.02865, b = 0.6 and c = 3.52681.

Cavitation has been taken into consideration using the Rayleigh-Plesset model. Vapour pressure has been considered constant and equal to 50 kPa. Both rotor and stator have been considered to be made of steel, whose thermophysical properties are presented in Table 6.

Single degree of freedom parameters

Regarding the single degree of freedom model, it has been considered that the rest position is at a distance of hrest = 50 μm from the pad and that the moving mass m weighs 4 kg. The model is set at an initial fixed value of h for a total time (tsteady) of 0.01 s followed by a transient loading condition (ttransient) that lasts 0.03 s. Two initial conditions and four transient loads have been examined. Regarding the initial conditions, the first considers that the bearing is at the nominal operating condition minimum film thickness of hmin = 20 μm, while the other considers that both designs are set at the same loading force F = 350.41 N. As for the transient part, a step impact load (Fig. 6a), a simple sinus (Fig. 6b) and a sum of two sinusoidal functions (Fig. 6c) have been coupled with the first initial condition to evaluate the methodology proposed. Finally, a complex load (Fig. 6d), defined in Eq. ( 13) and Table 7, is coupled to both initial conditions to compare the dynamic performance of the two thrust bearing designs.

𝑑𝑑𝐹𝐹(𝑡𝑡) = 𝐴𝐴 1 • 𝑑𝑑𝑖𝑖𝑖𝑖 � 𝜔𝜔 2 𝑡𝑡 + 𝜑𝜑 1 � + 𝐴𝐴 2 • 𝑑𝑑𝑖𝑖𝑖𝑖(𝜔𝜔𝑡𝑡 + 𝜑𝜑 2 ) + 𝐴𝐴 3 • 𝑑𝑑𝑖𝑖𝑖𝑖(2𝜔𝜔𝑡𝑡 + 𝜑𝜑 3 ) + 𝐴𝐴 4 • 𝑑𝑑𝑖𝑖𝑖𝑖(4𝜔𝜔𝑡𝑡 + 𝜑𝜑 4 ) + 𝐴𝐴 5 • 𝑑𝑑𝑖𝑖𝑖𝑖(8𝜔𝜔𝑡𝑡 + 𝜑𝜑 5 ) + 𝐴𝐴 6 • 𝑑𝑑𝑖𝑖𝑖𝑖(16𝜔𝜔𝑡𝑡 + 𝜑𝜑 6 ) (13)

Results

In the following figures, the dynamic coefficients are displayed as a function of minimum film thickness hmin, considering the nominal rotational velocity of n = 6,000 RPM. All the results that follow have been calculated for one pad. The tapered-land geometry generates higher load values for a specific minimum film thickness, in comparison to a textured geometry of similar principal dimensions at the same operating conditions. As a consequence, the stiffness coefficient, calculated by Eq. ( 1) is higher for the tapered-land geometry bearing, as shown in Fig. 7. Regarding the damping coefficient, the textured bearing has a significantly higher damping effect and is less linear than the tapered-land's coefficient, both in the approaching and distancing directions, as it can be seen in Fig. 8 and9, respectively.

(a) (b) The dynamic coefficients that resulted from the CFD analysis have been used in the single degree of freedom model for the two designs optimised at hmin = 20 μm.

To evaluate the methodology, the transient loading scenarios displayed in Fig. 6a, 6b and 6c have been coupled to the initial condition of 20 μm initial film thickness for the tapered-land design bearing and applied to the SDOF model. The resulting minimum film thickness time function has then been given as input to a transient CFD model of the bearing and the resulting force has been compared to the initial, as shown in Fig. 10.

As it can be seen from the comparison, while there are some differences, the initial load is predicted with good accuracy. In the step load scenario, there is only a 1% steady load difference between the applied load and the predicted through the method. In the simple sinus and sum of two sinus function scenarios, the major differences are noted at the peaks and valleys. In both cases, the maximum error is only a 2% of the initial load. Therefore, the methodology proposed can be trusted to give accurate results for the dynamic behaviour prediction. As shown in Fig. 15, the textured design, due to the higher damping coefficient, is less affected by the rapid load changes in both scenarios. However, in the case of same initial minimum film thickness, the maximum displacement difference is higher since the stiffness coefficient is lower. On the other hand, when subjected to the same load, the tapered-land design operates at a higher minimum film thickness, resulting in a reduced stiffness coefficient. Consequently, the tapered-land's maximum displacement difference is higher in that case.

(a) (b) The results indicate that while both designs follow the load's trend, the tapered-land bearing's results are more accurate compared to the textured one. More specifically, in the first initial condition, the tapered-land's results are almost identical to the imposed load, while in the second initial condition there are some differences, which may be caused by the denser data set collected near the optimum minimum film thickness. Regarding the textured pad bearing, the results are not as accurate. However, the figures represent the results without the steady part of the load. For instance, the peak displayed in Fig. 16a between 0.02 s and 0.025 s is approximately only 6% higher compared to the total load.

Conclusions

In the present work, a procedure has been proposed to calculate the transient thrust bearing response using CFD analysis, comprising the following steps:

• Calculation of the dynamic coefficients of the bearing using a combination of steady-state and transient CFD-based THD simulations, at different values of minimum film thickness, as well as rotor direction and acceleration levels. • Development of a SDOF dynamic model, capable of calculating the transient bearing response (rotor motion), for a given time-varying thrust load. • Feeding of the rotor motion, predicted by the SDOF model, back to the CFD bearing model, and calculations of bearing performance using a time-dependent mesh deformation computational approach.

The proposed methodology has been applied to both a tapered-land and a parallel textured fluid film thrust bearings, permitting the prediction of bearing response under transient thrust loading. This includes a steady-state simulation, in which the rotor has been set at a specific minimum film thickness, and a perturbation consisting in a constant acceleration motion to either approach or move away from the pad. Three acceleration values have been considered, according to the anticipated operational profile of the bearings. The results of the analysis yielded a map for each coefficient, at different values of rotational speed, thrust load and acceleration levels.

Subsequently, a single degree of freedom model has been developed to predict the response of the bearing under a transient loading condition. The single degree of freedom model results have been validated by setting the displacement derived, as a predefined motion of the rotor in a new transient CFD analysis. The load function that resulted showed good agreement with the initial loading force applied to the single degree of freedom model.

Comparing the dynamic performance of the two bearing designs, namely the tapered-land and the textured pad thrust bearings, the following conclusions have been drawn:

(a) The tapered-land bearing is affected more significantly by fast load changes while the textured bearing displays a higher maximum film thickness difference when both bearings are set at the same initial minimum film thickness.

(b) If the same load is applied to both bearings, the textured pad bearing exhibits smaller minimum film thickness variation, because it operates at smaller values of minimum film thickness.

It is noted that operating minimum film thickness values of thrust bearings become smaller over time, for improved efficiency, therefore the dynamic behaviour analysis is getting more important to avoid unpredictable response and failures. Further, experimental procedures to evaluate transient behaviour are both expensive and time consuming, whereas direct transient simulations for a given thrust load are practically impossible, due to the constant remeshing requirement, for the calculation of the correct position (film thickness) at each time-step. Thus, the present work proposes an alternative feasible procedure to calculate the dynamic performance of thrust bearings, utilising a SDOF bearing model, with dynamic coefficients appropriately calculated with state-of-the-art CFD approaches.

Figure 1 :

 1 Figure 1: Single degree of freedom model representation.

Figure 2 :

 2 Figure 2: Bearing's response calculation procedure.

Figure 3 :

 3 Figure 3: Geometrical configurations: (a) Tapered land bearing, (b) Textured pad bearing.

Figure 4 :

 4 Figure 4: Mesh details: (a) tapered-land, (b) textured

Figure 5 :

 5 Figure 5: Boundary condition faces.

Figure 6 :

 6 Figure 6: Transient load part scenarios: (a) Impact step, (b) Sinusoidal function, (c) Sum of two sinusoidal functions, (d) Complex load.

Figure 7 :

 7 Figure 7: Stiffness coefficient k to minimum film thickness hmin, 6,000 rpm.

Figure 8 :

 8 Figure 8: Distancing damping coefficient c to minimum film thickness hmin, 6,000 rpm, (a) tapered-land bearing, (b) textured bearing.

Figure 9 :

 9 Figure 9: Approaching damping coefficient c to minimum film thickness hmin, 6,000 rpm, (a) tapered-land bearing, (b) textured bearing.

Figure 10 :

 10 Figure 10: SDOF-CFD model evaluation.

Figure 11 :Figure 12 :

 1112 Figure 11: Tapered-land response -Initial condition Hmin = 20 μm, (a) Minimum film thickness and stiffness coefficient, (b) Velocity and damping coefficient.

Figure 13 :

 13 Figure 13: Tapered-land response -Initial condition F = 350.41μm, (a) Minimum film thickness and stiffness coefficient, (b) Velocity and damping coefficient.

Figure 14 :

 14 Figure 14: Textured response -Initial condition F = 350.41N, (a) Minimum film thickness and stiffness coefficient, (b) Velocity and damping coefficient.

Figure 15 :Figure 16 :

 1516 Figure 15: Tapered-land -Textured response comparison, (a) Initial condition Hmin = 20 μm, (b) Initial condition F = 350.41 N.

Table 1 : Perturbation analysis results.

 1

	Perturbation Δx [μm] Damping coefficient [N•s/m]	Difference [%]
		Distancing	Approaching Distancing Approaching
	0.05	18,673	18,771	3.6	5.4
	0.1	19,370	19,849	-	-

Table 2 : Thrust bearing geometrical characteristics.

 2

	Parameter	Symbol	Unit	Value
	Geometry			Tapered-land	Textured
	Optimum Hmin		[μm]	20	20
	Number of pads			8	8
	Outer diameter	Douter	[mm]	90	90
	Inner diameter	Dinner	[mm]	50	50
	Rotor height	Hrotor	[mm]	10	10
	Pad height	Hpad	[mm]	6	6
	Groove height	Hgroove	[mm]	4	4
	Groove extent	φgroove	[°]	4.92	4.92
	Taper extent	φgroove	[°]	33.75	-
	Taper inclination	H1 = Hmax -Hmin	[μm]	65	-

Table 3 : Hydrodynamic boundary conditions details.

 3

Table 4 : Thermal boundary conditions details.

 4

	Solid domains	Rotor		Pad
	Inner	Adiabatic boundary		Heat transfer coefficient:	750 W/mK
				Ambient temperature:	40°C
		Heat transfer coefficient:	25 W/mK	Heat transfer coefficient: 25 W/mK
	Outer
		Ambient temperature:	25°C	Ambient temperature:	25°C
		Heat transfer coefficient:	25 W/mK	Heat transfer coefficient:	1,000 W/mK
	Top / Bottom
		Ambient temperature:	25°C	Ambient temperature:	40°C

Table 5 : Lubricating oil parameters.

 5

	Description	Symbol	Unit	Value
	Material			Vapour	Oil
	Thermodynamic state			Gas	Liquid
	Molar mass	M	kg/kmol	383.92	383.92
	Density	ρ	kg/m 3	1.185	870
	Specific heat capacity	Cp	J/kgK	1004.4	2100
	Thermal conductivity	k	W/mK	0.0261	0.13
	Dynamic viscosity	μ	N•s/m 2	2.10E-05	-
	Kinematic viscosity, 40°C	ν40°C	cSt	-	47.26
	Kinematic viscosity, 100°C ν100°C	cSt	-	7.17

Table 6 : Thermophysical properties of the pad and rotor materials.

 6

	Description	Symbol	Unit	Value
	Thermodynamic state		Solid
	Molar mass	M	kg/kmol	55.85
	Density	ρ	kg/m 3	7854
	Specific heat capacity	Cp	J/kgK	434
	Thermal conductivity	k	W/mK	60.5

Table 7 :

 7 Load parameters.

	Amplitude	Value	Phase	Value
	A1	30	φ1	5.353297
	A2	40	φ2	3.878171
	A3	20	φ3	0.085219
	A4	15	φ4	2.010225
	A5	15	φ5	1.152626
	A6	15	φ6	3.710464