H El Aakib
email: hind.elaakib@gmail.com

J F Pierson

M Chaik

C Samba Vall

H Ait Dads

A Narjis

A Outzourhit

Evolution of the structural, morphological, optical and electrical properties of reactively RF-sputtered cobalt oxide thin films with oxygen pressure

Keywords: Cobalt oxide thin films, reactive RF-sputtering, optical properties, structural properties, electrical properties. I

Cobalt oxide thin films were deposited on glass and silicon substrates using reactive cathodic radio-frequency sputtering from a pure cobalt target in oxygen/argon gas mixtures with different percentages of oxygen (10, 30, 40 and 50 %). The X-ray diffractograms showed that deposited thin films crystallized in a cubic Co 3 O 4 phase. The crystallite size was found to vary with the oxygen pressure, reaching a maximum using 40 % of O 2 . Scanning electron microscopy (SEM), AFM, energy dispersive X-ray spectroscopy and Raman scattering analyses confirmed the formation of Co 3 O 4 films with nano-sized particles. In addition, it is found that the optical band gap of the films increased with increasing oxygen pressure in the sputtering atmosphere. Finally, the electrical resistivity of the sputtered thin films was shown to increase with oxygen-percentage.

Introduction

There is an increasing interest in the preparation and the investigation of materials with novel structural features and properties. Among these, transition metal oxides crystallizing in the spinel-type structure offer various optical, electrical and magnetic properties stemming from the multi-valence of the metal cations. In addition, nanostructuration is a promising approach to achieve unique materials properties due to specific surface increase and the effects of confinement [START_REF] Kodama | Magnetic nanoparticles[END_REF]- [START_REF] Gleiter | Nanocrystalline materials[END_REF].

Among transition metal oxides, cobalt oxide has received a particular interest because of its wide range of properties. Indeed, cobalt oxide Co 3 O 4 is an antiferromagnetic material that crystallizes in the Fd m space group with a mixed valence for cobalt ions: divalent Co 2+ and trivalent Co 3+ . The oxygen anions O 2-are distributed in a compact cubic stack with centered faces whereas the Co 2+ cations occupy the tetrahedral sites and the Co 3+ cations occupy half of the octahedral sites [START_REF] Zhu | Linear and nonlinear optical properties of Co_3O_4 nanoparticle-doped polyvinyl-alcohol thin films[END_REF], [START_REF] Oleksenko | Adsorption-semiconductor hydrogen sensors based on nanosized tin dioxide with cobalt oxide additives[END_REF].

Owing to it unique properties, Co 3 O 4 thin films have been potentially used in various applications such as thin film supercapacitors [START_REF] Kandalkar | Preparation of cobalt oxide thin films and its use in supercapacitor application[END_REF], selective gas sensor [START_REF] Wöllenstein | Cobalt oxide based gas sensors on silicon substrate for operation at low temperatures[END_REF], selective absorber for photothermal conversion [START_REF] Patil | Preparation and characterization of spray pyrolysed cobalt oxide thin films[END_REF], [START_REF] Avila | Cobalt oxide films for solar selective surfaces, obtained by spray pyrolisis[END_REF], fuel cells [START_REF] Benel | Synthesis and characterization of nanoparticulate La0.6Sr0.4CoO3-δ cathodes for thin-film solid oxide fuel cells[END_REF], catalysis [START_REF] Guyon | Deposition of cobalt oxide thin films by plasma-enhanced chemical vapour deposition (PECVD) for catalytic applications[END_REF], batteries [START_REF] Do | Preparation and characterization of CoO used as anodic material of lithium battery[END_REF], efficient visible light photocatalysts [START_REF] Roy | Ligand-assisted soft-chemical synthesis of selfassembled different shaped mesoporous Co 3 O 4 : efficient visible light photocatalysts[END_REF] and magnetic devices [START_REF] Makhlouf | Magnetic properties of Co3O4 nanoparticles[END_REF].

Consequently, several studies were devoted to low cost and environment-friendly techniques for the synthesis of cobalt oxides thin films with improved properties. In this context, different coating techniques were used to prepare the Co 3 O 4 thin films, such as electrodeposition [START_REF] Zheng | Electrochromic properties of nanoporous Co3O4 thin films prepared by electrodeposition method[END_REF], spin-coating [START_REF] Jayatissa | Fabrication of nanocrystalline cobalt oxide via sol-gel coating[END_REF], dip coating [START_REF] Drasovean | Optical properties of cobalt oxide films by a dipping sol-gel process[END_REF], chemical vapor deposition [START_REF] Cheng | Electrical conductivity of Co3O4 films prepared by chemical vapour deposition[END_REF], chemical spray pyrolysis [START_REF] Manickam | Cobalt oxide thin films prepared by NSP technique: Impact of molar concentration on the structural, optical, morphological and electrical properties[END_REF] and RF magnetron sputtering [START_REF] Liao | Structural characterization and electrochemical properties of RF-sputtered nanocrystalline Co3O4 thin-film anode[END_REF].

The main disadvantages of the wet methods [START_REF] Li | Preparation and characterization of nano-TiO2 powder[END_REF] lie in their incompatibility with device processing techniques, and limited flexibility in substrate choice, the inability to control growth of the films, and their limited capabilities for large-scale production. Physical processes such as reactive RF sputtering technique overcome these limitations and offer many advantages such as process reliability, production of high quality and uniform thin films with the desired chemical composition at industrial scale [START_REF] Yudar | Surface and optical properties of indium tin oxide layer deposition by RF magnetron sputtering in argon atmosphere[END_REF] and has a very high level of reproducibility unlike to the others mentioned methods. In addition, reactive sputtering allows to easily control the stoichiometry of the oxides and to deposit multi-layers of different compositions that can be tailored for specific applications (i.e., selective absorbers for solar applications).

In this work, reactive RF-sputtering method was used to make a further study the effect of oxygen pressure on structural, morphological, electrical and optical properties of sputtered cobalt oxide thin films. To the best of our knowledge, X-ray reflectometry is used for the first time in this work to follow the evolution of the density of the films with oxygen percentage in the sputtering chamber

EXPERIMENTAL DETAIL

Deposition of cobalt oxide thin films

Cobalt oxide thin films were deposited at ambient temperature by reactive RF sputtering on glass and silicon substrates that were beforehand ultrasonically cleaned in acetone and distilled water during 15 min. Before the sputtering process, the vacuum pressure in the chamber was lower to 10 -6 mbar. A pure metallic cobalt target (12 cm diameter and 0.2 cm thick with a purity of 99.99%) was pre-sputtered for 15 min and used to deposit the Co 3 O 4 thin films. The substrate-target distance was fixed at 7 cm. The reactive sputtering was carried out in an Ar/O 2 gas mixture with different O 2 percentages (10, 30, 40 and 50 %). The total pressure was 1 Pa while the RF-power was fixed at 200 W.

Characterization techniques

The crystal structure of obtained Co 3 O 4 thin films was investigated by Rigaku Smart lab SE X-ray diffractometer in the 2θ range of 20-80° with a CuKα radiation in the Bragg Brentano configuration ( = 1.5418 Å). Raman spectroscopy (Confotec -Model MR520 -3D Scanning Laser Raman Confocal Microscope) technique was used to identify all the obtained phases.

The surface topography of sputtered thin films was studied by an NTEGRA atomic force microscope (NT-MDT). The grain sizes were determined by the nova software that allows the visualization of data from atomic force microscopy. The elemental composition of films was determinate by using the Energy Dispersive X-Ray Spectroscopy (EDS -SSD Bruker). The optical properties were measured in the 500-2500 nm wavelength-range using an UV-Vis-NIR spectrometer (Varian Cary 5000). A profilometry technique was used to measure the thickness of the deposited films. The electrical resistivity of the obtained films was measured by the 4-point method using a Keithly 6430 sub-femtoamp remote source meter. and the increase in the full width at half maximum (FWHM) value of the peak (311). All the structural parameters such as lattice parameter, d-spacing, crystallite size and strain of the sputtered films are calculated from the XRD analysis and are listed in Table 1. The lattice parameter, a, was calculated by using the following formula [START_REF] Louardi | Characterization of cobalt oxide thin films prepared by a facile spray pyrolysis technique using perfume atomizer[END_REF]:

Results and discussion

X-ray diffraction study

(1)
Where d hkl is the lattice spacing and h, k, l are the Miller indices. The obtained lattice parameter varies from 8.09 to 8.13 Å by varying the oxygen percentage. This is in agreement with the results obtained in the JCPDF files (ref code: 80-1533). The crystallite size, D, of all the sputtered films was calculated by using Debye-Scherrer formula [START_REF] Bayansal | Effect of cobalt doping on nanostructured CuO thin films[END_REF]:

 (2)
Where  is the wavelength of Cu-Kα radiation (1.5418 Å) and β is the full width at the halfmaximum (FWHM). The crystallite size increases with increasing oxygen pressure reaching the maximum at 40% O 2 , which indicates the improvement of crystallinity. The crystallite size starts to decrease when the oxygen content exceeds 40%. On the other hand, the strain, , of the sputtered Co 3 O 4 thin films is calculated from the following equation [START_REF] Saleem | Simple Preparation and Characterization of Nano-Crystalline Zinc Oxide Thin Films by Sol-Gel Method on Glass Substrate[END_REF]:

 (3)
It was observed that the strain in the obtained films decreases with increasing pressure of O 2 until 40% O 2 and starts to increase. This shows that the increase in crystallite size is correlated with the decrease in strain, which indicates a decrease in lattice imperfections and formation of high-quality films.

According to XRD analyses, the 40% oxygen pressure is therefore considered optimal for the growth of Co 3 O 4 films with good crystalline quality under these conditions (200W) by increasing the deposition energy and consequently gives a large crystallite size.

3.2.X-ray reflectometry

The X-ray reflectance spectra of the all sputtered films are plotted in Fig. 2. These curves show a single stage of total reflection, which supports the homogeneity of the obtained films. The critical angle for total internal reflection (α c ) corresponds to the inflection point of the reflectance curve. This angle is related to the electronic density of the material, and then its density, by the following equation:

 (4)
Where r e is the classical electron radius = 2.8210 -15 m,  is the wavelength of Cu-Kα radiation (1.5413Å),  e is the electronic density, which is given by a weighted sum of the contributions of the various elements of the compound. The electronic density is related to the density  of the material by:

(

) 5
where N is the number of elements in the formula of the compound, c i is the stoichiometric coefficient of the element, M i is its molar mass and f 1i is the dispersive component of the Xray atomic scattering factor for Cu-Kα of element [START_REF] Henke | X-Ray Interactions: Photoabsorption, Scattering, Transmission, and Reflection at E = 50-30,000 eV, Z = 1-92[END_REF].

The uncertainty in density calculation was determined by: 

by neglecting the uncertainty in the other factors, we have:

(7)
Typically, the uncertainty in the measured angle is less than 0.005° since the resolution of the instrument is 0.001°, yielding a relative uncertainty of less than 4% [START_REF] Fujii | Recent Developments in the X-ray Reflectivity Analysis[END_REF].

The critical angle corresponds to the position of the minimum of the derivative of the reflectivity (Fig. 3). The results clearly indicate that the critical angle and consequently the density of the films increase with increasing oxygen pressure (Fig. 4). The density of the films tends towards the bulk value of Co 3 O 4 (6.1 g/cm 3 ) as the oxygen pressure is increased. The difference between thin film and bulk density is usually attributed to the grain boundaries in polycrystalline thin films. In addition, the presence of an amorphous will also reduce the effective density of the film as amorphous materials have lower densities (30 %) lower than their crystalline counterparts due the presence of voids [START_REF] Shiojiri | Desities of Amorphous Thin Films[END_REF]. The presence of a nonstoichiometric surface layers resulting from the interaction of the films with the atmosphere (O 2 , H 2 O ...) is not ruled out as this layer may also affect the reflectometry measurements and thus the density [START_REF] Gil | Limitations of x-ray reflectometry in the presence of surface contamination[END_REF].

3.3.Raman data

The structural properties of Co 3 O 4 were also investigated by Raman spectroscopy as shown in Fig. 5. Indeed, Co 3 O 4 has 7 vibrational modes represented as follows:

Γ = A 1g + E g + 3F 2g + 5F 1u + 2A 2u + 2E u + 2F 2u (8)
The Raman active modes are A 1g , E g and 3F 2g and Raman inactive modes are 2A 2u , 2E u and 2F 2u . In regard to 5F 1u 4 are infrared active and the other is an acoustic mode [START_REF] Ho | Tunable band gaps of Co 3-x Cu x O 4 nanorods with various Cu doping concentrations[END_REF], [START_REF] Hadjiev | The Raman spectra of Co 3 O 4[END_REF]. The Raman spectra of all the prepared thin films on glass substrates show clearly five dominant peaks located at 185, 470, 506, 601 and 670 cm -1 , that are attributed respectively to the fives Raman active modes E g , 3F 2g and A 1g of Co 3 O 4 [START_REF] Diallo | South Africa; b Nanosciences African Network (NANOAFNET)[END_REF]. The E g and 3F 2g mode are attributed to vibration of the Co 3 O 4 tetrahedral [START_REF] Jiang | Synthesis of sphere-like Co3O4 nanocrystals via a simple polyol route[END_REF] and the A 1g mode is assigned to the symmetric stretching of the Co-O bond in the octahedral sites [START_REF] Jiang | Synthesis of sphere-like Co3O4 nanocrystals via a simple polyol route[END_REF]. 2.

It is worth noting however, that these peaks are slightly down-shifted to lower wavenumbers in comparison with those of the bulk Co 3 O 4 reported in the literature (ex. A 1g located at 690 cm -1 ). This effect is attributed to heating of the films induced by the Raman excitation laser. This downshift is reduced when the intensity of the laser is reduced or when a silicon substrate is used instead of the glass as result of the high thermal conductivity of silicon. Also, it is noted that the intensity of all the peaks increase with increasing oxygen percentage until 40% O 2 and start to decrease which is in good agreement with the XRD results. The intensity of peaks starts to decrease at 50% of oxygen pressure. This is may be to the presence of defects in Co 3 O 4 structure at 50% O 2 . This behavior may also be attributed to several parameters such as crystallite size, strain and non-stoichiometry or changes in bond lengths [START_REF] Popović | Correlation between Raman wavenumbers and P?O bond lengths in crystalline inorganic phosphates[END_REF]. Furthermore, as-expected the measured FWHM of the Raman peaks of the films is still larger than that observed in bulk crystalline Co 3 O 4 , which is on the order of 5 to 10 cm -1 depending on the peak [START_REF] Hadjiev | The Raman spectra of Co 3 O 4[END_REF]. This is generally attributed the phonon confinement as results of the small crystallite size. However, narrowing and a shift of the Raman peaks to a higher wavenumbers are also observed as the oxygen pressure increases to 40% suggesting an improved crystallinity of the film in accordance with the XRD measurements. These results suggest that the films are partially crystallized which is in good agreement with XRD measurements.

Optical properties

Transmittance

The optical transmittance spectra of the obtained films were recorded in the 400-2500 nm wavelength range and plotted in Fig. 6. It is clear that the transmission of Co 3 O 4 films increases with increasing oxygen pressure and reached a maximum of 71.5% in the near IR region. This is due in part to the decrease in the thickness of the films sputtered at higher oxygen pressure, which is evident from the profilometer measurements as shown in Table 3.

The reduction of thickness with the increased oxygen pressure reflects the less efficiency of sputtering by oxygen plasma when compared with argon plasma. It may also be possible that more oxygen in the plasma causes the oxidation of the Co target which weakens the sputtering rate. Therefore, the oxygen percentage increases and the transmission of deposited films increase.

Furthermore, two absorption edges can also be clearly distinguished in the transmittance spectra.

Absorption coefficient

Optical band gap

The optical band gap values of the films are determined by the Tauc equation [START_REF] Barreca | Composition and Microstructure of Cobalt Oxide Thin Films Obtained from a Novel Cobalt(II) Precursor by Chemical Vapor Deposition[END_REF]:

(αhυ) n = A(hυ -E g ) (10)
Where, A is constant, hυ is the photon energy, E g is the band gap energy and n = 1/2 for indirect transitions and n = 2 for the direct transitions. Fig. 8 shows plots of (αhυ) 2 as a function of photon energy for the Co 3 O 4 thin films, which reveal two direct band gaps associated with the two observed absorption fronts. As also shown in table 3, the two optical gaps increase slightly with oxygen partial pressure during the deposition. The obtained values are in perfect agreement and close to those of bulk Co 3 O 4 [START_REF] Barreca | Composition and Microstructure of Cobalt Oxide Thin Films Obtained from a Novel Cobalt(II) Precursor by Chemical Vapor Deposition[END_REF].

The existence of two band gaps is well explained in Fig. 9 by the fact that the valence band has a great O 2p-character, whereas the principal contribution to the conduction band is given by the Co (II) 3d orbitals [START_REF] Lakehal | Structural, optical and electrical properties of Ni-doped Co3O4 prepared via Sol-Gel technique[END_REF]. [START_REF] Barreca | Composition and Microstructure of Cobalt Oxide Thin Films Obtained from a Novel Cobalt(II) Precursor by Chemical Vapor Deposition[END_REF]. One can also notice the increase in band gap with increasing the oxygen pressure (Table 3).

Electrical properties

The electrical properties of the obtained Co 3 O 4 thin films at different oxygen percentages were investigated by using the 4-point prob method. The obtained values are given in Table 3 and plotted in Fig. 10. It is found that the resistivity values increases from 131 to 912 Ω cm as the oxygen percentage is increased. This increase of the resistivity can be explained by an over-stoichiometry in oxygen of the compound which causes an increase of a part of Co 3+ in octahedral site to maintain the electroneutrality of the compound [START_REF] Cheng | Electrical conductivity of Co3O4 films prepared by chemical vapour deposition[END_REF]. It can also due to the grain boundaries. XRD results suggest an increase in the crystallite size as the oxygen pressure increases. The results are confirmed by AFM measurements which indicate that the grain size increase also which suggest that grain boundary contribution prevails. It is worth noticing that a simple thermoelectric analysis (not shown here) revealed that the prepared

Morphological and elemental analysis

The surface morphology of the sputtered cobalt oxide thin films was also investigated by atomic force microscope. Figure 11 increasing the oxygen pressure from 10 to 50 %, whereas the amount of O increased with increasing the oxygen pressure. However, the Raman and XRD spectra unequivocally showed that the films are Co 3 O 4 . We conclude that the RF sputtering allows to obtain uniform film thickness over large areas, a high degree of Co 3 O 4 films adhesion and makes a good controlling of the quality of Co 3 O 4 thin films in terms of structural, optical and electrical properties for better uses as absorbers in photothermal and photovoltaic applications.

Conclusion

Fig. 1

 1 Fig.1shows the X-ray diffractograms of thin films deposited on glass at different oxygen

Fig. 1 .

 1 Fig. 1. X-ray diffractograms of Co 3 O 4 thin films deposited using different oxygen percentages.

Fig. 2 .Fig. 3 .

 23 Fig. 2. X-ray reflectometry diagram of Co 3 O 4 films deposited at various oxygen percentages.

Fig. 4 .

 4 Fig. 4. Evolution of the density as a function of oxygen pressure.

  Its width is related to the Co-O bond lengths distribution and distortion occurring in Co 3 O 4. The spectra confirm the Co 3 O 4 spinel structure occurrence for all as-deposited films, where Co 2+ and Co 3+ occupy the tetrahedral and octahedral sites of Co 3 O 4 , respectively. The corresponding Raman parameters for all the sputtered films and the comparison of their characteristics peaks are summarized in Table

Fig. 5 .

 5 Fig. 5. Raman specta of Co 3 O 4 film as function of oxygen percentage.

Fig. 6 .

 6 Fig. 6. Transmission spectra of the CO 3 O 4 nanoparticles at different oxygen percentages.

Fig. 7

 7 Fig. 7 shows the variation of absorption coefficient, α, with photon energy, hυ, for the Co 3 O 4

Fig. 7 .

 7 Fig. 7. Absorption spectra of Co 3 O 4 thin films at different oxygen pressure.

Fig. 8 .

 8 Fig. 8. Plots of (αhυ) 2 as function of hυ for Co 3 O 4 thin films at different oxygen percentages.

Fig. 9 .

 9 Fig. 9. Schematic representations of the band structure of Co 3 O 4 .

Co 3 O

 3 [START_REF] Zhu | Linear and nonlinear optical properties of Co_3O_4 nanoparticle-doped polyvinyl-alcohol thin films[END_REF] thin films show positive Seebeck coefficient, which is well expected since this material is a p-type semiconductor at room temperature.

Fig. 10 .

 10 Fig. 10. Resistivity of Co 3 O 4 thin films as a function of oxygen percentage.

Fig. 11 .

 11 Fig. 11. AFM images of Co 3 O 4 thin films grown at (a) 10%, (b) 30%, (c) 40% and (d) 50% O 2 .

Co 3 O

 3 [START_REF] Zhu | Linear and nonlinear optical properties of Co_3O_4 nanoparticle-doped polyvinyl-alcohol thin films[END_REF] thin films were successfully deposited by reactive RF-sputtering at different oxygen percentages. The XRD and Raman results showed a pure Co 3 O 4 phase crystallized in cubic centered face structure with different particles sizes and high surface area. The X-ray reflectometry showed that the obtained density of films tends towards the theoretical value Co 3 O 4 (6.11 g/cm 3 ) as the oxygen pressure is increased. The optical properties indicated that the energy gap value depends on the oxygen pressure. Besides, two band gaps 1.46-1.50 eV and 1.87-2.17 eV were observed corresponding to edges of the charge transfer band Co 3+  Co 2+ and that of oxide ligand to metal O 2- Co 2+ and O 2- Co 3+ charge transfer band. The electrical conductivity of the sputtered films decreases with the increasing oxygen pressure. These studies confirm that the oxygen pressure strongly affects the properties of cubic structured Co 3 O 4 thin films. And according to XRD analysis, the 40% oxygen pressure is therefore considered optimal for the growth of Co 3 O 4 films with good crystalline quality.

Table 1 .

 1 Structural parameters derived from XRD measurements of Co 3 O 4 thin film films.

	Oxygen pressure (%)	Orientation	2θ (°)	FWHM (°)	Crystallite Size D (nm)	Strain (ⅹ10 -3 )	Lattice parameter (Å)	d-spacing (Å)
	10	(311)	36.36	0.95	8.74	4.14	8.13	2.45
	30	(311)	36.36	0.83	10.01	3.62	8.10	2.44
	40	(311)	36.36	0.69	12.04	3.01	8.08	2.44
	50	(311)	36.36	0.76	10.93	3.31	8.09	2.44

Table 2 .

 2 Summary of Raman analyses

	Oxygen	FWHM
	pressure	of A 1g
	(%)	(°)

F 2g (cm -1 ) E g (cm -1 ) F 2g (cm -1 ) F 2g ' (cm -1 ) A 1g (cm -1 )

	10% O 2	39.92	182	467	509	606	668
	30% O 2	37.45	183	468	509	607	668
	40% O 2	34.35	185	470	510	608	675
	50% O 2	35.85	185	472	512	608	675

  The presence of Co (III) centers in Co 3 O 4 gives rise to an intermediate band located inside the energy gap. According to Fig. 8 and 9, the band gap values close to and below 1.5 eV are generally associated with a charge transfer Co 3+ d (t 2g )  Co 2 + d (t 2 ) while the values close to and below 2 eV, are associated with O 2-(2p 6 )  Co 2+ (3d 7 ) and O 2-(2p 6 )  Co 3+ (3d 6 ) charge transfer. The obtained values are in perfect agreement and close to those of bulk Co 3 O 4

Table 3 .

 3 Optical properties of sputtered thin films at different oxygen pressures.

	Oxygen pressures	Optical band gap	Resistivity	Thickness
	(%)	(eV)		(Ω.cm)	(m)
		Eg 1	Eg 2
	10	1.46	1.87	131.21	1.20
	30	1.47	1.94	197.58	1.15
	40	1.48	1.95	795.63	0.97
	50	1.50	2.17	911.85	0.80

  presents the atomic force micrographs of Co 3 O 4 thin films grown at different oxygen pressures. The formed films at 10% O 2 showed a fine grain structure with smooth surfaces. The average grain sizes are 120, 130, 150 and 140 nm for 10%, 30%, 40% and 50% O 2 respectively. The root mean square roughness increased from 1.08 to 1.28 nm with the increasing oxygen pressure reaching 50% O 2 . On the other hand, the average surface roughness of the sputtered thin films is 0.69, 1.00, 1.64 and 2.35 for 10%, 30%, 40% and 50% O 2 respectively. Such samples can be suitable as a good selective absorber.

Table 4

 4

shows the EDX analysis of the elemental composition of cobalt oxide thin films deposited on glass substrate. These measurements have been performed by taking into account only the Co and O elements. It was found that the amount of Co decreased with

Table 4 .

 4 Elemental composition of Co 3 O 4 thin films at various oxygen pressures.

	Oxygen pressure (%)	Element (weight %) Co O	Co/O
	10	62.1	37.9	1.64
	30	47.2	52.8	0.89
	40	46.6	53.4	0.87
	50	45.7	54.3	0.84

II.


C Droz

J.-P Lainé

M N Ichchou

G Inquiété

A reduced formulation for the free-wave propagation analysis in composite structures

Keywords: Wave propagation, composite waveguides, dispersion curves, reduction, wave finite elements

This paper presents an improved numerical strategy for the broadband analysis of wave propagation in composite or complex cross-sectional waveguides using the wave finite element method (WFE). Numerical analysis of such structures require highly discretized finite element models and leads to extensive computations. The proposed formulation relies on a projection of the cross-sectional transfer matrices on a reduced set of shape functions associated to propagating waves. Dispersion curves are then predicted only using a reduced number of eigenvectors. The performances and stability of this method are evaluated using the wavenumbers and wave shapes. Validations are provided for a sandwich composite beam and a cylindrical elasto-acoustic waveguide.

Introduction

Wave propagation in composite waveguides in a broadband frequency range is widely investigated in automotive and aerospace industry. A waveguide is a structure whose main dimensions exhibit a periodicity or homogeneity in such a way that the propagation of mechanical energy in the main direction (axis of a beam or plane of a plate) is privileged. Waveguide hypothesis can significantly reduce the size of the problem, since the behavior of a single sub-structure of the guide yields the response of the entire structure. Dynamical behavior of such structures is determined by evaluating the set of structural waves propagating through a cross-section [START_REF] Duhamel | Finite element analysis of the vibration of waveguides and periodic structures[END_REF][START_REF] Mead | The forced vibration of one-dimensional multicoupled periodic structures: An application to finite element analysis[END_REF]. One of the major interests of guided waves is their potential for travelling long distances at velocities governed by the dispersion phenomena. The knowledge of these dispersion properties for propagating waves is fundamental for an effective use in engineering, for example in the field of structural health monitoring (SHM). Numerical prediction of these different waves and their dispersion curves has been extensively studied in last decades. The semi-analytical finite element (SAFE) and wave and finite element (WFE) methods are, among others, very efficient tools for this purpose. In Email address: christophe.droz@gmail.com (C. Droz) the SAFE method, sinusoidal functions are employed to formulate the displacement field in the direction of propagation. Nevertheless, it is necessary to develop specific semi-analytical elements for each application, which can severely limit its interest for industrial purposes. In order to overcome these limitations, the WFE method combines periodic structure theory (PST) introduced in Mead [START_REF] Mead | A general theory of harmonic wave propagation in linear periodic systems with multiple coupling[END_REF] with a finite element method (FEM). Therefore conventional finite element software packages can be easily used to compute mass and stiffness matrices of the whole structure. The one-dimensional WFE method was successfully applied to a wide range of waveguides as beams-like structures [START_REF] Thompson | Wheel-rail noise generation part iii: rail vibration[END_REF][START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF][START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF], plates [START_REF] Ichchou | Multi-mode wave propagation in ribbed plates: Part i, wavenumber-space characteristics; part ii, predictions and comparisons[END_REF] and more complex geometries as thin-walled structures [START_REF] Houillon | Wave motion in thin walled structures[END_REF], tyres [START_REF] Waki | Free and forced vibrations of a tyre using a wave/finite element approach[END_REF], pipes [START_REF] Bareille | Distributed piezoelectric guided-t-wave generator, design and analysis[END_REF] and curved layered shells [START_REF] Chronopoulos | Predicting the broadband response of a layered cone-cynlindercone shell[END_REF]. As the application field of WFE method reaches structurally advanced composite structures, various numerical difficulties can appear, especially for one-dimensional formulation which involves larger cross-sections. Poor-conditioning of the transfer matrix can lead to numerical errors (see Zhong and Williams [START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF] for alternative formulations), aliasing effects and round-off errors can also appear if cross-section length is not sized carefully. However, for the determination of propagating waves in industrial waveguides involving a large number of degrees of freedom, major obstacle remains the large CPU time needed to solve the eigenproblem. Some numerical issues were investigated for example by Waki et al. [START_REF] Waki | Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides[END_REF] and a reduction strategy based on a contributing waves selection was proposed by Mencik [START_REF] Mencik | A model reduction strategy for computing the forced response of elastic waveguides using the wave finite element method[END_REF] to compute forced response of elastic waveguides [START_REF] Reno | On the forced response of waveguides using the wave and finite element method[END_REF]. Mencik et al. proposed a substructuring technique to compute the appropriate wave motion in multi-layered waveguides. Homogenisation techniques were also investigated [START_REF] Gopalakrishnan | Wave propagation, Diagnostics and control in anisotropic and inhomogeneous structures[END_REF] to apply the WFE to laminated composites.

Although the aforementioned reduction techniques are interesting to compute forced responses using a reduced number of wave modes, they do not reduce the numerical costs associated with the computation of the wave basis. Therefore, these techniques do not provide a general reduction strategy to compute the propagating waves in elaborated waveguides. This paper presents a method, based on classical formulation of the periodic structure theory, to calculate the dispersion curves of propagating waves for complex cross-section or composite waveguides involving an important number of degrees of freedom. This alternative formulation of WFE method relies on a frequency interpolation of the transfer matrix eigenvectors through a subset of eigensolutions. Thus, the propagating waves can be determined accurately solving a smaller eigenproblem, enabling the application of the WFE method to a wide range of sophisticated cross-sectional waveguide configurations.

The paper is organized as follows. In section 2, a brief overview of the classical WFE formulation is shown. Section 3 describes the proposed reduction strategy. The model order reduction is formulated for the spectral problem and the strategy of wave interpolation is described; the wave basis is defined next, using a reduced set of propagating waves computed at the cut-on frequencies, associated to the appearance of new propagating waves over the frequency band; a method is then proposed to improve the basis orthogonality and approximate eigenvectors between the cut-on frequencies. Numerical examples are brought in section 4. The first application concerns a three-layered sandwich beam; both the classical WFE formulation and an analytical low frequency solution described in [START_REF] Banerjee | Dynamic stiffness formulation and free vibration analysis of a three-layered sandwich beam[END_REF] and [START_REF] Mencik | A substructuring technique for finite element wave propagation in multi-layered systems[END_REF] are discussed, and the requirement for a refined FEM of cross-section is highlighted; the reduction strategy is then applied to a detailed FEM model. In the second example, the reduced WFE formulation is extended to an elasto-acoustic problem; dispersion curves are computed for a cavity filled with fluid and compared to the analytical solution.

Overview of the WFE

Free wave propagation in 1D-waveguides

In this section, a formulation of the WFE method is given for free wave propagation in a one-dimensional straight elastic and dissipative waveguide. The structure can be assimilated to N identical subsystems of length d connected along the main direction x. A unit cell of the waveguide is illustrated Fig. 1. Displacements and Figure 1: Illustration of a waveguide and the state vector of a unit cell [START_REF] Duhamel | Finite element analysis of the vibration of waveguides and periodic structures[END_REF].

forces are written as q and f, and subscripts 'L' and 'R' denote the left and right edges of a cell. Both edges have the same number n of degrees of freedom. Mesh compatibility is assumed between the N subsystems. The discrete dynamic equation of a cell at frequency ω is given by: (-

ω 2 M + jωC + K)q = f (1)
where M,C,K are the mass, damping and stiffness matrices, respectively. For periodic structures, condensation on the left and right cross-sections of the inner DOF's is required. Introducing the condensed dynamic stiffness operator D = -ω 2 M + jωC + K and reordering degrees of freedom, equation can be stated as follows:

D LL D LR D RL D RR q L q R = f L f R (2)
where D LL and D RR are symmetric and D t LR = D RL . Denoting λ = e -jκd the propagation constant describing wave propagation over the cell length d and κ associated wavenumber, considering force equilibrium λf L +f R = 0 in a cell and invoking Bloch's theorem [START_REF] Brillouin | Wave Propagation in Periodic Structures[END_REF], q R = λq L into Eq. (2) leads to the following quadratic spectral problem [START_REF] Duhamel | Finite element analysis of the vibration of waveguides and periodic structures[END_REF]:

(λ i D LR + 1 λ i D RL + D LL + D RR )Φ q i = 0 (3)
where ((Φ q ) i , λ i ) i=1,...,2n stands for the wave modes of the waveguide. The associated eigenvalue problem can be formulated by an appropriate state vector Φ q = [(Φ q ) t , (Φ f ) t ] t , leading to a symplectic transfer matrix T.

T Φ q Φ f = λ Φ q Φ f (4)
with

T = D -1 LR D LL D -1 LR D RL -D RR D -1 LR D LL -D RR D -1 LR (5)
where eigenvectors represents both nodal displacements and forces associated to a wave mode. The dynamical behavior of the global system can be expressed by expanding amplitudes of incident and reflected waves on a basis of eigenvectors. If the structure is undamped, solutions are divided into propagative waves, whose wavenumbers are real, and evanescent waves for which wavenumbers are imaginary. In dissipative case, complex wavenumbers are associated to decaying waves.

Computational issues for a complex cross-section

In practice, direct computation of the eigenproblem Eq. ( 4) can be prone to numerical errors when the meshed cross-section involves a large number of degrees of freedom. Indeed, the transfer matrix T requires to inverse matrix D -1 LR which can be poorly conditioned. To limit this issue, various formulations of the eigenproblem are available, for example:

-D RL -(D LL + D RR ) 0 -D RL q L λq L = λ 0 D LR -D RL 0 q L λq L (6)
However, numerical errors can become serious when eigenvalues are either very large (λ i ) or very small (1/λ i ), see Zhong and Williams [START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF] for a detailed discussion. When complex waveguides are considered, an insufficient discretization of the cross-section will produce significant errors, especially for eigensolutions associated to waves whose section shape have a short wavelength, whereas refined meshes exhibit numerous evanescent solutions, thus considerably increases computation time and worsen round-off errors due to the truncation of inertia terms, see Waki [START_REF] Waki | Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides[END_REF].

Yet, structures considered in this paper require a high degree of precision due to their geometry, their inner components or for high order wave shapes calculation. In these situations computation time grows exponentially with the number of nodes involved. Classical techniques based on modal basis reduction are not available, since a cross-section boundary conditions are arbitrary for a uniform waveguide or subjected to structural periodicity otherwise. Such an issue is addressed in the next section, introducing projection on a reduced set of shape functions.

Proposed reduction strategy

3.1. Spectral problem projection 3.1.1. Definition of a subspace for model order reduction The basic principle of Model Order Reduction (MOR) is to approximate Eq. ( 4) by creating a low order system over a specified range of frequencies. Spectral problem (3) can be rewritten:

S(λ, ω)u L (ω) = 0 ( 7
)
where S stands for:

S(λ, ω) = λD LR (ω)+(D LL (ω)+D RR (ω))+ 1 λ D RL (ω) (8)
Most MOR techniques can be interpreted as the formulation of two rectangular matrices B (r,n) and A (n,r) with r n, which are used for system (7) linear projection:

S(r,r) = B T S (n,n) A (9)
The quality of this approximation strongly depends on the subspaces spanned by these matrices. Although matrices A and B can be defined numerous ways, numerical stability commonly requires A to be real, equal to B and orthogonal [START_REF] Schilders | Model order reduction: theory[END_REF]:

A = B, B ∈ R n * r and B T B = I (10)
An effective way to obtain a reliable and low order reduced system is to make use of solution subsets from the initial system. This process is used, among others, for modal expansion or Krylov subspaces reduction [START_REF] Watkins | The Matrix Eigenvalue Problem: GR and Krylov Subspace Methods[END_REF].

Modal shapes -wave shapes

All reduced subspaces are not equivalent. In the field of structural analysis, modal analysis provides an efficient tool for numerous model order reduction strategies. Ritz method is based on the hypothesis that displacements field of a structure is spanned by its modal shapes {φ i } i=1,...,∞ . In a wave approach context, displacement of a bounded structure under harmonic excitation can be expressed as a sum of wave mode amplitudes [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF]. Each stationary mode being linked to a wave reflection on the boundaries, both approaches can theoretically lead to the same conclusions [START_REF] Langley | Some perspectives on wave-mode duality[END_REF]. Stationary modes are boundary dependent and modal approaches can lead to important system reduction if the structure has a low modal density, but suffer limitations in midfrequencies. Moreover, modal selection has numerous disadvantages, among other, a limited selection criteria: each mode is associated to a wave at a specific boundary dependent frequency. Consider the dispersion curves shown Fig. 2. A limited set of eigenvectors is retained at frequencies ω 0 , ω 1 and ω 2 , supposing an appropriate selection criteria can be found. They are denoted φ (0) i 1≤i≤r 0 , φ (1) i 1≤i≤r 1 and φ (2) i 1≤i≤r 2 . Considering that these eigenvectors represent the waves shapes through an elementary section, the wave shape φ(ω) shown Fig. 2 can be approximated as:

φ(ω) r 0 i=1 φ (0) i ν (0) i (ω) + r 1 i=1 φ (1) i ν (1) i (ω) + r 2 i=1 φ (2) i ν (2) i (ω) (11) where ν (0) i , ν (1)
i , ν (2) i are the reduced wave amplitudes. Eigenvectors φ i are associated to propagation constants (λ i , 1/λ i ) with the convention |λ i | ≤ 1 for positivegoing solutions and |λ i | ≥ 1 for negative-going solutions. Their norms are closer to unity for propagating waves and smaller for evanescent waves, in the dissipative case, |λ| < 1. These eigenvalues provide a simple criteria to sort and filter propagating waves (details are given in Section 3.2.2)). The following section describes the selection process and a method for computing an appropriate set of eigenvectors.

Definition of the wave basis

Selecting the appropriate waves

As noticed Section (3.1.2), accurate approximations require the exact solution to be a linear combination of the chosen eigenvectors. When a wave shape remains the same on a given frequency band, eigensolutions associated to this wave can be interpolated using an single, chosen eigenvector. This interpolation strategy is based on the hypothesis that major variations in the deformed shapes of propagating waves occur close to their cut-on frequencies. This is an important statement for eigenvectors interpolation. Therefore, computation of the exact solutions can be limited to a reduced set of chosen frequencies, where eigenvectors are prone to important variations. Singular deformed shape can be observed in several situations. First, when the wave enters the propagative domain the new mode of propagation is generally associated with a modification of its wave shape. In higher frequency, boundary effects can occur. Other phenomena as coupling effects, which are due to energetic exchanges between wave modes, can also contribute to waves variations. Waves shapes at cut-on frequencies are expected to be good candidates to enrich the reduced basis, supposing a effective determination of their frequencies can be done.

Identification of progressive waves appearance

Propagating waves have a positive wavenumber and low imaginary part, excepted for highly decaying waves. On the other hand, they are evanescent and have a negative wavenumber before they reach their cut-off frequencies. Considering wavenumbers are, among other, a continuous function of the frequency, the existence of a value ω A for which the wavenumber κ(ω A ) = 0 is assured. In terms of propagation constants, this means λ = 1. The frequencies ω A associated to eigenvalues λ = 1 can be determined by rewriting the quadratic eigenvalue problem Eq. ( 7):

S(1, ω A )Φ A = 0 ( 12
)
Wave cut-on frequencies and associated wave shapes Φ A are then solutions of the eigenproblem of size n/2:

(K * LL +K * RR +K * RL +K * LR )Φ A = ω 2 A (M LL +M RR +M RL +M LR )Φ A (13) where K * = (1 + jη)K is the complex stiffness consid- ering Rayleigh damping (α = 0, β = η/ω).
Let us denote n 0 the number of first-order waves, n 2 the number of higher-order waves in ∆ω and n A the number of cut-on frequencies determined using Eq. ( 13). As the first-order waves have the same cuton frequency (at ω = 0), the total number of cut-on frequencies is n A = 1 + n 2 . This means that the number of cut-on frequencies is lower than the total number: n 0 + n 2 of propagating waves in the frequency range.

Exact wave computation using classical WFEM

Denoting ω A(k) with 1 ≤ k ≤ n A the cut-on frequency of the kth wave, the exact eigensolutions can be computed for each new propagating wave from the system of n A independent quadratic eigenproblems:

                   S(λ, ω A(1) ) {Φ u L } = {0} ... S(λ, ω A(k) ) {Φ u L } = {0} ... S(λ, ω A(n A ) ) {Φ u L } = {0}                    (14)
Each line of system ( 14) can be solved using for example Eq. ( 6), and provides n eigensolutions (λ i , Φ i ):

Φ = Φ u L Φ u R (15)
Eigenvectors Φ u L i (ω A(k) ) stands for the deformed shape on the left side of a cell, associated with ith wave at kth wave cut-on frequency. Rectangular matrices A and B used for linear projection Eq. ( 9) can be build upon these eigenvectors.

Properties of wave functions

To summarize, the determination of a reduced set of eigenvectors is proposed to approximate full solutions over the whole frequency range, reducing the computation to one modal analysis and less than n A quadratic eigenproblems of size n/2, where n A is the number of cut-on frequencies in ∆ω. These vectors are independent of the cell length and a normalization can be performed to improve numerical conditioning.

Wave basis optimization

In this section, the reduced wave basis for linear projection is determined using selected eigenvectors. For a sake of clarity, unit-normalized deformed shape Φ u L i (ω A(k) ) will be written Φ i (k):

Φ i (k) = Φ u L i (ω A(k) ) , 1 ≤ i ≤ n (16)

Restriction to propagative waves

At this point, no selection criteria was applied to the n eigenvectors Φ i (k). If their wavenumbers are purely imaginary, they are evanescent. If a wavenumber is complex the wave is propagating and | (κ)| is the wave attenuation. In order to retain the least decaying waves and exclude highly dissipative waves, displacement can be expanded using only the positive-going propagative waves. This selection can be done for example retaining a number of eigenvectors for which (κ)

| (κ)|, say ñk n for the kth computed frequency. The solutions subspace will be spanned using ñ following vectors:

Φi (k) , 1 ≤ i ≤ ñk 1≤k≤n A and ñ = n A k=1 ñk (17)

A simple optimization process for wave functions

However, the concatenation of these wave shapes may involve redundant vectors, leading to a larger basis and numerical instabilities. This issue appears when two wave shapes ( Φi (k)) i and ( Φi (k + 1)) i are computed at close frequencies. In this case, vectors associated to the same wave can be correlated. Denoting . and . T the complex conjugate and transposition, respectively, we can define correlation between to complex vectors using modal assurance criteria (MAC):

Mac(φ, ψ) = (φ T ψ) (ψ T φ) (φ T φ) (ψ T ψ) (18)
Let us define the orthogonality relation, in the sense of MAC correlation, between two vectors φ and ψ as:

φ ⊥ ε ψ ⇔ Mac(φ, ψ) ≤ ε (19)
where ε ∈ [0 , 1] is a numerical criteria describing maximal tolerance for vector correlation in the reduced basis. We extend this definition to a vector's orthogonality to a set φ ⊥ ε {ψ} , as the vector φ being orthogonal to all the elements of {ψ}.

In order to retain a limited number of sufficient vectors an iterative filtering procedure can be done, defining:

Ψ(1) = Φi (1) , 1 ≤ i ≤ ñ1 (20)
then building the reduced basis from the recursive form:

Ψ(k + 1) = Ψ(k) , Ψ(k + 1) (21)
where Ψ(k + 1) contains the vectors of ( Φi (k + 1)) 1≤i≤ñ k+1 which are orthogonal to the elements of {Ψ(1), ..., Ψ(k)}.

The retained elements Ψ(k + 1) are defined as:

Ψ(k + 1) = χ ∈ ( Φi (k + 1)) 1≤i≤ñ k+1 \ χ ⊥ ε Ψ(k) (22)
The normalized and numbered vectors ( Ψ j (k)) j,k are denoted Γ i . It is also possible to evaluate wave correlation between two vectors Ψ i (k) and Ψ j (k + 1) at frequencies ω A(k) and ω A(k+1) using Zhong and Williams [START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF] method, in which case the following product is maximized:

(Ψ i (k)) T [J n ] (Ψ j (k + 1)) where [J n ] = 0 I n -I n 0 (23)

Building the reduced basis

Summarizing, the basis for wave expansion is defined using Eq. ( 13) to compute a reduced set of fre-quencies used to form system Eq. ( 14). Thus the reduced basis can be build from iterative procedure described Eqs. ( 20)- [START_REF] Langley | Some perspectives on wave-mode duality[END_REF]. Considering hypothesis on numerical stability Eq. ( 10), all the vectors of matrix A are real. This means that for an undamped waveguide vectors can be retained as they are defined earlier, whereas in dissipative cases it is preferable to display complex vectors Γ i as:

Γ * i 1 Γ * i 2 = (Γ i ) (Γ i ) ( 24
)
where and denote real and imaginary parts. However, the use of complex vectors for the reduced basis theoretically leads to the same conclusions, providing hermitian product for matrices projection. In most cases, wave dispersion is required on a frequency range starting at 0 Hz. Then, a set of undeformed motions, computed using system Eq. ( 25), is enough to describe first-order shape functions.

K LL K LR K RL K RR q L q R = 0 0 (25)
These real vectors, say n c , represent cell rigid motions of the waveguide under first-order waves. 3), matrix B (n,r) defined Eq. ( 9) can be expressed using notations Eq. [START_REF] Mencik | Wave finite elements in guided elastodynamics with internal fluid[END_REF] as:

B = Γ 1 , ..., Γ n c , Γ * n c +1 , ..., Γ * n c +r ( 26
)
where r is the sum of n c and the number of retained wave functions selected among the ñ vectors defined Eq. [START_REF] Banerjee | Dynamic stiffness formulation and free vibration analysis of a three-layered sandwich beam[END_REF]. Spectral problem Eq. ( 7) can now be rewritten:

λ DLR (ω) + ( DLL (ω) + DRR (ω)) + 1 λ DRL (ω) µ L = 0 (27)
where Di j stands for the projected dynamic stiffnesses:

Di j = B T D i j B (28)
The r × r eigenproblem can now be written:

-DRL -( DLL + DRR ) 0 -DRL µ L µ R = λ 0 DLR -DRL 0 µ L µ R (29)

Desired eigensolutions and error evaluation

The reduced eigenproblem can also be formulated:

S Φ = λ Φ (30)
where

S = 0 I D-1 LR DRL D-1 LR ( DLL + DRR ) (31)
and the (r × 1) reduced eigenvectors are denoted:

Φ = µ L µ R (32)
Reduced eigenvectors are written in spatial coordinates Φ = B Φ, as:

B = B 0 0 B (33)
The expanded (n × n) transfer matrix is defined as well: S = B SB T . The error on the reduced system eigensolutions ( λ, Φ), induced by the aforementioned basis reduction can be expressed as:

ε = (S -S) Φ Φ (34)
As S Φ = λ Φ, we obtain:

ε = (S -λI) Φ Φ (35)
Then, the error induced by the wave base truncation can be defined by injecting the reduced eigensolutions in the full eigenproblem. Thus no knowledge of the exact solutions (λ, Φ) is required.

Discussion

The size, r, of the reduced eigenproblem described above depends on the number of propagating waves in the frequency band ∆ω. The number of resolutions, n A , of the full eigenproblem is independent of the number of frequency samples N f s = ∆ω /δω + 1. As a result, the computational cost will be greatly reduced, solving n A N f s full eigenproblems of size (n × n) and N f s reduced problems of size (r × r) (n × n).

Algorithm for the reduced WFE method

The proposed reduction strategy is illustrated figure 3, each step is detailed as follows:

Build the mass and stiffness matrices of the waveguide cross-section and reorder as described Eq. ( 2).

Identify the cut-on frequencies over the frequency band using Eq. ( 13).

Solve the n A eigenproblems defined Eq. ( 14).

Select eigenvectors associated with propagating waves using section 3.3.

Build the wave basis Eq. ( 26) using the optimization process proposed Eqs. ( 20)- [START_REF] Langley | Some perspectives on wave-mode duality[END_REF].

Solve the N f s reduced eigenproblems of size (r ×r) given Eq. ( 29).

Derive eigensolutions using projection defined Eq. ( 33).

Numerical applications

Example 1: Sandwich beam

In this section, three solutions for free wave propagation in the sandwich beam illustrated Fig. 4, are compared. An analytical formulation depicted in Ref. [START_REF] Banerjee | Dynamic stiffness formulation and free vibration analysis of a three-layered sandwich beam[END_REF], and results obtained from WFE method in Ref. [START_REF] Mencik | A substructuring technique for finite element wave propagation in multi-layered systems[END_REF] are compared to a WFE formulation using a refined FEM for the cross-sectional discretization. The limitations of a coarse mesh are discussed and the requirement for a refined mesh for the WFE computation is highlighted. Then, the proposed reduction strategy is applied to the refined FEM model and results are discussed. The rectangular cross-section three-layered sandwich beam is composed of a rubber core surrounded by two steel layers. Detailed description of the layers is given Table 1. The width of the beam is 40 mm and the total thickness is 25 mm. A constant loss factor η = 0.01 is assumed for the numerical modeling while no damping is considered in the analytical formulation. and refined WFE formulations provide similar results for the longitudinal (1), torsional (2) and transverse flexural (3) waves. The flexural (4) wave slightly differs from the coarse model, while the complex thickness contractional ( 5) is clearly different between the two models. Note that the cut-on frequency for this last wave is 500 Hz, instead of 470 Hz for the coarse model, and a second-order longitudinal wave (6) cuts-on at 300 Hz in the refined model only.

After 500 Hz: the analytical solution (A1) and both WFE formulations appears to be consistent for the longitudinal wave [START_REF] Duhamel | Finite element analysis of the vibration of waveguides and periodic structures[END_REF]. Fig. 7 shows a zoom of the propagating wavenumbers previously depicted in Fig. 6. Three orders of longitudinal waves are highlighted, compared to the analytical solution and illustrated. Only the refined model differentiates the 3 longitudinal wave modes described Fig. 7. It can also be seen that several forms of longitudinal waves can coexist at the same frequencies. L2 and L3 waves exhibit the impedance mismatch between the skins and the soft core of the beam. The coarse model accuracy for these waves can be explained by possibility to approximate them using 4 linear elements in the beam thickness. The coarse model fails to compute the transverse flexural wave (A2) predicted by the analytical formulation. This wave mode appears to require a finer meshing for the inner rubber core. However, the analytical solution, based on a thickness homogenization does not entirely describe the waveguide behavior. Indeed, between 800 Hz and 1000 Hz three localized transverse flexural waves, presented Fig. 8, are propagating. They are associated to the transverse wave localization in the incorrect after 500 Hz. At 750 Hz, artificial veering occurs around k = 15 when the wave shapes turns more warped, leading to possible confusions in the wave modes matching.

The importance of a finely meshed model for the application of WFE method to a sandwich beam has been confirmed. The coarse model presented above fails to predict the wave propagation in this waveguide, especially over 500 Hz, and most of the wave shapes computed are erroneous. However, accurate dispersion curves and wave shapes can be obtained using the WFE method, under the condition that the model is highly discretized and the sampling frequency is small enough to match properly the wave modes.

Reduced WFE formulation

The proposed reduction strategy is applied to the refined model described above and compared to the classical WFE formulation. Real wavenumbers are presented Fig. 9. Vertical lines denote the cut-on frequencies, where the the full eigenproblems are solved. Note that all of the propagating waves starts at the predicted cuton frequencies. Other cuts-on are linked to important veerings associated to a propagating wave appearance, such as the thickness contractional wave at 500 Hz, or complex waves turning into propagating waves. There are 4 primary waves at 0 Hz and 8 cut-on frequencies, which is consistent with the 12 final branches identified on the dispersion curves Fig. 9. This result shows that the number of propagating waves was correctly predicted using Eq. ( 13). Very good correlation is ob- served between the classical formulation and the reduced model, built only using 10 eigensolutions: at the 2 frequency range boundaries and the 8 cut-on frequencies. The reduced wave basis is computed using a correlation criteria ε = 0.6, and contains 124 vectors instead of 1530 for the full system. Veering phenomena, where important variations of the wave shapes can occur, are accurately estimated using the reduced wave basis, see Fig. 9 at 650 Hz for the wavenumber k 13. The error on the wave shapes, defined Eq. ( 35) is presented Fig. 10 for the main propagating waves. Wave matching is provided in order to compare the errors associated with each propagating wave. A color scatter is also proposed to identify mostly propagating wave modes. It shows a good accuracy between full and reduced WFE, the maximal error being around 5×10 -4 . Cut-on frequencies are identified by vertical dashed lines. At 230 Hz, 300 Hz and 650 Hz, it is noticeable that each appearing wave shapes retained provide accurate solutions around the cut-on frequency. Significant reduction of CPU time is observed for the proposed WFE formulation, a detailed description is given table 2. This method provides a 97% reduction of CPU time in comparison with classical WFE formulation, and provides accurate dispersion curves and wave shapes. Note that most of the computation time for the reduced model is due to the 10 full computations, and similar accuracy was obtained by increasing the correlation criteria ε = 0.95 and using only 6 full computations, leading to a 148 vectors basis.

� 1�� 2�� 3�� 4�� 5�� 6�� 7�� 8�� 9�� 1��� 1� - 9 1� -8 1� -7 1�

Conclusions

In this example, the proposed reduced WFE formulation was applied to a finely discretized cross-section of a sandwich beam. First, the requirement for a refined FEM model was demonstrated by comparing the WFE solutions with an analytical solution and a coarsely meshed WFE model. Higher orders for propagating waves were identified. Then, the proposed reduced formulation was compared to a full computation and showed a very good agreement for both wavenumbers and wave shapes. Note that the proposed criteria for the choice of full computation frequencies can be excessive when numerous waves are propagating and lead to unnecessary rich wave basis lowering the numerical reduction efficiency.

Example 2: Elasto-acoustic waveguide

Description

As mentioned earlier, guided waves have numerous applications in structural dynamics. Elasto-acoustic waveguides are also a topic of interest, especially in SHM for crack and leakages detection. The theory of the propagation of vibro-acoustic waves in fluid-filled, cylindrical elastic shells has been studied in detail by Fuller and Fahy 1982 [23]. This example provides an extension of the reduced WFE formulation for a coupled fluid-structure system, where the inner fluid involves an important number of degrees of freedom to ensure decent discretization with FEM. The WFE formulation has been proposed by Mencik and Ichchou [START_REF] Mencik | Wave finite elements in guided elastodynamics with internal fluid[END_REF] to predict the forced response of a fluid-filled pipe shown Fig. 11a. The reduced WFE formulation can be applied as detailed Section 2 by using the dynamic stiffnesses described in [START_REF] Mencik | Wave finite elements in guided elastodynamics with internal fluid[END_REF], which exhibit quadratic dependence with respect to ω. A formulation is proposed in Appendix A to determine the cut-on frequencies by solving a linear eigenproblem in ω 2 instead of a quadratic problem in ω. In this example we consider a similar elastic pipe filled with an homogeneous, barotropic and compressible fluid whose viscosity is neglected. The cross-section shown Fig. 11b is defined as R e -R i = h with R e = 80 mm and h = 5 mm, the pipe has a Young modulus E s = 190 G pa, density ρ s = 7800 kg.m -3 and Poisson ratio ν = 0.29. The fluid has a sound velocity c a = 1484 m.s -1 and a density ρ a = 1000 kg.m -3 . The finite element model of a representative section of length d = 5mm involves three-dimensional linear ele- ments and n = n s +n a = 1706 degrees of freedom where superscripts s and a refer to solid and fluid domains.

Results

The dispersion curves are computed at frequencies between 0 Hz and 20 kHz and wavenumbers real parts are presented Fig. 12. A very good correlation is observed between the two WFE formulations over the whole frequency band. The 8 cut-on frequencies are identified and 10 full computations are performed. The analytical solution is given for torsional and longitudinal waves at low frequencies and is consistent with numerical solutions. The reduced wave basis is computed using a correlation criteria ε = 0.8 and contains 148 vectors. Significant reduction of CPU time is observed with this reduced formulation, a detailed description is given Table 3 ber and frequencies of the full computations retained. A reduced basis containing only low frequency eigenvectors is expected to produce significant errors as the frequeny induces changes of the waves shapes. An example is given Fig. 13, where the wavenumbers are computed from 0 Hz to 20 kHz, while only 4 cut-on frequencies are retained. The correlation criteria is ε = 0.8 and propagating waves are retained if | (k)| ≤ 0.5, leading to a (60 × 60) reduced system. A good accuracy is shown for the waves whose cut-on frequencies are retained. Wavenumbers are correct until 20 kHz while the last eigenvector is computed at 10,5 kHz. This result indicates that the reduced system solutions are consistent with the number of wave shapes retained for the reduced basis. However, no solutions are found for the higherorder waves. Noteworthy that the number of eigenvectors needed for the reduced basis can be overestimated if an important number of complex waves are retained. This issue can occur in highly dissipative structures.

Conclusions

In this example, the proposed reduced WFE formulation is applied to a cylindrical elastic shell filled with fluid. Dispersion curves are accurately determined with a reduction of 94% of the computation time. As expected, a truncated wave basis may not accurately predict the system solutions, the errors are consistent with the frequencies retained for the reduced wave base computations.

Concluding remarks

The main limitation of the WFE method for the broadband analysis of composite waveguides is due to their cross-sectional discretization. Coarsely meshed model produces inaccurate solutions and restrains the analysis to low-order propagating phenomena, especially within the field of composite structures where internal components can enable numerous localized wave modes. In this paper a reduction strategy was presented to study wave propagation in such waveguides. It provides a significant reduction of the computation time and enables broadband WFE analysis for a wide range of composite waveguides with complex cross-sectional geometries. The method was successfully applied to a refined model of a laminated composite beam. An error evaluation on the reduced eigenvectors was proposed and high-order localized wave modes were accurately predicted although a 97% reduction of the computation time. An elasto-acoustic pipe was considered. A adapted formulation was provided to identify the cut-on frequencies from a linear eigenproblem. The influence of a truncated basis was evaluated and showed good agreement with the proposed wave selection. Furthermore, the proposed method might be used to study highorder waves in composite waveguides such as wave localization and couplings. Further applications will be reported.

Figure 2 :

 2 Figure 2: Illustration of the dispersion curves used for the eigenvectors approximation.

3. 4 .

 4 Formulation of the reduced eigenproblem 3.4.1. Projection of dynamic stiffnesses Considering sections (3.1) and (3.

Figure 3 :

 3 Figure 3: Procedure for the reduced wave finite element method.

Figure 4 :

 4 Figure 4: Illustration: sandwich beam

  mesh used in Ref.[START_REF] Mencik | A substructuring technique for finite element wave propagation in multi-layered systems[END_REF] has 8 linear block elements having 3 degrees of freedom (DOF) per node and a total of 45 DOFs. The proposed refined mesh has 224 linear elements and 1530 DOFs, both cells have a length d = 2 × 10 -3 m in the direction of propagation, their meshes are shown Fig. 5.

Figure 5 :

 5 Figure 5: Cross-section discretization using a coarse mesh (left) and a refined mesh (right).

4. 1 . 2 .

 12 Refined FEM model for the WFE method Real wavenumbers obtained using the analytical, coarsely meshed WFE and finely meshed WFE for-mulations are presented Fig. 6 over the frequency range ∆ω = [0, 1000 Hz]. Below 450 Hz: coarse

Figure 6 :

 6 Figure 6: Dispersion curves for the global waveguide. Refined WFE: (red ---), Coarse WFE: (black --), Analytical solution: (black ----).

FrequencyFigure 7 :

 7 Figure 7: Deformed shapes associated with 3 orders of longitudinal waves and detailed wavenumbers.

Figure 8 :

 8 Figure 8: Deformed shapes associated with three transverse waves. T1 and T2: 3D view.

Figure 9 :

 9 Figure 9: Comparison between reduced (red ) and full (black -) wavenumbers for the sandwich waveguide. Vertical blue lines stand for predicted cut-on frequencies.

Figure 10 :

 10 Figure 10: Error ε for the propagating waves eigenvectors (|Im(k)| < 1).

Figure 11 :

 11 Figure 11: (a) Illustration of an elastic pipe filled with acoustic fluid. (b): Finite element model of a representative section of the elastoacoustic pipe.

Figure 12 :

 12 Figure 12: Comparison between reduced (red ) and full (black -) WFE formulations for the elasto-acoustic waveguide. Analytical low frequency solutions for axial quasi-longitudinal extensional and torsional waves (---).

Figure 13 :

 13 Figure 13: Wavenumbers for an elasto-acoustic waveguide. Reduced WFE formulation: (red ). Full WFE (gray -).

Table 1 :

 1 Material properties for the sandwich beam layers

	The coarse

Table 2 :

 2 Sandwich beam: CPU time for the full and reduced WFE formulations.

  . The reduction provides a 94% reduction of CPU time in comparison with classical WFE formulation, and gives accurate dispersion curves and wave shapes. The accuracy of this method relies on the num-

	Classical WFE	Number of Time per	CPU time
		iterations	iteration
	Resolution			Total:
	[0, 20000Hz] per 100Hz	×200	215s	11h 56m 40s
	Reduced WFE	Number of Time per	CPU time
		iterations	iteration
	Appearance ω A(k)	×1	28s	28s
	Full solution at ω A(k)	×10	215s	2150s
	Building reduced basis	×1	196s	196s
	Resolution
	[0, 20000 Hz] per 100 Hz	×200	0.08	16s
				Total: 39m 50s

Table 3 :

 3 Pipe filled with fluid: CPU time for the full and reduced WFE formulations.

Acknowledgement

The authors would like to gratefully acknowledge Airbus Helicopters for their financial support.

Appendix A. An extension for the fluid -structure interaction

Considering the (U,p) formulation for the elastoacoustic problem [START_REF] Mencik | Wave finite elements in guided elastodynamics with internal fluid[END_REF]:

where M s , K s , M a , K a stand for the structural and acoustic masses and stiffnesses, C is the fluid -structure coupling matrix, and F s , F a are the structural and acoustic forces vectors. Following the steps described section 3 and using same notation and D s,a = K s,aω 2 M s,a , the spectral problem can be written:

Then denoting Σ. = . LL + . LR + . RL + . RR and finding unknown cut-on frequencies ω A for wich λ = 1 leads to the eigenproblem:

Henceforth, the generalization of the reduced formulation is straightforward and leads to the eigenproblem:

with the dynamic stiffness operator [START_REF] Morand | Interactions fluids-structures[END_REF]:

) Solution µ defined Eq. 29 is written:

where u stands for the spatial displacement and Ψ is the acoustic velocity potential.