S Vassant

A Archambault

F Marquier

F Pardo

U Gennser

A Cavanna

J L Pelouard

J J Greffet

Epsilon Near Zero mode for active optoelectronic devices

Keywords: numbers: 73.20.Mf, 78.68.+m, 42.79.Ta, 85.60.-q, 78.20.Ci

The electromagnetic modes of a GaAs quantum well between two AlGaAs barriers are studied. At the longitudinal optical phonon frequency, the system supports a phonon polariton mode confined in the thickness of the quantum well that is called Epsilon Near Zero (ENZ) mode. This ENZ mode can be resonantly excited through a grating resulting in a very large absorption localized in the single quantum well. We show that the reflectivity can be modulated by applying a voltage. This paves the way to a new class of active optoelectronic devices working in the mid and far infrared at ambiant temperature.

A cornerstone of nanophotonics is the enhancement of electron-photon interaction by combining field confinement and resonances as in microcavities or nanoantennas. This enhancement can also be obtained by using surface waves, e.g. surface plasmon or surface phonon polaritons (SPhP) [1,2]. Examples of applications include the enhancement of Raman scattering [3], the enhancement of fluorescence [4,5], the enhancement of the radiative decay rate [6][7][8], localized heating [9]. More recently, a further step was done by introducing active plasmonic devices [10][11][12][13]. In the IR, SPhP allows designing resonant absorbers as demonstrated with polar semiconductors in the infrared [14,15] and THz [16] ranges. However, it is more difficult to use surface waves to enhance interaction with electrons in a quantum well (QW), as their spatial extension is generally more than one order of magnitude larger than a typical QW thickness and therefore the electron-photon spatial overlap is far from optimum. Here, we overcome this issue by the Epsilon Near Zero (ENZ) effect, which has been discussed in a series of recent papers [17][18][19][20]. Here, we introduce an ENZ SPhP mode (hereafter called ENZ mode), which allows enhancing the incident intensity by more than two orders of magnitude in a single QW. We first discuss the ENZ mode properties in the case of a single AlGaAs/GaAs/AlGaAs QW. We then demonstrate experimentally a resonant absorption due to the ENZ mode excitation. We finally report an electrical control of this absorption that allows modulating the reflectivity at ambient temperature.

Let us first introduce the ENZ mode. If one considers an interface between two media (1 and 2), the electric displacement ǫ zz E z normal to the interfaces is continuous so that E z2 /E z1 = ǫ zz1 /ǫ zz2 , where ǫ ij is the diagonal complex dielectric tensor and E the electric field. Thus if |ǫ zz2 | ≪ |ǫ zz1 |, the field in medium 2 is enhanced compared to the field in medium 1. We can define an enhancement factor

K EN Z = |E z2 /E z1 | 2 = |ǫ zz1 /ǫ zz2 | 2 .
We now consider the GaAs/Al 0.33 Ga 0.67 As system and show that it naturally exhibits ENZ enhancement. Our goal is to combine this effect with the resonant excitation of a SPhP to achieve a resonant absorption. In Fig. 1a, we plot the real part of the dielectric function of GaAs. The green (dark gray) shaded area corresponds to K EN Z ≥50, whereas the gray shaded area corresponds to the frequency range where GaAs and AlGaAs have the real part of their dielectric constants with opposite signs so that a SPhP can exist at a GaAs/AlGaAs interface. The ENZ and SPhP frequency range do not overlap for this geometry. However, if we consider a slab of GaAs embedded in AlGaAs (Fig. 1c), SPhP propagating at both interfaces can couple leading to a splitting of the SPhP dispersion relation in two branches [21]. As can be seen in Fig. 1b, the splitting increases when the slab thickness (d) decreases.

For small thicknesses (here 22nm), the frequency of the higher energy branch enters the ENZ frequency range but remains smaller than the longitudinal optical phonon frequency of GaAs (ω GaAs,LO ≈ 291.55 cm -1 [21]), so that the condition Re(ǫ GaAs (ω)) < 0 is fulfilled. We refer to this mode as ENZ mode. Its field distribution is depicted in Fig. 1c, for a GaAs thickness of 22 nm and a 7000 cm -1 wavevector. Two important features are seen: i) the field is confined in the QW; ii) the amplitude in GaAs is more than one order of magnitude larger than in Al 0.33 Ga 0.67 As. The system bears some resemblance with two surface plasmons propagating along a thin metallic film. In both cases, we have a metal-like medium embedded in a dielectric. The ENZ-mode has the same field symmetry as the long range surface plasmon (symmetric normal electric field) [22,23]. Despite these similarities, the ENZ mode and the long range surface plasmon have significant differences: i) a different field distribution: |E| is maximum in the film for the ENZ mode whereas it is minimum in the case of the long-range surface plasmon, ii) the ENZ mode dispersion relation is almost flat for a large range of wavevectors, with ω ≈ ω GaAs,LO .

We now show how the ENZ mode can be used for optoelectronic devices. Here, we limit the discussion to the 1a. When d=∞, the system is similar to a single interface geometry. As the thickness is decreased, the mode separates into two branches (symmetric and antisymmetric normal electric field distribution [22,23]). When d=22nm, the symmetric branch frequency enters the ENZ regime of GaAs. (c) left: Al0.33Ga0.67As/GaAs/Al0.33Ga0.67As slab geometry, right: spatial distribution of the z-electric field component of the ENZ mode at 34.3 µm for a wavevector of 7000 cm -1 and a GaAs thickness d of 22 nm.

electrical control of THz reflectivity. Other THz modulation schemes have already been demonstrated [24][25][26][27], but they encounter certain difficulties when working at higher frequency [28] and do not work in the optical phonon frequency range of their constitutive materials. Using optical phonons for optoelectronic devices can be surprising as they are usually considered to be detrimental. Indeed, the electron/optical-phonon interaction is a relaxation process much faster than the electron/photon interaction at work in most optoelectronics devices. Here, we use a different scheme based on the resonant absorption by the ENZ mode. Electrons are not used to absorb light but to modulate the ENZ effect. The principle is summarized here: incident light is coupled through a grating to the ENZ mode leading to a resonant absorption. The GaAs layer thickness is adjusted to 22 nm to form a QW where one intersubband transition (ISBT) energy is close to that of the ENZ mode. Thanks to the ISBT, the z-component of the QW dielectric tensor (see [21] for details) depends on the electron density in the well and can be tuned by a gate voltage. This gives an electrical control of the ENZ enhancement and thus of the absorption, as a lower ENZ enhancement leads to a smaller absorption.

The device is depicted in Fig. 2a. The incident THz field is coupled to the ENZ mode by a grating made of gold stripes which also serves as a Schottky gate for electrical control of the electron density in the single 22 nm thick GaAs QW. The substrate is a doped GaAs wafer, which behaves as a back-mirror at THz frequencies. The grating features and spacer thickness were chosen to optimize the absorption in the depleted QW through numerical simulations [21]. Details about the fabrication can be found in [21]. Reflectivity spectra[21], with a p-polarized incident wave (Fig. 2b), show that a 76% absorption is achieved at 34.3 µm (corresponding to the longitudinal optical phonon frequency 8.74 THz) when a negative gate voltage is applied. This bias depletes the QW, so that the ISBT contribution to the dielectric constant vanishes. No absorption is observed for s-polarization [21]. Numerical simulations, using bulk dielectric constants agree well with the reflectivity data (Fig. 2b). Moreover, they show that 45% of the incident energy is absorbed within the QW by the ENZ mode. The remaining 31% is absorbed in the rest of the structure. The calculated near field distribution at resonance confirms the strong |E z | 2 enhancement confined in the QW (Fig. 3). Fig. 2c shows that the resonant absorption does not depend on the angle, as expected from the almost flat dispersion relation of the ENZ mode. It is worth comparing this 45% absorption by the ENZ mode in a single QW with the intersubband absorption of a Quantum Well Infrared Photode- tector (QWIP) consisting of 50 QWs, which is typically less than 25% [29]. Furthermore, we did calculations on the same structure without a grating, and found a maximum absorption in the QW of 4%, only at the Brewster angle.

We now turn to the modulation of the reflectivity. Experimental reflectivity spectra at room temperature for different gate voltages are shown in Fig. 4a. It is seen that at the resonance (34.3 µm) the reflectivity decays from 0.38 to 0.25 when the applied voltage goes from +0.3 V to -0.3 V. The corresponding relative variation is as large as 54% as displayed in Fig. 4b. For larger voltages, the reflectivity change tends to saturate. The physical mechanism is based on controlling the ENZ effect with ISBT in the QW. As previously pointed out, we have an ISBT at the frequency of the ENZ mode (See [21] for more details). An increase of the electron density, through an increase of the gate voltage (-0.3 to +0.3 V) between the metallic grating and the doped substrate, increases the imaginary part of the intersubband contribution to the dielectric constant at the interface mode frequency. Hence ǫ zzGaAs is no longer close to zero, result- ing in a suppression of the ENZ field enhancement. By fitting the experimental reflectivity spectrum we obtain the values for the dielectric function for various electronic densities[21] (Fig. 4c) and calculate the corresponding K ENZ factor that strongly decreases as the electronic density in the QW increases (Fig. 4d). For practical applications, the time response of the device is an important figure of merit. We estimate [21] the maximum bandwidth, mainly limited by the capacitance of the device, to be on the order of 1 GHz. Note that the maximum electrical power flowing through the device is 425 nW at +0.3V so that thermal effects are negligible[21].

(c) (d) 0.3 V 0.2 V 0.1 V 0.0 V -0.1 V -0.2 V -0.
To summarize, we have introduced an ENZ mode that exhibits strong field enhancement. It can be used to engineer large absorption in a layer as thin as 22 nm (≈ λ i /1560). The strong interaction between the ENZ mode and the bidimensional electron gas allows controlling electrically the absorption. A relative change in the reflectivity up to 54% is experimentally obtained at 8.74 THz and room temperature with a single Al-GaAs/GaAs/AlGaAs QW. The ENZ effect is limited to a narrow range of wavelengths for a given material system, but is present in all polar crystals [START_REF]two different atoms per unit cell are polar and have at least one infrared active optical phonons[END_REF]. As a result this effect can be used at various wavelengths such as 11 µm in SiC, or 13.5 µm with GaN, where existing modulation scheme become challenging. The control of the ENZ effect appears to be a very promising mechanism to enhance light-matter interaction at the nanoscale. Due to Kirchhoffs law, enhancing and controlling the absorption means that emission can also be enhanced and controlled. Enhancing the ENZ mode electron interaction suggests that electrical transport measurements in the QW should allow the detection of incident radiation. Hence, using the ENZ mode opens new avenues for the design of modulators, sources and detectors.

  FIG.1:(Color online) (a) Real part of GaAs and Al0.33Ga0.67As bulk dielectric functions as a function of the frequency. The gray area corresponds to the frequency range where an SPhP at the GaAs/Al0.33Ga0.67As single interface can exist. The green (dark gray) shaded area correspond to the ENZ regime of GaAs, defined by KENZ = |ǫ zzAlGaAs /ǫzzGaAs| 2 ≥ 50. (b) Dispersion relation for the slab geometry and different GaAs thicknesses d. The grey and green (dark gray) shaded areas have the same meaning as in Fig.1a. When d=∞, the system is similar to a single interface geometry. As the thickness is decreased, the mode separates into two branches (symmetric and antisymmetric normal electric field distribution[22,23]). When d=22nm, the symmetric branch frequency enters the ENZ regime of GaAs. (c) left: Al0.33Ga0.67As/GaAs/Al0.33Ga0.67As slab geometry, right: spatial distribution of the z-electric field component of the ENZ mode at 34.3 µm for a wavevector of 7000 cm -1 and a GaAs thickness d of 22 nm.

FIG. 2 :

 2 FIG. 2: Excitation of ENZ mode. (a) Scheme of the device made of a gold lamellar grating with 4 µm period, 200 nm thickness, a 0.66 filling factor, a 22nm GaAs quantum well buried in Al0.33Ga0.67As barriers, a 1300 nm Al0.5Ga0.5As spacer[21] and a doped GaAs substrate. (b) Reflectivity spectrum for a p-polarized incident wave incident at 15 • , with an applied gate voltage of V=-0.3 V, data and theory. An absorption peak appears at 34.31 µm. The incident wave is p-polarized. Note that the second absorption peak near 36 µm is due to a surface plasmon propagating on the top of the metallic grating [21]. (c) Reflectivity for various angles of incidence for an applied voltage of V=-0.3 V. The absorption peak is almost independent on the angle of incidence. The reflectivity changes as a function of the angle of incidence above 34.6 µm are discussed in [21].

FIG. 3 :

 3 FIG. 3: Normalized near field distribution |Ez| 2 at the resonance (34.3 µm) for a p-polarized plane wave incident at 15 • on the structure. The white lines represent the geometry. Only one period of the grating is represented

2 FIG. 4 :

 24 FIG. 4: Electrical modulation of the reflectivity. (a) TM reflectivity for different applied voltages for an angle of incidence of 15 • . The inset shows the evolution of the reflectivity at 34.31 µm as a function of the gate voltage. (b) Corresponding relative reflectivity modulation between -0.3 and +0.3 V, data and theory. (c) Dielectric function of GaAs without (electronic density Ne = 0 cm -2 ) and with ISBT (Ne = 1.4× 10 11 cm -2 , γ isb = 330 fs) the ENZ enhancement factor KENZ as a function of the electronic density in the QW.


Daniele Barettin

A V Platonov

Alessandro Pecchia

V N Kats

G E Cirlin

A D Bouravleuv

Lucien Besombes

Henri Mariette

Matthias Auf Der Maur

Aldo Di Carlo

E Cirlin

D Bouravleuv

A F Ioffe

Model of a GaAs quantum dot embedded in a polymorph AlGaAs nanowire

Keywords: Nanowire, quantum dot, polymorph materials

I. INTRODUCTION

Since the discovery of self-assembled quantum dots (QD), the field of optoelectronics has moved interest on nanoobjects of different dimension, size and shape in the search of novel useful properties. The challenge of the QD fabrication technology is to obtain highly homogeneous arrays with little dispersion in size, composition and shape in order to maximize quantum efficiency of LEDs and lasers. This requirement has been solved only partially with self-organization [1] that allows to achieve a rather uniform distribution of the dot sizes with dispersion of 10% [2]. However, the technique has a number of essential ineradicable lacks. For instance, because of the presence of a wetting layer, it is impossible to reach high quantum efficiencies. Furthermore, this growth process offer little to no control over the dot size, an essential parameter in order to tune the emission wavelength.

Daniele. Barettin@uniroma2.it Thanks to the fast development and improvement of control in the growing conditions of nanostructures research has moved to new solutions in nano-engineering quantum confined states. The growth of free standing semiconductor nanowires (NW) and nanocolumns (NC) is one of the most promising solution for application in nanoelectronics [3],and nanophotonics [4] thanks to a better relaxation of the strain field that allows more freedom in combining III-V materials. The technology of manufacturing such NWs is now intensively developing [5]- [15]. These NWs can include one or several quantum dots (or quantum disks). Modern epitaxial methods allows to control the sizes and location of NW at the nanometer scale, that favorably distinguishes such QDs from planar and is very attractive for applications in optoelectronic devices.

However, GaAa/AlGaAs NWs also possess growth mechanisms that might have strong influences on their electronic properties. It is known [16] that the crystal structure of a <111> nanowire may spontaneously switch during the growth from zinc-blende (Zb) to wurzite (Wz) phases. The effect leads to a strong modification of the time resolved photoluminescence spectra in the range of the optical transitions within the barrier material [13]. Similarly, it is expected that optical transitions in the QDs are sensitive to the crystal phase of the dot itself and its surrounding media. For example the energy position of the ground and excited states in the QD may be influenced by the electric field originating from the crystal phase transition. To observe the effect of this transition it is necessary to make a comparison between experimental data obtained by time resolved photoluminescence (PL) and photoluminescence excitation (PLE) measurements and theoretical simulations.

The electronic structures and the variation of the fundamental gap of the Wz phase of these crystals have been only quite recently studied. Knowledge of the full band structure of the these polymorph materials is fundamental in order to possibly extract parameters for tight-binding [18] (TB) and k • p [17] model calculations. However, at present time only a relative small amount of parameters for bandstructure calculations are available in the literature.

Quasiparticle calculations of newly observed Wz polymorphs of InAs and GaAs has been reported in Ref. [START_REF] Zanolli | [END_REF] using the GW approximation. De et al. have presented calculations of the bulk electronic band structures of non-nitride III-V semiconductors in the Wz phase using empirical pseudopotentials including spin-orbit coupling[20], while Jancu et al. have studied the electronic structure of the Wz phases of GaAs and InP using TB model and calculate the optical properties of crystal phase quantum dots formed by thin Zb layers inserted in a Wz-InP NW [21]. Since the Zb and Wz structures differ only by the arrangement of the second neighbors of a given atom, they used the same set of parameters for the two phases, while the atomic positions are given as input data and lattice constants are precisely determined by x-ray diffraction.

Very recently quasiparticle self-consistent GW calculations have been used to analyze the band structure of GaAs in Wz phase. [22]. Here, also the Kohn-Luttinger and Rashba-Sheka-Pikus effective Hamiltonian band structure parameters of Wz GaAs have been determined from these first-principles band structure calculations.

In this work, we have simulated an AlGaAs NW with an embedded GaAs QD, considering two different cases, a pure Zb case and a polymorph Wz-Zb case, respectively. Dimensions, structural and geometric parameters have been strictly derived from experimental data. We have implemented an 8-band k • p model [START_REF] Lew Yan Voon | The k • p Method: Electronic Properties of Semiconductors[END_REF] for bandstructure calculations, and the parameters for the Pikus-Bir Hamiltonian of the Wz materials have been derived from the correspondent Zb materials using the cubic approximation as described by Chuang et al.[32], while the stiffness parameters of Wz crystals have estimated from Martin's transformations [27], [28] using the Zb parameters.

The paper is organized as follows. In the first section we present the experimental results that have inspired our numerical investigation. In section III we describe the models and the way we have derived all the parameters needed in the simulations, while all the numerical results are discussed in section IV. Finally, we give our conclusions.

II. EXPERIMENT

In these studies we have considered an AlGaAs NW with an embedded GaAs QD. AlGaAs NWs were grown on GaAs (111)B semi-insulating substrate using an EP1203 MBE system equipped with solid sources supplying Ga and Al atoms, and an As effusion cell to produce tetramers. The Al content along the NW and within the wire cross-sections was monitored by Raman spectroscopy in order to exclude possible inhomogeneities. We have found Al concentration in the range c = 0.24 -0.26. The growth was completed with deposition of GaAs to avoid possible oxidation of the AlGaAs shell layer. We did not perform any growth interruption at the heterointerfaces. The diameter of NWs in our samples is in the range 20-50 nm. Scanning electron microscopy (SEM) and transmission electron microscopy (TEM) characterizations were performed before the optical measurements.

The samples were excited by either second harmonic of Nd:YAG laser or doubled Ti:Sa laser with 100 fs pulse duration and 70 MHz frequency repetition. The average power of the lasers was about 10 µW, the spot diameter was 1.5 microns. Jobin-Yvon 30 cm monochromator and CCD detector were used for registration of the spectra. We used microscope to focus the laser beam and to collect the signal in the micro PL measurement of the single NW. The minimal diameter of the laser spot was 1.5 microns that allow us to extract the signal from a single NW. Additionally the set-up allows measure spectra with time resolution. The TEM image of the single AlGaAs NW with a slice of GaAs is shown in the insert of Figure 1. The dark area corresponds to the single GaAs QD. It is clearly seen that the interfaces along the axis of growth and along the radius are rather sharp. Considering cylindrical symmetry, the GaAs QD looks like a disk surrounded by the wide-gap semiconductor AlGaAs. The diameter of the NW in our samples was 25-50 nm, the length 1000 nm, and the thickness of the QD was of several mono-atomic layers, i.e., 2-5 nm.

Figure 1 shows PL spectra taken from a single QD embedded into a NW. At low excitation intensities of ≈ 1W/cm 2 a single narrow line is observed at the energy of 1.6525 eV. The width of this line is 0.2 meV and mainly limited by the spectral resolution of our spectrometer. It is natural to connect this line with the annihilation of an exciton in the single QD. This figure shows also the spectrum of photoluminescence excitation (PLE) taken from the same dot. The distance between the absorption (PLE) lines corresponds to the "transverse" carrier quantization in the dot of 15 nm in diameter. The energy of the "longitudinal" quantization is ≈130 meV, which corresponds to the QD thickness of 2.5 nm. In the PLE spectrum one can see very broad absorption lines at high energy; we connect these lines with carrier absorption in the NW itself.

In Figure2 a TEM image of one of our AlGaAs NW clearly shows the evidence of phase transitions from Zb to Wz crystal structures.

Here, layers in which the brighter lines apper parallel to the growth direction corresponds to ABAB lattice repetitions of Wz, whereas the layers in which the bright lines appear to zigzag correspond to the ABC repetitions of the Zb phase. Interestingly, the darkest layers between the two crystal phases suggest the presence of a strain field given by a lattice mismatch between Zb and Wz AlGaAs.

III. MODELS

A. Strain calculations

The continuum model we have used for our calculations of the electromechanical fields has been derived following the theory in Ref. [23]. The total mechanical and electrical free energy density change dU of a piezoelectric medium can be written as [24]:

dU = dU mech + dU elec = T dS + σ ik dε ik + E i dD i , (1)
where T , S, σ ik , ε ik , E i , and D i are the temperature, entropy, stress tensor, strain tensor, electric field, and electric displacement, respectively.

With isentropic conditions (dS = 0) and using crystal symmetry considerations the elastic energy for a crystal with Zincblende (Zb) symmetry read [25]:

U mech Zb = 1 2 C 11,Zb (ε 2 xx + ε 2 yy + ε 2 zz ) + 2C 12,Zb (ε xx ε yy + ε xx ε zz + ε yy ε zz ) + 4C 44,Zb (ε 2 xy + ε 2 xz + ε 2 yz ) , (2)
where C ij,Zb are the linearly independent stiffness constant parameters for a Zb structure.

Similarly, for a crystal with Wurtzite (Wz) symmetry we have for the elastic energy [25]:

U mech W z = 1 2 C 11,W z (ε 2 xx + ε 2 yy ) + C 33,W z ε 2 zz + 2C 12,W z ε xx ε yy + 2C 13,W z ε zz (ε xx + ε yy ) + 4C 44,W z (ε 2 xz + ε 2 yz ) + 2(C 11,W z -C 12,W z )ε 2 xy , (3)
where C ij,W z are the linearly independent stiffness constant parameters for a Wz structure. While the stiffness parameters for Zb Gas, AlAs and AlGaAs alloy can be found in literature [26], we do not have available the correspondent Wz parameters, so we have estimated them from Martin's transformations [27], [28] using the Zb parameters:

C 11,W z = 1 6 (3C
where

∆ ≡ √ 2 6 (C 11,Zb -C 12,Zb -2C 44,Zb ). (5)
The piezoelectric field is included directly in our fullycoupled continuum model [23], and the not available Wz piezoelectric parameters have been approximated using parameters of well known Wz materials, as GaN, AlN and AlGaN alloy.

Expressions from 1 to 5 have been used to derive the governing equations of the electromechanical fields of the crystal [23].

B. k • p method

We have used for band structure calculations 8-band k • p model [START_REF] Lew Yan Voon | The k • p Method: Electronic Properties of Semiconductors[END_REF], which describes electron, heavy-hole, light-hole, and spin-orbit split-off bands around the Γ point of the Brillouin zone and treats all other bands as remote bands. An 8×8 effective-mass Hamiltonian has been implemented, based on Foreman's application of Burt's exact envelope function theory to planar heterostructures [START_REF] Foreman | [END_REF].

The wave function of a state n with energy E n is given by a linear combination of the eight Bloch parts weighted by the respective envelope functions,

ψ n (r) = 8 i=1 φ i (r)u iΓ (r), (6)
where φ i are the envelope function and u i are the Bloch states at k=0 [31].

The Luttinger-Kohn parameters the Zb k • p Hamiltonian (L 1 , L 2 , M 1 , M 2 , M 3 , N 1 , and N 2 ) have been taken from Ref. [26], while the A i parameters in the Pikus-Bir model for Wz have been derived from the correspondent Zb materials using the cubic approximation as described in Ref. [32]:

A 1 = -γ 1 -4γ 2 , A 2 = -γ 1 + 2γ 2 , A 3 = 6γ 2 , A 4 = -3γ 2 , A 5 = -3γ 3 , A 6 = -3 √ 2γ 3 , (7)
GaAs Pikus-Bir parameters are in a good agreement with the corresponding parameters given by GW calculations in Ref. [22]. For the band offset of the GaAs Wz material we refer to Ref. [16], and we have assumed an equivalent offset for the AlAs Wz material.

Results for electromechanical fields with the continuum model are included in the framework of the k• p as it is usually done via deformation potentials [START_REF] Bir | Symmetry and Strain-Induced Effects in Semiconductors[END_REF].

The complete set of parameters for the Wz materials we have used for our simulations is given in TableI, with the same usual notation as in Ref. [26].

Parameters

GaAs Optical properties are calculated in the dipole approximation, starting from the dipole matrix elements µ nm , between states n and m, respectively. These are calculated from the momentum matrix elements which is represented in k-space as [START_REF] Lew Yan Voon | [END_REF]

p nm = m 0 ∂H nm ∂ k . (8)
Finally, oscillator strengths for the optical transitions are given by [36]:

P nm = 2πe 2 2 ǫ 0 m 2 0 V |ê • p nm | 2 , ( 9
)
where ê is the direction unit vector of the electric field of the linearly polarized incident light, V is the quantum dot volume, e is the electron charge, and ǫ 0 is the vacuum permittivity.

C. Geometry and structure

In our simulations we have considered an Al c Ga 1-c As nanowire (c = 0.25) with a diameter of 40 nm and height of 50 nm, with an embedded GaAs QD with a diameter of 20 nm and an hight of 5 nm. Two different cases have been studied and compared: a pure Zb case, where the whole structure is assumed to be grown in pure Zb crystal phase, and a mixed structure as shown in Figure 3, where the GaAs QD is sandwiched by two adjacent 8 nm high Wz layers and surrounded by a Zb AlGaAs shell. The mixed structure is completed by Zb layers at the top and at the bottom. We will refer in the following to this structure as the mixed (Mx) case.

IV. RESULTS

A. Electromechanical fields

The electromechanical fields (strain and piezoelectric fields) are evaluated as described in section III-A. In Figure 4 we plot the magnitude of the absolute value of the strain field for the pure Zb and the Mx structures, respectively.

It is interesting to notice that in the Mx case, besides the usual strain field inside and in the proximity of the QD, given by the lattice mismatch between the two materials, we can observe a weaker strain field at the interfaces between the two crystal phases of AlGaAs. This field is also observed in the TEM image shown in Figure 2. Despite its weakness, which was indeed expected from the small lattice mismatch between the GaAs and the Al c Ga 1-c As (c = 0.25) and the even smaller crystal mismatch between Zb-and Wz-AlGaAs, the strain field is non zero practically everywhere in the nanowire. This also affects the electric properties of the structure.

The electric field originating from the piezoelectric effect in the Mx case is relatively strong and presents an unusual shape, different from the typical results obtained for Zb QDs. This can be observed in Figure 5, showing the z component of the electric field for both Zb and Mx structures.

Instead of the usual quadrupole shape obtained for Zb QDs, the Mx case gives a Wz-like shape of the electric field close and inside the dot, with its typical dipole characteristics, and opposite signs inside and just outside of the top and the bottom of the dot. This occurs even though the central volume of the column is in a pure Zb crystal phase and owes to the stronger magnitude of the piezoelectric effect induced by the Wz layers, which leads to stronger a electric field. A stronger magnitude of the piezoelectric effect in Wz heterostructures compared to Zb heterostructures has been already reported in the literature [25].

B. Results for confined states

The quantum confinement in the QD is mainly given by the band offset between the two materials and between the two different crystal structures of the AlGaAs in the Mx case, since the impact of the observed strain fields on the confinement via the deformation potentials [START_REF] Bir | Symmetry and Strain-Induced Effects in Semiconductors[END_REF] is not particularly relevant. In Figure 6 we plot the conduction (E c ) and valence (E v ) band along the z direction for the considered structures, i.e., Zb and mixed. The usual case of a Zb QD embedded in a Zb matrix has been considered as a useful comparison in order to better understand the results for the energy levels. For the Mx conduction band, the well of the QD is surrounded by the energetically higher barrier of the Wurtzite layers, higher as compared to the lower potentials of the Zb AlGaAs. This results in an increased quantum confinement for the electron states. On the other hand, the hole bands of the Mx structure exhibit a sort of double well with the potential profile of the QD and the two adjacent Wurtzite layers acting as a second wider well, that can contain several confined states. The form of the potential has a notable impact on the hole wave functions.

In Figure 7 we plot all the confined states for the electrons and the holes for both Zb and mixed situations. All states are doubly degenerate for spin inversion due to time reversal symmetry. As far as electron states are concerned, the main difference we observe between the Zb case and the Mx case is a higher degree of quantum confinement in the latter, which slightly increases the energy levels and gives raise to two further confined states, i.e., states 11 and 12 are delocalized in the Zb model, while confined within the QD in the Mx case.

The shape and the symmetry of the electron wave functions do not significantly differ between the two models, implying that they are mainly dependent on the ZB crystal structure of the QD, and are barely affected by the electric field and its symmetries. However, the actual lack of knowledge of the real values of piezoelectric coefficients for the Wz phase of these materials makes uncertain the impact of the electric field on the symmetry of the wave functions.

The higher electron states are completely delocalized in the Zb layers adjacent to the Wz layers. These states are not optically active, since they have almost zero overlap with the hole wave functions confined in the QD.

The situation for the hole states is notably different and quite interesting. Comparing the curves plotted in the bottom part of Figure 7, it is clearly possible to see that, while the first confined states of the two models are very close in energy (also showing a similar trend), higher energy states (in absolute value) display different trends. The states of the Mx model are systematically lower in energy. The onset of this change corresponds to states that start to delocalize in the Wz layers, adjacent to the QD. The increased available volume explains their lower energy compared to the case of a pure Zb structure. This is confirmed by the contour plots of the probability density which we show in Figure 8. Here we present a slice of the xz plane, with x = 30 nm and z = 40 nm, showing basically the center of the nanowire with the dot and the two adjacent Wz layers. Figure 8 (left) shows the probability density of the holes ground state, which, as expected, is completely confined in the QD. Figure 8 (right) shows the probability density of the hole state number 32, according the bottom y labeling of Figure 7. This state is mainly confined in the Wz layers. A similar situation, with states confined in the volume formed by the QD and the two adjacent Wz layers, is observable for a considerable number of the higher energy hole states. All these states are optically active, since they have a non-zero overlap with the electron states confined in the dot.

C. Oscillator strengths results

We finally present results for the oscillator strengths, according dipole approximation. In Figure 9 we show the optical transitions.

Despite some resemblances in the distributions of the peaks, we can notice an effect which confirms what we have discussed in the previous section. We observe a blue shift of the lowest transitions of the Wz instance, namely for the ground state -ground state transition and the transitions between the more confined energy levels. This is clearly due to the higher energy of the electrons in the Mx case. But, on the other hand, for the peaks at higher energy we notice a red shift of the lines correspondent to the Mx models, due to the difference in energy between the less confined hole states, up to 20 meV as observable in Figure 7.

Finally, in Figure 10, the lowest range of the spectrum is compared with the positions of three peaks given by experimental results from a specific nanowire sample with approximately the same dimensions of our model.

Here, the peak correspondent to the lowest transition, namely the ground state -ground state transition, has been measured by a PL experiment, while the other two positions we show are given by a PLE data set from the same sample. Although our Mx model cannot be completely validated due to the fluctuation of measurement data from different samples, the agreement for this particular sample is very good.

V. CONCLUSION

We have presented numerical simulations for polymorph Wz-Zb cylindrical AlGaAs NW with an embedded GaAs QD. Firstly, the connection with experimental structures has been shown, and consequently suitable continuum elastic and 8 band k • p models have been derived. All the necessary numerical parameters have been determined by appropriate approximations, where not available in the previous literature.

The results for the strain fields are in a good agreement with the experimental data, showing, besides the usual field given by the lattice mismatch between different material, a not-zero field at the interfaces between the two crystal phase of AlGaAs.

The analysis of the confined energy levels, wave functions and optical transitions clearly show the differences between an uniform crystal phase and a polymorph heterostructure. The comparison of numerical results and experimental spectra for optical transitions comfort the validity of our polymorph numerical model.

Figure 1 .

 1 Figure 1. PL spectrum taken from a single QD at the excitation by 533nm laser line with intensity of 0.1 W/cm 2 . PLE spectrum taken from the same QD. Insert: TEM image of the GaAs QD in the AlGaAs NW.

Figure 2 .

 2 Figure 2. TEM image of an AlGaAs NW with wurtzites and zincblende phase alternations.

  Figure 3. A scheme of the simulated AlGaAs nanowire with the embedded GaAs quantum dot (QD) for Mx case. The different crystal phases of the correspondent layers are also indicated. All the dimensions are given in the text.

Figure 4 .

 4 Figure 4. The magnitude of the absolute value of strain field for a pure Zb structure (top) and for the Mx structure (bottom).

Figure 5 .

 5 Figure 5. The z component of the electric field in the xy plane in the middle of the QD (top), and in the xz planes (bottom), respectively. Left column Zb structure, right column Mx structure. All values in GVolt/m.

Figure 6 .

 6 Figure 6. Conduction (top part of the figure) and valence (bottom part of the figure) band along the z direction for both the Zb and the Mx cases The energy is expressed in eV, where the 0 is given by the Fermi level.

Figure 7 .

 7 Figure 7. The calculated states for electrons (top) and holes (bottom). The x scales enumerates the number of the correspondent state (different scale for electrons and holes), while the energy in the y axis is referred to the Fermi level as in Figure6.

Figure 8 .Figure 9 .

 89 Figure 8. Holes ground state (left) and holes state number 32 according the bottom y labeling of Figure 7 (right) in xz plane (x = 30 nm and z = 40 nm). Continuous lines indicate the edges of the QD (below) and the edge between the Wz and the Zb phase (top).

Figure 10 .

 10 Figure 10. Oscillator strengths in x directions for the Mx model compared with the positions of experimental transition results: PL for the first transition, PLE for the second and third ones.

  11,Zb + 3C 12,Zb + 6C 44,Zb )

-3∆ 2 /(C 11,Zb -C 12,Zb + C 44,Zb ) C 12,W z = 1 6 (C 11,Zb + 5C 12,Zb -2C 44,Zb ) + 3∆ 2 /(C 11,Zb -C 12,Zb + C 44,Zb ) C 13,W z = 1 6 (2C 11,Zb + 4C 12,Zb -4C 44,Zb ) C 33,W z = 1 6 (2C 11,Zb + 4C 12,Zb + 8C 44,Zb ) C 44,W z = 1 6 (2C 11,Zb -2C 12,Zb + 2C 44,Zb ) -6∆ 2 /(C 11,Zb -C 12,Zb + 4C 44,Zb ),

ACKNOWLEDGMENT

The authors would like to thank...

Daniele Barettin received the Master degree in Physics from the University of Rome "La Sapienza", Rome, Italy, and the Ph.D. degree in Applied Mathematical Models from the University of Southern Denmark. Since 2011 is Postdoctoral Researcher with the Department of Electronic Engineering of the University of Rome "Tor Vergata". He is been visiting researcher at the Ioffe Institute of Saint Petersburg, Russia. His research interests include the study of optoelectronic properties of nanostructures, analysis and optimization of organic electronic devices, and the study of complex systems in statistical mechanics. He currently is a researcher with the Department of Electronic Engineering of the University of Rome "Tor Vergata". His research interests are in electronic transport optical properties of nanostructured devices and implementation of multiscale/multiphysics simulation models.

PLACE PHOTO HERE

Aldo Di Carlo received the Laurea degree in physics (cum laude) from the University of Rome "La Sapienza", Rome, Italy, and the Ph.D. degree from the Technical University of Munich, Munich, Germany. In 1996, he became a Research Assistant in the Department of Electronic Engineering, University of Rome "Tor Vergata", Rome, where he is currently an Associate Professor, since November 2001 and a Coordinator of the Center for Hybrid and Organic Solar Energy. His current research interests include the study of electronic and optical properties of nanostructured devices, their analysis and optimization and the fabrication of organic electronic devices, the charge transport in nanostructured devices, high electron mobility transistors, organic thin-film transistors, molecular devices, and carbon nanotube-FETs. Dr. Di Carlo has been an Invited Speaker at the International Conference for the development of the nonequilibrium theory for the microscopic description of the transport process in nanostructured devices.


Nicolas Lesaffre

Jean-Jacques Sinou

Fabrice Thouverez

J-J Sinou

Contact Analysis of a Flexible Bladed-Rotor

Keywords: Rotor dynamics, flexible bladed-rotor, stability analysis, Campbell diagram, rotor/stator contacts, rotating loads

may come

Contact analysis of a flexible bladed-rotor

I. Introduction

Nowadays, more than ever, turbo machinery designers seek for increase the efficiency of their machine what could result in a loss of weight and a decrease in fuel consumption. The efficiency of turbo-generators depends strongly on the clearance between the rotating and stationary parts (Song and Martinez-Sanchez, 1997-a, Song and Martinez-Sanchez, 1997-b, Childs, 1993[START_REF] Ehrich | Handbook of Rotordynamics[END_REF][START_REF] Yamamoto | Linear and Nonlinear Rotordynamics: a Modern Treatment with Applications[END_REF]: the wider the clearance, the less efficient the machine. Hence it is desirable to reduce the clearance by as much as possible. Unfortunately, it cannot be reduced below a certain minimum level due to uncertainties in the dynamics, fluid dynamics (Muszynska andBently, 1996, Ehrich, 1993), thermal loads [START_REF] Goldman | Rotor-to-Stator, Rub-related, Thermal/Mechanical Effects in Rotating Machinery[END_REF], machining tolerances and material properties [START_REF] Maozhong | Erosion Wear Behaviour and Model of Abradable Seal Coating[END_REF][START_REF] Marscher | A Phenomenological Model of Abradable Wear in High Performance Turbomachinery[END_REF]. The more accurately these parameters are known, the smaller the clearance that can be made. However, reduced clearances increase the possibility of violent contacts between rotor and stator. Under certain operating conditions, depending for instance on the dynamic properties of the rotor or/and the stator and on the rotational speed, the contacts may lead to an unstable rotor-stator system [START_REF] Feng | Rubbing Phenomena in Rotor-Stator Contact[END_REF]Zhang, 2002, Edwards et al., 1999). A common interest of the rotating-machinery industry is to better understand the dynamic behaviour of rotating structures and in particular that of flexible bladed-rotors whose dynamic behaviour differs from those usually studied in rotor-dynamics especially because of the dynamics of the flexible blades that can lead to instabilities [START_REF] Crandall | Rotordynamics[END_REF]Dugundji, 1981, Genta, 2003). The study of such structures is quite complex and must be performed using numerical models (Gmür, 1997, Lalanne and[START_REF] Lalanne | Rotordynamics Prediction in Engineering[END_REF]. To obtain a closed form solution suitable for stability analysis, it is possible to resort to simplified models [START_REF] Sinha | Dynamic Characteristics of a Flexible Bladed-rotor with Coulomb Damping Due to Tip-rub[END_REF][START_REF] Al-Bedoor | Stability Analysis of Rotating Blade Bending Vibration Due to Torsional Excitation[END_REF][START_REF] Turhan | Dynamic Stability of Rotating Blades (Beams) Eccentrically Clamped to a Shaft with Fluctuating Speed[END_REF]. It is fundamental for turbo machinery industries to understand the leading phenomenon in case of contacts between the blades of the rotor and the casing and some researches have been recently investigated [START_REF] Schmiechen | Travelling Wave Speed Coincidence[END_REF][START_REF] Arnoult | Modelisation Numerique et Approche Experimentale du Contact en Dynamique: Application au Contact Aubes/Carter de Turboreacteur[END_REF][START_REF] Legrand | Modèles de prédiction de l'interaction rotor-stator dans un moteur d'avion[END_REF]. This paper presents firstly a simplified model of a flexible fully-bladed rotor constructed in the rotating frame by using an energetic approach and taking into account the basic effects of rotating machines. In a second part, an elastic ring having only one mode shape is modelled also in the rotating frame by using an energetic method. Its behaviour when excited by rotating loads is studied. Finally, the contact problem of the blades of the rotor with a casing modelled by a flexible ring is investigated. The results and the phenomenon underlined are compared to the latter ones.

II. Model formulation

The rotor considered in this study is shown on Fig. 1. The model is inspired from an existing model by Sinha [START_REF] Sinha | Dynamic Characteristics of a Flexible Bladed-rotor with Coulomb Damping Due to Tip-rub[END_REF] but developed in the rotating frame by using an energetic method. The rotor is composed of a shaft modelled by an Euler-Bernoulli beam, connected to a rigid disk modelled by a concentrated mass with rotational inertia. This shaft is set on bearings at multiple locations. On the rigid disk is clamped a full set of flexible blades also modelled by Euler-Bernoulli beams. In the rotating frame, two degrees of freedom are considered for the shaft: two orthogonal translations ) , ( t z u and ) , ( t z  in the disk's plane (see Fig. 2), and one degree of freedom for each blade defining its deflection ( , ) j s t  .

A Rayleigh-Ritz approximation is used to express the degrees of freedom of these different parts. Thus, they are expressed by a sum of shape functions multiplied by time-dependent coefficients:

1 ( , ) ( ) ( ) ( ) tot m o m m m u z t U t U t W z     (1) 1 ( , ) ( ) ( ) ( ) tot m o m m m z t V t V t W z      (2) 1 ( , ) ( ( )) ( ( )) tot n j n j n j n s t X t Y s     (3)
where z is the axis of the shaft and s , the axis along the blade. In these expressions, ) (t Uo and ) (t Vo are rigid-body translations of the shaft. tot m and tot n are the number of modes considered to express its motion and the flexure of the th j blade respectively. In this Rayleigh-Ritz approach, the shape functions have only to verify the geometric boundary conditions of the problem. Due to the fact that all blades are clamped in the rigid disk, the shape functions of these latter ones must verify

0 ) 0 (  n Y and 0 ) 0 '(  n Y
, where prime denotes differentiation versus space coordinates. Thus, the chosen expression for the shape function of the blade deflection is:

( ) sin( ) n n n Y s a s s    (4) with ) 2 ( ) 1 2 ( L n n    
, and n n a    . Concerning the shaft, since it is supported by bearings, its shape function has no geometric boundary conditions to verify and its chosen expression is:

) cos( 1 ) ( z z W m m    (5) with ) 2 ( ) 1 2 ( l m m    
. However, it could be noticed that these functions verify

0 ) 0 (  m W and 0 ) 0 '(  m W
thus no motion is permitted at the right end of the shaft (i.e. at z = 0). Since an energetic method is used to develop this model, energies and potentials have to be defined for the shaft and for the blades as well. Thus, the kinetic energy blade T of the th j blade, located at an angle

tot j N j   2 
in this frame (as indicated in Fig. 2) where tot N is the total number of blades, can be fully expressed by the following relation:

0 0 0 2 / / / 0 0 1 1 ( ) 2 2 L L T blade b b S R S R S R T S V G ds ds        I    (6)
where b  and b S are the density and the area of a blade cross section respectively.

) ( 0 / G V R S  , 0 / R S  
and I are the speed, the rotation of the mass center of a blade cross section (in relation to the fixed frame) and its inertia matrix, defined in its inertial frame. It may be noticed that the complete expression of the kinetic energy should be considered because in certain extreme geometric conditions, an approached kinetic energy can lead to unrealistic instabilities of the rotor. The potential energy blade

  L b b b d ds t s I E F blade 0 2 )) , ( ' ' ( 2 1    (8)
Next, a pre-stress potential taking into account the centrifugal stiffening effects is defined by:

2 2 2 2 0 1 ( ) ' ( , ) 2 2 blade L g b b R s r S s t d s               (9)
In the same way, the kinetic energy shaft T (Lesaffre et al., in press), the potential energy shaft int  associated with the elastic deformation of the shaft:

2 2 int 0 1 ( '' ( , ) '' ( , )) 2 shaft l s s E I u z t z t dz      ( 10
)
where s E and s I are the elastic Young's modulus of the shaft material and the shaft area moment of inertia for flexure respectively, are defined, as well as a function of dissipation (Lesaffre et al., in press). Then, Lagrange's equations are used to obtain the system of equations of the dynamic behaviour of the full flexible bladed-rotor. This system of equations can be written under the following form: and N rotor are the mass matrix, the damping matrix, the gyroscopic matrix, the stiffness matrix, the centrifugal stiffening matrix and the spin softening matrix, respectively. The vector rotor X defining the generalised degrees of freedom vector of the system contains 2

    M X D G X K K N X
2 tot tot tot m n N  
elements and has the following expression:

0 0 1 1 11 21 1 1 . . . . . . . . . . tot tot tot tot tot tot T rotor m m n N n N U V U V U V X X X X X      X (13)
Since the model has been fully developed in the rotating frame there are no time-dependent terms resulting from the periodicity of the rotating structure, in the analytical formulation and this is particularly useful to study such rotating structure. For instance, the eigenvalues and the stability of the flexible bladed rotor can be investigated by determining the solutions ib a   of the characteristic equation: are the generalised damping and stiffness matrices for the stator respectively. The rotor becomes unstable if any one or more of the eigenvalues' real part a is positive whereas the eigenvalues' imaginary parts correspond to the eigen frequencies of the system. So, the evolution of these eigenvalues' imaginary parts as a function of the rotational speed of the rotor give a Campbell diagram in the rotating frame. Figure 3 . In the rotating frame, the negative slope curves are forward modes (precession motions in the same sense of rotation as the own rotation of the shaft) and the positive slope ones are backward modes. The evolutions of the shaft's eigen-frequencies, shown on Fig. 3(a), are mainly due to the gyroscopic effects. The evolution of the blades' eigen-frequencies (starting at 31 Hz) i.e. their stiffening, is due to the centrifugal effects included in the pre-stress potential blade g  . Two loci separation phenomena [START_REF] Pierre | Mode Localization and Eigenvalue Loci Veering Phenomena in Disordered Structures[END_REF] can be seen as illustrated by the points A and B: when two eigenvalue loci approach each other, they either cross or do not cross; often in the latter case, even though the loci nearly intersect, in fact they do not but rather veer away from each other. Fig. 4 represents the mode shapes associated with the veering phenomenon A. Before this phenomenon, the mode shape associated with the lower frequency (among both concerned mode shapes) consists mainly in a shaft mode having one node (see Fig. 4(a)) and the one associated with the higher frequency consists also mainly in a shape mode but, having two nodes (see Fig. 4(b)). After the veering, the mode shape associated with the lower frequency consists in a shaft mode having two nodes (see Fig. 4(d)) where as the mode associated with the higher frequency consists in a shaft mode having one node (see Fig. 4(c)).During these veering phenomena, mode shapes and sense of rotation are switched between the eigenvalues that veer away from each other. Moreover, the stability of the system can be known by correlating these curves to these of the Fig. 3(b) representing the evolution of the eigenvalues' imaginary parts but as a function of the associated eigenvalues' real parts. In this case, it can be seen that the system is perfectly stable because there is no eigenvalues having a real part positive.Such stability analysis can be performed for several rotating systems and stability maps can also be obtained [START_REF] Lesaffre | Stability Analysis of a Flexible Blade-Rotor[END_REF].

det( ² ) 0 rotor rotor rotor      M C K   (14) R r L  Z -axis l tot j N j   2  L ) , ( t s  s r Rigid disk

III. Rotor-stator contact a. Simple model of rotating loads on the elastic casing

Prior to the study of the contact between the blades of the rotor and the casing, a simple model of an elastic ring excited by rotating loads is developed in the rotating frame by using also an energetic method. This model consists in p loads, composed of a spring (with stiffness j k ) and a mass j m located at the angles j  in the rotating frame, rotating permanently on an elastic stator without rubbing (Cf. Fig. 5b). The only radial degree of freedom of these spring-masses are thus the same as the radial displacements s u of the ring at j  . Their kinetic energy and their potential energy are defined:

    2 2 2 1 1 ( , ) ( , ) 2 p loads j s stat s j j T m u t R u t                  (15)   2 1 1 ( , ) 2 p loads j s j j k u t          ( 16
)
The casing is modelled by a flexible circular beam (see Fig. 5a). The in-plane flexural vibrations of the ring are considered. Thus, two degrees of freedom are considered in the rotating frame: the radial displacement ( , )

s u t 
and the tangential displacement ( , ) w t  (see Fig. 5). This latter one can be expressed by [START_REF] Love | A Treatise on The Mathematical Theory of Elasticity[END_REF]

: 2 ( , ) ( )cos ( )sin tot k n n n w t A t n B t n        (17)
where the rigid body motion has been eliminated. The considered stator is supposed to be inextensible, implying thus that the radial displacement can be expressed from the tangential one by:

    , , s w t u t       (18)
The stator is also fully described in the rotating frame. Its kinetic energy stat T , as well as its potential energy stat  are defined:

2 2 2 1 ( , ) ( , ) ( , ) ( , ) 2 t s stat stat stat s stat t u w T S u t t w t t R d                                               (19) 2 2 3 1 ² ( , ) ( , ) 2 ² t stat stat s stat stat t E I u w t t d R                       (20)
In these expressions, stat  , stat S , stat E and stat I are the density, the cross section area, the elastic Young's modulus of the ring material and the area moment of inertia for flexure of the casing respectively. It must be noted that the present ring is not damped. For simplicity reason, only one mode shape of the casing (the th n one) is considered at a time thus the dynamic behaviour of the system is governed by a 2 2  matrix equation. The matrix equation of the dynamic behaviour of the ring excited by p rotating spring-masses, without rubbing, can be obtained by using Lagrange's equations and has the following expression: The stability of the casing can be studied by determining the eigenvalues of this matrix equation. The influence of the stiffness, of the mass and of the number of the moving loads can be studied. For instance, in the particular case of only one rotating radial springmass on the casing, the matrix equation of the dynamic behaviour of the system is the following one:

    2 2 2 2 1 1 2 2 2 2 1 1 1 sin ( ) sin( )cos( ) sin( )cos( ) 1 cos ( ) p p stat stat stat j j j j j j j n p p n j j j stat stat stat j j j j S R n m n n m n n n A B m n n n S R n m n n                                                       2 2 0 2 1 2 1 0 stat stat stat n n stat stat stat S R n n A B S R n n                                 2 2 2 2 2 2 3 2 2 1 2 2 2 1 2 2 2 2 2 2 3 2 2 1 ( 1) 1
( ) sin( ) cos( ) ( ) s i n( ) ( 1) 1 (
E I n n S R n n R k m n n n k m n n E I n n S R n n R k m n n n k                                      2 2 2 1 ) cos ( ) n n p j j j A B m n n                                      (21) 2 1 2 1 sin( ) cos( ) p j stat j j p j stat j j m R n n m R n n                         (
    2 2 2 1 0 0 1 stat n n stat M n A B M n mn                            2 2 0 2 1 2 1 0 stat n n stat M n n A B M n n                               2 2 2 2 2 2 2 2 2 2 2 2 2 2 ( 1 ) 1 0 0 ( 1 ) 1( ) stat stat n n stat stat K n n M n n A B K n n M n n k m n                           2 0 stat mR n          (22) with stat stat stat stat M S R    and 3 stat stat stat stat E I K R  
It can be seem that the stiffness matrix can become singular and thus the system is likely to experience some unstable phenomenon [START_REF] Crandall | Rotordynamics[END_REF] as shown on Fig. 5 representing a stability analysis of the two nodal diameter mode shape of the stator excited by one rotating spring-mass. Two kinds of instabilities appear: divergence (instability at zero frequency) of the forward mode shape and mode coupling (the forward and the backward mode shapes of the stator become two mode shapes having the same eigen-frequency but one of them being X-axis

(rotating frame) stat R  w  s u  Y-axis (rotating frame) s u  j m j k  1 j m  1 j k  1 j m  1 j k  w 
stable and the other, unstable). The critical rotational speeds where the system diverges can be expressed analytically through the Routh-Hurwitz criterion applied to the characteristic polynom 2 ( )

P x Ax Bx C     instead of   4 2 2 ( ) det M G K P s As Bs C s s      
, in this simple case of only one rotating spring-mass (in the latter equation, M , G and K are the matrix written in Eq. 22). According to this criterion, the polynom 2 ( )

P x Ax Bx C    
has all its roots with real parts negative if A , B and C have the same sign. Thus, it turns out that :

for a light mass i.e.

     (23) with     2 2 2 2 2 2 1 1 1 1 c c stat stat m M n M n m          (24) and   2 2 2 2 4 2 2 1 1 n stat stat stat c stat stat stat n E I S R n n        (25) in the case of 2 2 c k m     or, if 2 2 2 2 c c      in the case of 2 2 c    .
if the spring-mass has a heavy mass i.e. 2 stat M m  , the system can also experience divergence if :

2 2 3 c    (26) with           2 2 2 2 3 2 2 2 2 2 1 1 2 1 1 2 1 stat c c stat stat M n n m k m n M n m n M n                     (27)
In these expressions, c  is the wave propagation speed in the stator for its th n mode shape and so, its th n critical rotational speed in the rotating frame. Thus, it appears that both the mass and the stiffness of the spring-mass can make the system diverge. In a general manner, the bigger the stiffness is, the larger the speed range where the system can experience divergence is. For strong stiffness, this speed range begins at c  and closes far away at 2 c  .

Concerning the mass, the heavier it is, the earlier the speed range where the system can experience divergence is. The critical speed where the system experiences mode coupling cannot be determined with the Routh-Hurwitz criterion as for the divergence instabilities (in the general case, the real part of 2 x s  is not the same as s ). However, this critical speed can also be determined analytically in the particular case of only one rotating spring-mass. As a matter of fact, a sufficient condition for the apparition of flutter, in this particular case of an undamped structure, is to have two eigenvalues with real parts null and opposite imaginary parts. These conditions can be satisfied if the discriminant of ( ) P x  is null. Thus, the critical speed of mode coupling for the system satisfy 2 4 0 B AC   . In the particular case of only one rotating spring-mass, the rotational speeds satisfying the latter equation can be expressed by :

  1 2 4 mc m mn             (28)
providing that the square roots can be defined, with:

    2 2 2 2 8 1 1 stat M n m n            (29)        2 2 2 2 2 2 4 1 2 4 1 stat stat stat c M n K n k mn K n k n m              (30)           2 2 2 2 2 2 2 2 2 1 1 2 1 stat stat stat stat c M n n m M n n k K K mn mn k n                   (31)   2 3 3 stat stat stat stat E I K n n R    (32)
The following figures illustrate these theoretical results. Figure 7 shows a stability analysis of the system in the case of a light mass 2 stat M m  . Figure 7b shows that only the mass is responsible of mode couplings. Figures 7c and7d show that the stiffness seems to move back the critical speed where mode couplings appear. Figure 8 shows also a stability analysis for the same system but in the case of a heavy mass . However, this latter divergence cannot be seen on Fig. 8 because mode couplings occur before. Thus, this shows (with Fig. 7b and7c) that the heavier the mass is, the earlier the system experiences mode couplings (at least in the mass range studied here).

The last parameter study considers the number of moving loads. Figures 9a,9b, 9c and 9d represent the real parts and the eigen-frequencies of the two nodal diameter mode shape of the stator with respectively one, two, three moving loads separated from 60° to each others (like blades , ,  on Fig. 11) and three moving loads, two separated from 180° and the third one at 60° from one of the latter two (like blades ,  and  on Fig. 11). These pictures have been made in the case of light masses 10 m kg  and stiffness 1 1 7 . k e N m   . In a general case, the influence of the number of moving loads on the stability of the stator is not obvious. However, in some particular configurations of the moving loads in comparison to the mode shape considered, the divergence instability can disappear as shown on

                           
) to avoid divergence instability is:

                                       (33)
Thus, for instance, in the case shown on Fig. 9c, all the conditions are satisfied: the masses are light, all the stiffness and masses are the same and:

2 2 2 2 2 2 2 3 2 3 6 6 6 6 6 6
3 sin (2 ) sin ( 2) sin ( 2) cos ( 2) cos ( 2) cos ( 2) 2 2 2 3 3 sin(2 )cos( 2) sin( 2)cos( 2) sin( 2)cos(2 ) 0 6 6 6 6 6 6

                            (34)
what is not the case of the Fig. 9d

                               (35)
The Routh-Hurwitz criterion shows also that if

  2 2 2 2 2 2 s i n( ) 2 1 j j j stat j j stat n m n M m M n                          
  then, the system can experience divergence for :

              2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 s i n( ) s i n( ) 2 1 2 s i n( ) 2 s i n( ) 1 2 1 1 1 1 stat j j j j j stat j j j j j c j j j j j j j s t a t stat j j j stat stat M n k n k n m n M n n m n m n n m n n m M n M n m M n M n                                                                           2 2 stat M                (36)
Figure 10 shows the same stability analysis as Fig. 8c but with three masses of 1000kg . The system would have experience divergence from 786RPM   but since it has be seen that the heavier the mass is, the earlier the system experiences mode couplings, the system is unstable from 622RPM   but does not diverge. Since some phenomenon occurring when an elastic structure is excited by rotating loads have been analysed and understood, the study of the contact between the blades of the rotor and the flexible casing of a turbo-machine can be performed. The contact between the flexible blade tips and a flexible casing is introduced by assuming no rubbing between the two structures.

In this section, the casing considered is the same as in the latter section but with internal viscous damping (with coefficient stat  ) introduced through a function of dissipation

stat d F : 2 2 3 2 2 3 3 2 2 1 ( , ) ( , ) ( , ) ( , ) 2 stat t stat stat s s d s t a t stat t E I u u w w F t t t td R                                      (37)
and with more than one mode shape considered. If tot k modes shapes are considered to express the ring displacements, its dynamic behaviour in the rotating frame can be described, as for the rotor, by a set of equations that can be written under the following form:

stat stat stat stat stat stat    M X C X K X 0     (38)
where stat M , stat C  and stat K  are the mass matrix, the generalised damping matrix, including the gyroscopic matrix, and the generalised stiffness matrix, including the spin softening matrix, for the stator respectively. The vector stat X which defines the generalised degrees of freedom vector of the stator contains 2 2 tot k  elements and has the following expression: 2 2 . . .

tot tot T stat k k A B A B      X .
The contact problem is treated by using Lagrangian multipliers j  . A contact energy j E  is defined for each contact at the th j blade located at the angle j  in the rotating frame:

  ( , )cos ( , )sin ( , )

j j j j s j E u l t l t u t clearance             ( 39
)
This energy is proportional to a Lagrangian multiplier (one per contact) and to the clearance between the blade tip and the casing. Thus there is one more equation per contact and also one more second member per contact in the behaviour equation of the system. Since the contact energy depends also on the radial displacement of the ring, the degrees of freedom of both structures are coupled. The contact between the rotor and the casing is supposed to be due to an unbalance strength. The dynamic behaviour of the whole structure can be described through the following matrix equation:

(40)

The rotational speed where the first contact appears has been placed under both the critical rotational speed of the rotor and of the stator by increasing the unbalance mass. In the rotating frame, this non linear problem is a static one. The balanced static contact configurations of the system are found by iterations: the first configuration of the rotor is found by resolving:

rotor rotor unbalance  K X F  (41)
If the clearance is consumed, the contact blades/casing exists and a new balanced static position is calculated by resolving the static part of Eq. 40:

(42) This non linear system has several balanced contact configurations and the evolution of the load acting on the system will influence the balanced configuration adopted. Thus, by increasing the rotational speed of the rotor, the evolution of the clearance between each blade tip and the casing can be followed as a function of the its rotational speed as shown in Fig. 12, in the particular case of a rotor having six blades (see Fig. 11). It can be seen that the first blade to touch the casing is the blade , then successively the blades ,  and , when the rotational speed increases. The associated deformed shape, at 164 RPM, is plotted on Fig. 13(a). When the rotational speed increases again, the blade  touches the stator. Thus, all the blades are in contact with the casing apart from the blade . The deformed shape associated, at 286 RPM, is plotted on Fig. 13(b). The system keeps this configuration while the rotational speed increases, until 310 RPM corresponding actually to the three nodal diameter mode shape critical speed of the casing. The associated deformed shape, at 310 RPM, is plotted on Fig. 13(c). From this point, the algorithm diverges. It appears that up to this rotational speed, all the balanced static contact configurations found are stable ones but after this speed, the system has only unstable balanced static contact configurations. Moreover, from about 700 RPM, the system has no longer balanced static contact configurations. Figure 14 shows a synthesis of these three rotational speed ranges. It could have been thought that between 310 RPM and 700 RPM, the system might have stable dynamic contact configurations but, all the time integrations performed diverged on the three nodal diameter mode shape of the casing. This phenomenon can be explained by using the simple model of rotating loads on an elastic structure studied in the previous section.

Effectively, it has been shown that an elastic ring excited by rotating loads could experience two kind of instabilities: divergence instabilities and mode couplings. In order to know if the divergence of the rotor coupled with the flexible casing is due to this phenomenon, the latter model of rotating loads on the elastic ring can be used. To compare the results of both models, the loads exciting the ring must have the same modal parameters than the rotor hence, its modal parameters have to be evaluated. For simplicity, this will be done at 0RPM  

. If the rotor is submitted to a constant strength, its first modal mass an modal stiffness can be expressed by: 1

1 1 X M X T rotor m   and 1 1 1 X K X T rotor k  
, 1 X being its first eigen vector and 1 X  , its first eigen vector normalized so that the displacement of the rotor at the point of application of the strength and in its direction could be unity. The modal parameters for the rotor estimated like this are . The latest stable balanced contact configuration (i.e. just before the system diverges) obtained statically and by time integration consists in five contacts: at blades , , , , and . Figures 15a andb represent the comparison between respectively the results (real parts and eigen-frequencies) obtained with the three nodal diameter mode shape of the simplified model having five moving loads configured like blades , , , , and  with 1 m and 1 k as parameters and, the Campbell diagram and the associated decay rate plot for the flexible casing in contact at blades , , ,  and . The importance as well as the influence of the others mode shapes of the stator on each others appears through all the others instabilities that can be seen on Fig. 15b and not on Fig. 15a. However, these figures confirm that the stator experiences divergence through its three nodal diameter mode shape just after 310 RPM which is the critical speed c  of this mode shape. So it can be said now that the instability of the three nodal diameter mode shape of the stator results from rotating loads on this elastic structure. Moreover, with this configuration of rotor i.e. with six blades separated by 60° from each others, there is no contact configuration that satisfies the four conditions to avoid the divergence of the three nodal diameter of the stator. Thus, until the dynamics of the blades is not introduced during contact i.e. until the balanced contact configurations consist in steady state contact, the blades will act as rotating loads on the stator and this one will always experience divergence on its three nodal diameter mode shape.

Conclusion

A flexible fully bladed-rotor has been modelled in the rotating frame using an energetic approach. The gyroscopic effects, the spin softening effects and the centrifugal stiffening effects have been taken into account so from a phenomenological point of view the system is almost complete. Moreover, the formulation obtained in the rotating frame appears to be useful to study such rotating structures, especially in case of contacts between the flexible blades and the casing. This very important problem has been first studied through a simple model of an elastic ring having only one mode shape and excited by rotating loads. This study has underlined that the considered structure can experience two kinds of instabilities: divergence instabilities and mode couplings. The influences of the loads' parameters have been analysed and in some simple cases, the critical speeds associated have been determined analytically. Then, the casing of a turbo-machine has been modelled by an elastic ring too and the contact between the rotor of the machine and this casing has been studied. This problem has been solved, in part, by calculating the succession of the balanced static contact configurations of the system during an increase in its rotational speed, by using a Lagrangian multipliers method. It turns out that in the particular case of the rotor considered in this study, the system has several balanced configurations and adopts stable static configurations as a function of its rotational speed until the three nodal diameter mode shape critical speed of the casing where it always diverges. From this rotational speed, the system has neither stable balanced static contact configurations nor stable dynamic contact configurations. The latter simple model of an elastic ring excited by rotating loads used with the same parameters as the rotor has shown that this phenomenon results indeed from the action of the blades of the rotor acting like rotating loads on the casing.

F

  elastic Young's modulus of the blade material and the blade area moment of inertia for flexure respectively, is defined, as well as a function of dissipation blade d associated with the internal damping (viscous damping with coefficient b  ) of the th j blade:

Figure 1 Figure 2

 12 Figure1Model of flexible bladed-rotor[START_REF] Sinha | Dynamic Characteristics of a Flexible Bladed-rotor with Coulomb Damping Due to Tip-rub[END_REF]

  (a) represents a Campbell diagram for a flexible twenty eight-bladed-rotor of one meter in length, supported by three isotropic bearings having radial stiffness and viscous damping coefficients

  Fig. 4. Mode shapes of the rotor associated to the veering phenomenon A on Fig. 3, (a) and (b): before the veering, (c) and (d): after the veering.

Figure 5

 5 Figure 5 (a) Flexible casing, (b) model of moving loads on an elastic stator without rubbing

  Fig. 8c. In our case and according to the Routh-Hurwitz criterion, a sufficient condition for the stator excited by several rotating loads with light masses (

Figure 6

 6 Figure 6 Stability analysis of the two nodal diameter mode shape of the stator excited by one rotating spring-mass with 10 m kg  and

Figure 7 Figure 8 Figure 9

 789 Figure 7 Stability analysis of the two nodal diameter mode shape of the stator excited by one rotating spring-mass with 2 stat M m  a)

Figure 10 Figure 13 Figure 15

 101315 Figure 10 Stability analysis for the same system as shown by Fig. 8c but with masses of 1000kg

D rotor , G rotor , K rotor , K rotor geom

  0

	rotor rotor 		rotor		rotor		rotor		rotor		rotor geom		rotor	rotor		(12)
	where rotor

M

,