Igor S Nefedov
email: igor.nefedov@tkk.fi

Andrey S Soloviev

Anne-Claude Tarot

Wafa Abdouni

S Zouhdi

Surface waves in a magnetized ferrite slab filled with a wire medium

Keywords: . 41.20.-q Applied classical electromagnetism, 42.70.Qs Photonic bandgap materials, 75.30.Ds Spin waves

Novel metamaterial, based on wire medium embedded into magnetized ferrite, is studied. Waves in unbounded ferrites filled with wire media, surface wave at the interface of this metamaterial and the air as well as waves in a ferrite slab adjacent to a wire medium are considered. Different geometries of wires arrangement and different magnetization directions are discussed. Effective permeability was introduced for the case where both plasma and magnetic properties take place. Dispersion diagrams and applicability of the Drude model for the description of the wire medium in a host matrix, possessing high permittivity and permeability, are discussed.

Introduction

In this paper we discuss electromagnetic properties of new metamaterial, based on the magnetized natural magnetics (ferrites) and the wire medium (WM) (the medium formed by a lattice of ideally conducting parallel wires, see [1,2]) If the wires are thin (compared with the lattice period), such the medium can be homogenized and its effective relative permittivity dyadic can be written as (wires are oriented in z direction) ǫ = ǫ h (u x u x + u y u y ) + ǫ z u z u z (1) where ǫ z is expressed by the "plasma" (the Drude) formula

ǫ z = ǫ h 1 - ω 2 p ω 2 ǫ h µ h = ǫ h 1 - k 2 p k 2 ǫ h µ h . (2)
Here ǫ h and µ h are permittivity and permeability of the host medium, respectively, k = ω/c √ ǫ h = k o √ ǫ h µ h , c is the speed of light, the constant ω p (or corresponding k p ) is the equivalent "plasma frequency" for the case when the host medium is vacuum that gives grounds to call a wire medium as an "artificial plasma". (There exist different plasma models and we will use the formula for k p obtained in [3]). If the wavevector in WM has nonzero component along the wires, the plasma model Eq. 2 should be corrected introducing the spatial dispersion, see [START_REF] Shvets | Advanced Accelerator Concept: Tenth Workshop[END_REF][START_REF] Belov | [END_REF]. It follows from Eq. 2, effective permittivity of WM can be made smaller than the permittivity of the host medium, even negative. Thus, embedding WM into ferrite, we can obtain material exhibiting high permeability and low effective dielectric constant. In next sections we will study wave properties of the ferrite -wire medium (FWM) as a PBG structure and examine the applicability of the effective medium theory for the description of FWM.

The paper is organized as follows: in Sec. 2 we present the effective medium model for waves propagating in unbounded transversally magnetized ferrite host matrix, filled with WM. In Sec. 3 the full-wave model, based on Green's function method, for the same structure is described. Dispersion diagram of FWM, considered as PBG structure, is discussed. Applicability of the effective medium model is demonstrated using results of Green's function method. In Sec. 4 surface waves in several structures containing ferrites and WM, whose geometries allow description in the framework of effective medium theory, are discussed.

Effective medium model

Let us consider a two-dimensional periodic (in x and y directions) lattice, formed of conductive, infinitely long wires, embedded into a magnetized YIG host matrix. Directions of the magnetic bias field H 0 and the wave vector k are shown in Fig. 1. The relative permeability tensor takes form:

µ =     µ -iκ 0 iκ µ 0 0 0 1     (3)
where

µ = 1 + ω m (ω 0 + iαω) (ω 0 + iαω) 2 -ω 2 , κ = ω m ω (ω 0 + iαω) 2 -ω 2 ω m = 2πγM 0 , ω 0 = 2πγH 0 , α = 2πγ∆H 2ω 0 .
Here γ = 2.8 GHz/KOe is the gyromagnetic ratio, M 0 is the saturation magnetization, H 0 is the biasing magnetic field and α is the damping factor, where ∆H is the ferromagnetic resonance linewidth. The value ω 0 is referred to as the ferromagnetic resonant frequency. For an extraordinary wave, propagating transversely to H 0 (in y direction for our case), the effective permeability can be introduced:

µ eff = µ 2 -κ 2 µ . (4)
One can show that µ eff → ∞ if ω → ω ⊥ = ω 0 (ω 0 + ω m ) in the lossless case (α = 0), where ω ⊥ is referred as the transversal ferromagnetic resonance frequency. High frequency magnetic field of this wave is parallel to H 0 . For an ordinary wave, having orthogonal polarization with the extraordinary one, effective µ is equal to static permeability which is unit in our case. Substituting µ eff from (4) into Drude formula (2) one can find effective permittivity of FWM.

Figure 1: A wire medium embedded into ferrite host matrix.

Numerical analysis and discussion

First consider FWM as a PBG structure. Numerical analysis was implemented using Green's function in the spectral domain, derived in [6] and adopted to the 2D case, and Galerkin's method.

For the waves, propagating in the xy plane we use the following representation of electric field, induced by current in wires:

E z (x, y) = η ikǫ h µ h L G(x, y, x ′ , y ′ )J z (x ′ , y ′ ) dl ′ ( 5
)
where η is the impedance of free space, L is a contour of a perfectly conducting metal wire, µ h is the effective permeability of the ferrite, given by formula (4) and the scalar Green function has the form

G(x, y, z|x ′ , y ′ , z ′ ) = -1 d 2 ∞ m=-∞ ∞ n=-∞ e -i[αm(x-x ′ )+βn(y-y ′ )] k 2 ǫ h µ h -α 2 m -β 2 n ( 6
)
where α m = k x + q m , β n = k y + q n , q s = 2πs/d. The values ψ x = k x d and ψ y = k y d have sense of phase shifts per lattice period d. Green's function ( 6) can be transformed to the form [6]:

G(x, ȳ) = 1 2d ∞ m=-∞ e -iαm
x tm e tm ȳ e -iψy -tmd

1-e -(iψy +tmd) + e -tm ȳ e +iψy -tmd

1-e +iψy -tmd + e -tm|ȳ| ,

where

x = x -x ′ , ȳ = y -y ′ , t m = α 2 m + β 2 n -k 2 ǫ h µ h . The term G 1 (x, ȳ) = 1 2d ∞ m=-∞ e -iαm xe -tm|ȳ| t m (8)
in (7) contains the logarithmic singularity. The series, asymptotic to (8), is

G1 (x, ȳ) = e -iα 0 x 2π ∞ m=1 cos ms m e -mq = = e -α 0 x 4π {q -ln [2(cosh q -cos s)]} (9)
where s = 2π d

x, q = 2π d |y| and

G 1 (x, ȳ) → G 0 (x, ȳ) ≡ - e -α 0 x 4π ln [g 2 (x 2 + ȳ2 )], (10)
when x → 0, ȳ → 0.

Equating the left-hand side of ( 5) to zero at the surface of the conductor, one obtains the integral equation (IE):

L G(x, y, x ′ , y ′ )J z (x ′ , y ′ ) dl ′ = 0, x ∈ L, y ∈ L. ( 11
)
This first kind Fredholm's IE can be identically transformed to the second-kind IE:

J(x, y) L G 0 (x, y|x ′ , y ′ )dx ′ dy ′ + L [G(x, y|x ′ , y ′ )J(x ′ , y ′ ) -J(x, y)G 0 (x, y|x ′ , y ′ )] dx ′ dy ′ = 0, x ∈ L, y ∈ L. (12)
The integrand in square brackets of (12) does not contain singularity, so this IE can be solved replacing the integral by a quadrature sum. Turning to the cylindrical coordinate system

x = r cos ϕ x ′ = r ′ cos ϕ ′ y = r sin ϕ y ′ = r ′ sin ϕ ′ , (13)
and using the rectangle summation rule, the contour integral can be expressed as

2π 0 f (r, ϕ)dϕ ≈ 2π N N -1 j=1 f (r cos ϕ j , r sin ϕ j ) (14)
where ϕ j = 2π(j -1)/N and N is the number of azimuth nodes. The contour integral over G 0 (x, y|x ′ , y ′ ) can be evaluated analytically:

I 0 = 2π 0 G 0 (r, ϕ | r, ϕ ′ )dϕ = -r ln (rg) - r 4 ln 2. ( 15
)
Replacing integrals (12) by sums (14), one obtains the linear system of algebraic equations:

B k J k + N -1 j=1, j =k A kj J j = 0, k = 1, ..N -1, (16)
where unknowns are the values of the current at the nods ϕ j ,

A kj = 2π N [G(r, ϕ k | r, ϕ j ) -δ kj G 0 (r, ϕ k | r, ϕ j )] , (17)
B k = I 0 . ( 18
)
Dispersion diagram for a square array of metal cylinders inside YIG is shown in Fig. 2. Following parameters are chosen: M 0 = 1.75 KOe, H 0 = 1.8 KOe, ∆ H = 5 Oe, ε h = 15 (relative permittivity of YIG), period of WM lattice d = 5 mm. In opposite to metal PBG structure with isotropic host matrix, dispersion of FWM structure is considerably determined by ferromagnetic resonant properties of the ferrite. It follows from Eq. 4, that there are two special frequency points, determining behavior of µ eff , namely, f ⊥ = ω ⊥ /2π = 7.078 GHz and f b = (ω m + ω 0 )/2π = 9.94 GHz. At the frequency f ⊥ µ eff → ∞ and

µ eff = 0 if f = f b .
Four pass bands are shown in Fig. 2. Grey area marks the frequency range, occupied by a countable spectrum of modes with the accumulation point f ⊥ , where µ eff → ∞. The stop band between f ⊥ and f b is caused by the negative value of µ eff . The value µ eff is positive if the frequency satisfies F > f b but it grows slowly that causes large distances between pass bands in this frequency range. Figure 3 illustrates the applicability of the Drude formula for the description of FWM in the region of positive permittivity. Here ε d is the permittivity calculated using Eq. 4. More accurate ε eff is found from the numerical calculation of the propagation constant for the first-order pass band and is expressed as

ε eff = ψ 2 x d 2 µ eff . ( 19
)
One can see that even for a small ratio r/d = 0.02 the quite good agreement is observed only for small values ε eff and ε d . For the ratio larger than r/d = 0.06 the Drude formula can be used only for qualitative estimations of the effective permittivity. Note, that this consideration refers to the range of the small effective permeability µ ef f ∼ 2.

Surface waves at an interface air-FWM

In this section we will consider several geometries which allow relatively simple description in the framework of effective medium theory.

Geometry 1

The first geometry is shown in Fig. 4. It follows from Maxwell's equations that for this geometry the waves can be split into ordinary wave with and extraordinary wave with

E x = 0, H y = 0, H z = 0, E y = 0, E z = 0, H x = 0
E x = 0, H y = 0, H z = 0, E y = E z = H x = 0.
Fields of the ordinary wave do not interact with both wires and magnetic system. For the extraordinary wave the known dispersion equation for surface waves at the interface air-tangentially magnetized ferrite [7] can be used:

µµ ⊥ β ± γµ a + αµ = 0, ( 20
)
where γ is the propagation constant of the surface wave,

α = γ 2 -k 2 ǫ eff µ ⊥ , β = γ 2 -k 2 , ǫ eff = ǫ xx = ǫ f 1 - k 2 p k 2 ǫ f µ ⊥ .
Thus WM in this geometry can be described in the framework of the effective medium theory where the spatial dispersion does not take place. Dispersion diagram is shown in Fig. 5. Propagation of the surface wave demonstrates the non-reciprocity which is caused by a shift of the maximum of the high-frequency field in the direction, orthogonal to the magnetization. There is the analogy with waves in a tangentially magnetized ferrite slab, located between air on one side and metal on the other side. The surface wave, propagating along the interface air-ferrite has the limiting frequency point f air = F 0 + F m /2 and the wave propagating at the metal interface has f met = F 0 + F m . Dispersion curve calculated at H 0 = 2.9 kOe illustrates the tunability of the structure.

Geometry 2. Inverted structure

In this subsection the wave propagation in an inverted structure which is a ferrite adjacent to WM (see Fig. 6), is discussed. Here we observe that the reduction of thickness of the ferrite slab shifts up dispersion curves. In opposite to the previous case, the dispersion is negative at small γ. Note, areas with negative dispersion are usual in finite-thickness ferrite films [8]. The curves 1,2 and 3 correspond to the wave, whose field is localized at the interface ferrite-WM. In the limit γ → ∞ the cutoff frequency tends to f air , however the effective medium model fully fails at large wavenumbers γ. For the wave, propagating in the opposite direction, whose field is localized at the metal, we have shown the low-γ part of the dispersion curve. At larger γ the dispersion is the same as for the surface wave propagating along the metal boundary of a semi-infinite ferrite.

Geometry 3

Let us consider the next geometry, see Fig. 8. Maxwell's equations for the fields in this structure are split also into two subsystems, corresponding to the TE wave with

E x = 0, H y = 0, H z = 0; E y = E z = H x = 0
H x = 0, E y = 0, E z = 0; H y = H z = E x = 0.
The TE wave is usual surface Damon-Eshbach wave [9] and wires do not affect this wave. Regarding the other polarization, propagation of two waves, TM and TEM is possible [START_REF] Belov | [END_REF]. They are coupled at the interface and additional boundary conditions are necessary for the solution of boundary-value problems [10]. There are no purely surface waves because TEM mode will carry away energy from the interface and it is reasonable only to consider waves in a finite-thickness slab. Waves in this structure, see Fig. 9, are not split into subsystems and all field components are present. Homogenization is needed.

Geometry 4

Conclusion

Novel metamaterial, composed of the magnetized ferrite and wire medium, is described. It can be used as a low-loss tunable metamaterial in microwave range. Analysis made for an unbounded FWM is valid for a finite thickness slab placed between two metal planes (patch and ground plane) if the field does not depend on the z coordinate. It was demonstrated, that for high values of the permittivity of the host medium the Drude formula gives essential error for the calculation of ǫ eff .

Figure 2 :

 2 Figure 2: Dispersion diagram, calculated for r = 0.1 mm (solid curves) and r = 0.3 mm (dashed curves).

Figure 3 :

 3 Figure 3: Effective µ (dotted curve), ε ef f , calculated using Eq. 19 (solid curves) and ε d (dashed curves).

Figure 4 :

 4 Figure 4: Geometry 1. Interface air -wire medium, inserted into the ferrite, magnetized along the wires.

Figure 5 :

 5 Figure 5: Dispersion of surface waves propagating in opposite directions along the interface.

Figure 6 :

 6 Figure 6: Geometry 2. Ferrite slab adjacent to the metal plane (on top) and the WM in air matrix (on the bottom)

Figure 7 :

 7 Figure 7: Dispersion diagram for surface waves propagating at the interface ferrite-WM. The curve number 1 corresponds to the semi-infinite ferrite, the curve 2 corresponds to the slab with the thickness d = 0.3 cm, the curve 3 -d = 0.01 cm. The curve 4 is the low-wavenumber part of the dispersion diagram for the surface wave propagating in the opposite direction.

Figure 8 :

 8 Figure 8: Geometry 3. Wires are normal to the interface.

Figure 9 :

 9 Figure 9: Geometry 4.

Acknowledgment

The authors would like to thank the CNRS (Centre National de la Recherche Scientifique) for its financial support.


Pascal Picart
email: pascal.picart@univ-lemans.fr

Julien Leval
email: julien.leval@univ-lemans.fr

Francis Piquet
email: francis.piquet@univ-lemans.fr

Jean Pierre Boileau
email: jean-pierre.boileau@univ-lemans.fr

Jean Pierre Dalmont
email: jean-pierre.dalmont@univ-lemans.fr

Analysis of Clarinet Reed Oscillations With Digital Fresnel Holography

Keywords: digital holography, holographic interferometry, interferometry, instrumentation measurement and metrology, optical inspection, acoustics, musical acoustics PACS numbers: 07.60.Ly, 42.40.-I, 42.40.Kw 1

This paper describes optical instrumentation devoted to vibration analysis.

Two strategies based on digital Fresnel holography are presented. The first one, called time-averaging consists in the numerical reconstruction of the hologram after recording with an exposure time much greater than the vibration period. Thus, the holograms have their amplitude modulated by a Bessel function. The last strategy, using a pulsed regime, allows the reconstruction of the full movement of the vibration, even if it exhibits very high amplitude. Experimental results presented in this paper concern the study of the vibrations of a clarinet reed under free and forced oscillation regimes.

INTRODUCTION

Digital holography became properly available since its confirmation was established in 1994 [1]. In few years, a lot of wonderful applications have been demonstrated such as digital holographic microscopic imaging [2], securing information [3], surface shape measurements [4], remote metrology [5] and vibrations analysis [6][7][8][9]. In vibration analysis, time averaging [7] was demonstrated and it appears that it is a very useful tool to investigate the properties of objects under sinusoidal excitation. However, the vibration is observed in the sense of its amplitude since its phase is lost during the process. So, pseudo-pulsed digital holography was developed for reconstructing both amplitude and phase of the vibration [8]. An extension of this principle into simultaneous 2D full field vibration analysis by combining digital holography and spatial multiplexing of holograms was also proposed [9]. However such a strategy is only valid for sinusoidal vibrations.

Pulsed lasers are also used for reconstructing vibrations [6]. This allows freezing the object at the time of the recording. Generally, experimental results are limited to double pulse experiments, that means vibration is not fully reconstructed in the sense of its amplitude and phase. Considering that, this paper aims to present a thorough investigation of vibrations by using digital Fresnel holography. The investigation is applied to a clarinet reed. The reed is excited following two modes: a pure monochromatic sinusoidal regime and a pure autooscillation regime. Two interferometers were constructed allowing studies with the time averaging principle and the full reconstruction of the vibration using a pulsed light beam. The paper is organized as follows: Sections 2 describes the basic principle of digital Fresnel holography, in Section 3 different strategies for vibration analysis will be discussed, Section 4 presents the experimental setups developed for the study of forced and free oscillations of a clarinet reed. Section 5 draws the conclusions of the study.

THEORETICAL BASICS

Consider a reference set of coordinate {x,y,z} attached to the surface of an object of interest, z being perpendicular to the surface and considered as the propagation direction of the diffracted light beam (see figure 1 for example). The object surface illuminated by a coherent beam of wavelength λ produces an object wave front noted A(x,y) = A 0 (x,y)exp[jψ 0 (x,y)]. Phase ψ 0 is random because of the roughness of the surface and it will be considered that it is uniform over [-π,+π]. In what follows, it is considered that the object wave propagates through a distance d 0 at which it interferes with a plane and smooth reference wave having spatial frequencies noted {u R ,v R }. At distance d 0 , the reference set of coordinate is chosen to be {x',y',z}. The diffracted field produced by the temporally varying object is given in the Fresnel approximations by [10,

( ) dxdy yy xx d i y x d i t y x A y x d i d d i i t d y x O       + -       + ×       + = ∫ ∫ ∞ + ∞ - ∞ + ∞ - ' ' 2 exp exp , , ' ' exp / 2 exp , , ' , ' 0 2 2 0 2 2 0 0 0 0 λ π λ π λ π λ λ π . ( 11] ( ) ( ) ( ) ( ) ( )
) 1
The reference wave can be simply described using a spatially biased phase:

( ) [ ] ' ' 2 exp ) ' , ' ( y v x u j a y x R R R R + = π . (2)
In the interference plane, the temporally varying hologram H is written:

( ) ( ) ( ) ( ) ( ) ( ) ( ) t d y x O y x R t d y x O y x R y x R t d y x O t d y x H , , ' , ' ' , ' , , ' , ' ' , ' ' , ' , , ' , ' , , ' , ' 0 0 2 2 0 0 * * + + + = . (3)
In digital Fresnel holography, the interferogram is recorded with a matrix of pixels thus discretizing the spatial coordinate in the recording plane. Generally, the detector includes M×N pixels of pitches p x and p y , each of them sized ∆ x ×∆ y . Furthermore, in the case of a time-varying hologram, the hologram effectively recorded by the sensor during the exposure time T is given by:

( ) ( ) ∫ + = T t t dt t d y x H T d y x H 1 1 , , ' , ' , , ' , ' 0 0 . ( 4
)
The object wave can be digitally reconstructed by using a discrete version of the Fresnel integral (equation ( 1)). The numerically reconstructed object will be in focus for a reconstruction distance equal to d R =-d 0 in the reconstruction algorithm. Thus, the diffracted field in the +1 order, noted 1 + R A , at distance d R =-d 0 from the recording plane can be computed with (K,L) ≥ (N,M) data points by evaluating the following equation [9] :

( ) ( ) ( ) ( ) ( ) ( )       + + ×       + - ×       + - - - = - ∑ ∑ - = = - = = + Y lp X kp d j p l p k d j T d lp kp H Y X d j d d j j d Y X A y x K k k y x L l l y x R 0 1 0 2 2 2 2 0 1 0 0 2 2 0 0 0 0 1 2 exp exp , , , exp / 2 exp , , λ π λ π λ π λ λ π . ( 5
)
Note that in the case where (K,L) ≥ (N,M) the digital hologram will be padded with (K-N,L-M) zeros [START_REF] Ferraro | [END_REF]. In reference [7] was demonstrated that the diffracted field in +1 order is given by (* means convolution):

( ) ( ) [ ] ( ) ( ) ( ) ( ) 0 0 0 2 2 0 2 0 2 0 1 , , , , , exp , , d v Y d u X d Y X W Y X A Y X R v u d j d d Y X A R R NM T R R R λ λ δ λ π λ - - * - * × + - = - * + . ( 6
)
where

( ) 0 , , ~d Y X W NM - is the filtering function of bi-dimensional discrete Fourier transform, given by [7] ( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 0 0 0 0 0 0 / sin / sin / sin / sin 1 1 exp , , ~d yp d Myp d xp d Nxp d yp M j d xp N j d y x W y y x x y x NM λ π λ π λ π λ π λ π λ π ×       - - - - = - . ( 7
)
Note that A T (X,Y) is the real object temporally integrated by the detector. Two regimes may be considered for the recording: time averaging regime and pulsed regime. If the hologram is recorded in time averaging regime when the object is excited with a pure sinusoidal excitation (frequency f 0 , T >> 1/f 0 ), the reconstructed object appears to be the real object whose amplitude is modulated by the Bessel function J 0 according to [7]:

( ) ( ) ( ) ( ) Y X A Y X J Y X A m T , , , 0 0 ϕ ∆ = , ( 8
)
where ∆ϕ m is the maximum amplitude of the sinusoidal optical phase variation due to the oscillation of the object. The Bessel modulation indicates the modal structure of the object under the sinusoidal excitation.

In the opposite case, ie the hologram is recorded with a pulse light shorter than the oscillation period of the object, then one gets H(x',y',d 0 ,T) ≈ H(x',y',d 0 ,t j ) with t j being the time at which is started the recording, then we have [8] ( ) ( )

( ) [ ] j T t Y X i Y X A Y X A , , exp , , 0 0 ψ = . ( 9
)
In this last case, the variations of the optical phase ψ 0 between several recordings allow to reconstruct the full movement of the excited object, knowing the sensitivity vector of the interferometric setup. The next section discusses the investigation of vibrations using digital Fresnel holography.

INVESTIGATION OF VIBRATIONS

For the analysis of vibrations, one must distinguish between free and forced oscillations.

In forced oscillation regime, the amplitude of vibrations can be controlled by the excitation source whereas in the free oscillation regime the object imposes its own rule on the observer. So, the digital holographic setup must be configured according to the amplitude of the oscillation. In the time-averaging regime one particularity of the setup is that the vibration is under a monochromatic excitation. This means that the maximum amplitude of oscillation to be measured can be controlled and chosen so that the number of fringes due to the Bessel modulation in reconstructed holograms is low (equation ( 8)).

This way for analysing vibrations is quite simple and doesn't need a complicated setup although the experimental results are partial because of the disappearance of the vibration phase in the reconstructed object [7]. Nevertheless, as a general rule, acoustics engineers are often confronted with free oscillations which are multi-chromatic and have high amplitude regarding the optical wavelength. This kind of vibration is generally not controlled since it is often initiated with nonlinear phenomena leading to an autooscillation phenomenon. The study of such oscillation needs a more elaborated strategy than that used for forced oscillations. Indeed, the oscillation must be finely temporally sampled with a temporal resolution allowing the fringe resolution in the phase data, especially is the amplitude is high. The full movement of the object can then be reconstructed by adding all the partial results obtained with the temporal sampling. Note that the greater the oscillation amplitude is, greater must be the number of recorded

holograms. An example of such vibrations is found in a single-reed woodwind musical instrument. The clarinet is typically such an instrument. Its reed is a small slice of natural cane or synthetic material put into vibration when the musician is blowing air into the instrument. The clarinet is usually considered as the association of a linear resonator, the pipe, and a nonlinear exciter, the reed, subject to the air flow from the mouth. The role of the lowest reed resonance frequency in sound production by single-reed instruments was established in the past [13]. The understanding of the physical mechanisms that lead to the production of sound in the clarinet needs an in-depth investigation of the reed oscillation under monochromatic excitation and under the free oscillation regime. So the experimental investigation of the clarinet reed behaviour was performed with the development of two digital holographic interferometers allowing the visualisation of the reed movement with time averaging holograms and with pulsed holograms.

EXPERIMENTAL SETUPS AND RESULTS

The holographic interferometers are based on Mach-Zehnder architecture with off-axis recordings [7][8][9]. The laser is chosen according to the oscillation regime: a continuous HeNe laser (λ=632.8nm) for time-averaging and a frequency doubled NdYAG pulsed laser (λ=532nm) for the free oscillation regime. The reference wave is issued from a spatial filter in order to be smooth, plane and with uniform amplitude. The sensor is a digital camera with {M×N}={1024×1360} pixels of pitches p x =p y =4.65µm (PCO Pixel Fly). The numerical reconstruction of the holograms is performed according to equation ( 5) with a two-dimensional discrete Fourier transform. The clarinet reed is placed at a distance which depends on the size of the inspected zone. For a rectangular zone of size (a×b), larger size is vertical, a>b, the distance at which must be placed the object is given by [14]:

λ a p d y 4 0 = , (10)
and the couple of optimal spatial frequencies for the reference wave are given by ( )

        ± = 0 , 8 3 , y R R p v u ( 11
)
The optimal values of these parameters allow an optimal recording in which the three orders do not overlap and the object is fully included in the diffracted field.

CASE OF TIME AVERAGING

For the forced oscillation regime the inspected zone is about (a×b)=(35mm×9mm). Thus the recording distance is set to be d 0 =1030mm. In forced oscillation regime, the clarinet reed is submitted to an acoustic wave emitted by a loudspeaker placed behind the reed. In the set-up, the clarinet reed is fixed with a ligature as it is the case in usual conditions.

The amplitude and the frequency of the loudspeaker can be controlled so that the number of fringes in the Bessel modulation remains relatively low, ie the fringes are resolvable.

The laser power is adjusted such that the exposure time of the detector is T=1s, ensuring the condition T>>1/f 0 for f 0 >100Hz. Figure 1 shows the experimental setup for excitation and recording in the time averaging regime. Figure 2 and3 show the vibration modes of the clarinet reed in the frequency range 1kHz-10kHz. It can be observed that the first and the second pure bending mode at 1880Hz and 4500Hz, the first torsion mode at 3260Hz

and other modes which are characteristic of the vibration behaviour of this reed. These modal structures are indicators of mechanical properties of the reed; it can be expected that future work will allow to qualify the aptitude of any reed for the musical playing.

CASE OF FREE OSCILLATIONS

As pointed out previously, the case of free oscillations for the clarinet reed corresponds to playing conditions, i.e. when the clarinet reed is in the autooscillation regime. In order to produce suitable oscillation conditions, an artificial mouth modeling that of a musician playing clarinet was built in the laboratory [15]. This artificial mouth is constituted mainly with an air-proof caisson and an artificial lip modeling the tight-lipped going of the musician. The schematic diagram of the artificial mouth is depicted in figure 4. In the setup, the clarinet reed is fixed with a ligature as it is the case in usual conditions. The artificial lip is made with a rubber membrane filled with water so that it has quite the same consistency as the real musician's lip. The mouth includes an optical window which allows illumination of the reed. To initiate free oscillation of the reed, flow is blown in the mouth. The overpressure in the mouth is of the order of 60 hPa. It induces a flexure of the reed and introduces an air flow into the instrument. Hereafter a threshold pressure, the equilibrium position of the reed becomes unstable and the reed oscillates at a frequency inversely proportional to the effective length of the clarinet. The sound generated by this oscillation is periodic but multi-chromatic and in the setup the fundamental frequency was about 162 Hz. The particularity of the observed autooscillation regime, classically called beating reed regime, is that the amplitude is imposed by the characteristics of the mouthpiece. This amplitude can vary from a few tenths of millimetre to one millimetre. A pressure sensor made with piezo-resistive gages (microphone) is placed in the mouthpiece of the clarinet, as indicated in figure 4 (top view). This sensor measures the pressure fluctuations induced by the mechanical vibration of the reed. The pressure signal delivered by the mouth is of primary use for the optical setup described in figure 5. The laser beam is issued from the double frequency pulsed NdYAG laser with a pulse width of 20ns and pulse energy about 7mJ. The inspected zone is a square region of about (9mm×9mm), thus, the mouth is placed in such a way that the reed is at about d 0 =315mm

from the detector area. In the set-up, illuminating angle θ is set to 22° for experimental convenience. Figure 6 Especially, high amplitude reconstructed with a large number of holograms exhibit less noise than small amplitude reconstructed with a small number of holograms. Note that the signal of figure 7 is not exactly that of the laser sampling since the effective acquisition rate is 5.63Hz. Consequently the reconstructed movement is not the real displacement of the reed but it corresponds to a synthetic displacement constituted with sampling data points temporally spaced from 177.62ms. On figure 7, points A to E correspond to open valve, this means the reed is in its higher position. Points E and F correspond to the closing phase, and then from F to L one observes the phase during which the reed goes to the beak to close the valve. Figures 8 to 13 show the deformation of the reed and the inspected zone included 138490 data points. The color map gives the relative displacement in micron units. Figures 11,[START_REF] Ferraro | [END_REF] exhibit shocks of the reed. These shocks induce a strong deformation of the reed. To the best of our knowledge, these results constitute a first attempt to reconstruct a synthetic reed deformation on a pseudo-period of its autooscillation regime in an artificial mouth.

CONCLUSION

This paper presents methods for vibration analysis using digital Fresnel holography.

Methods are based on time averaging for forced oscillations and pulsed recording devoted to free oscillations. The two methods are applied to a clarinet reed. In the forced oscillation regime, the reed is excited by an acoustic wave with controlled frequency and amplitude such that the reconstructed holograms exhibit resolvable Bessel fringes. These

Bessel fringes are a signature of the mechanical behaviour of the reed. In the case of the autooscillation regime, the reed is placed in an artificial mouth. Deformation of the clarinet reed can be extracted with the recording of 3150 pulsed digital holograms virtually temporally spaced of 2 µs. The method allows high spatial resolution on the reconstructed image. Experimental results show the vibration behaviour of the clarinet reed under autooscillation regime, exhibiting high amplitude shocks during the closed phase. These works give opportunity to acoustics researchers to better understand sound production and musical quality of a clarinet.

6.

Captions

Figure 13

  figures, some shapes (see for example figure12) feature high irregularity while others

Figure 1 .

 1 Figure 1. Digital holographic interferometer for the time averaging regime

Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .Figure 13 .

 2345678910111213 Figure 2. Time averaged holograms of the reed in the frequency range 1kHz-5kHz

ACKNOWLEDGMENT

The authors are very grateful to the reviewers for their careful reading and their useful remarks.

7.


M Moret

S Ruffenach

O Briot

B Gil

Metal Organic Vapour Phase Epitaxial growth of indium-rich InGaN alloys with robust photoluminescence properties

We report the growth of indium-rich InGaN alloys by Metal-Organic Vapour Phase Epitaxy, using ammonia, trimethylgallium and trimethylindium as precursors. Compared to indium nitride, our alloy samples present substantial photoluminescence robustness with temperature. The analysis of the optical properties of these samples versus temperature indicates the existence of two non radiative recombination channels: one with a thermal activation temperature of 77 K and another one with an activation temperature ranging from 300K for InN up to 640K for In 0.72 Ga 0.28 N. The latter competes with the former at low temperatures whilst the former rules the optical properties at ambient conditions.

I. INTRODUCTION

Indium nitride and related heterostructures are promising materials for realizing arsenic-free optoelectronic devices operating in the near infrared portion of the electromagnetic spectrum.

The band gap of indium nitride sits at 650 meV 1 and recent demonstrations using the Molecular Beam Epitaxy technique of InN based quantum well heretostructures with quantum confined stark effect [2][3][4] indicate potential applications of these materials in various field as optical modulators, optical switches and other light emitting /detecting devices. In this paper we report the growth of indium-rich InGaN layers by Metal Organic Vapour Phase Epitaxy with photoluminescence properties robust with temperature as a probe of a fairly reasonable density of non radiative photoluminescence killing centres.

II. SAMPLES GROWTH

The samples were grown by low pressure metal-organic vapour phase epitaxy, using an AIXTRON AIX200/4RF-S reactor. The growth temperature was 550°C and reactor pressure 200 mbar for all samples, with a molar V/III ratio of 15000, in order to prevent any indium droplet formation. Trimethylindium (TMIn) and trimethylgallium (TMGa), provided by SAFC Hitech Limited, were used as group III element precursors, ammonia (NH 3 ) was the nitrogen source. The InGaN layers were grown using N 2 as carrier gas, since hydrogen damages InN, but hydrogen was used as carrier gas for GaN. The purification of H 2 and N 2 , as well as NH 3 was ensured by a complete SAES setup of corresponding purifiers. The different Indium contents were obtained by changing the TMIn/TMGa ratio without changing the total molar flow. The growth rate of InGaN is insensitive to Ga content due to the similar low growth rate of both binaries, which are limited by the low decomposition rate of ammonia. The (In,Ga)N layers thicknesses were measured by SEM and are about 800nm for all the samples. The InGaN layers were deposited on top of a GaN epilayer (typically 1 µmin thickness), deposited onto a sapphire substrate. X ray rocking curve measurements were performed on GaN pseudo-substrate and indicate comparable defects densities in all of them. The samples were characterized by X-Ray diffraction to determine their composition and crystalline quality, the Theta/2theta measurement are displayed Figure 1.

III. PHOTOLUMINESNENCE RESULTS AND MODELISATION

In Figure 2 is reported the evolution of the photoluminescence intensity of an In 0.72 Ga 0.28 N layers for temperature ranging from 4 to 300K. An extended InGaAs Detector was used together with a grating blazed at 1.5 µm. The trend is a collapse of the photoluminescence intensity with increasing temperature, as a probe of the presence of non radiative recombination channels in this sample and all other of ours. The figure shows the robustness of the intensity of the photoluminescence which in case of In 0.72 Ga 0.28 N layers decreases by a factor 30 between 4K and room temperature. This collapse scales two decades for high quality InN films. We have recorded the evolution of the energy of the photoluminescence peaks with temperature. We find a red shift of 59 meV and 1 meV for InN and In 0.82 Ga 0.18 N layers respectively. In contrast to this, we report observation of a 20 meV blue shift for In 0.72 Ga 0.28 N layer. The InN value is comparable with the value reported for best samples in the literature 1 , to the best of our knowledge. Our InN data indicates an initial redshift up to 50K then a weak blue shift and a general trend that is easily fitted using a Varshni equation. This weak S-shape behaviour is an indication of the weakly localized nature of the excitons or e-h pair 5 which recombine at low temperature and a bound to free phase transition near 60K. Concerning the alloys, we are probing the more complex physics of recombination of bound or localized excitons or e-h pairs in strongly disordered systems. It is worthwhile noticing in particular that a simple line shape fitting of the photoluminescence lines performed in order to obtain the maximum of fluorescence energy indicates the onset of blue shift of the photoluminescence line for the 28-percent gallium-rich sample to occurs near 100-120 K ( ~ 10 meV activation energy) simultaneously when the photoluminescence intensity starts to decrease. We attribute this behaviour to simultaneous detrapping of excitons (or photo-created carriers) from shallow low energy potential fluctuations lakes surrounded by barriers of moderate heights and to their further retrapping to other higher energy regions of the crystal surrounded by higher potential barriers. This behaviour is different from the behaviour that is generally reported for gallium-rich InGaN alloys. There are several reasons to this. First, the value of the band gap bowing parameter is large for InGaN alloys (2.8 eV) [START_REF] Yu Davydov | [END_REF] . Fluctuations of the chemical alloy composition give deeper wells for gallium rich alloys than they give for indium-rich ones. The spontaneous production of shallow potential wells by strict chemical disorder is more favourable in case of indium rich alloys than for gallium rich ones. Second, the input parameters used to describe the physics of these materials in the context of quantum mechanics are very different. This reflects the competition between the average spatial extension of the chemical disorder induced potential fluctuations and the spatial extension of carrier (electrons and holes) wavefunction.

Let us restrict ourselves here to electrons, for the sake of our discussion. The effective mass is very small for InN (0.07) while it is substantially high (0.22) for GaN 1 . Let us consider twenty percent alloy disorders near both ends of the full composition range. Effective mass arguments tell us that an electron can probe a more extended crystal region in the case of an indium-rich alloy than it would do it in case of a gallium rich one. Therefore the carrier trappingdetrapping processes are a priori more complicated for indium-rich alloys than they are for gallium-rich ones. This will be inferred later in this paper. We now examine the behaviour of the emitted photoluminescence intensity with increasing temperature. In the most general way, the photoluminescence intensity of excitons or e-h pairs is given by the following proportionality relation:

rad I τ τ ∝ ( Equation 1)
, where τ is the photoluminescence decay time, τ rad is the radiative photoluminescence decay (Equation 2)

Statistical physics arguments tell us that τ rad is proportional to T n/2 where n is the dimensionality of the crystal. Integer n =3 in a bulk, 2 for a quantum well, 1 for a quantum wire and 0 for a zero dimensional confining centre (quantum dot or 3D localizing centre). The non radiative recombination channels (we supposed to have several ones) are thermally activated:

∑ τ = τ i i nr ) T ( ) T ( 1 1
.

(Equation 3)

In general, the evolution of photoluminescence intensity is fitted in the context of an Arrhenius plot, that is to say by plotting the logarithm of the intensity versus the reciprocal temperature in order to obtain in a more or less readable manner, from the intensity behaviour at high temperatures an activation temperature Ti that rules the non radiative recombination channel i.

This gives a rough estimate of what happens. A complementary approach consists in adding to

this a fit of the linear dependence of the photoluminescence intensity in order to correctly describe the complex statistics physics in the range of temperatures where radiative and non radiative recombination channels do compete. From the behaviour of the sample with 28 percents gallium we immediately remark the clear need to introduce two non radiative recombination channels in order to fit the experiment. This means that we have to handle the following equation: 4)

T T T T Be Ae ) T ( I 1 0 1 1 + + ∝ (Equation
, where A, B, T 0 and T 1 are fitting parameters.

Table 1 indicates the result of the full fitting procedure. The remarkable result is the constant value we get for T 0 = 77 K which correspond to an activation energy of some 6.6 meV. It is worthwhile noticing that this weak activation energy corresponds to the detrapping energy we always observe in our InN samples. Next we remark that T 1 increases with the alloy composition from 300 K (26 meV) for InN up to 640 K (55 meV) for In 0.72 Ga 0.28 N.

In Figures 3 and4 are reported the results of the fitting procedures in the context of the reciprocal and direct temperatures respectively. A complementary plot which gives the evolution of the ratio: (B/A)e (T1-T0)/T for these samples is shown in Figure 5 where one can clearly see that the robustness of the photoluminescence in InGaN layers is correlated with the high values reached by this ratio at high temperatures: the relative values of A and B are such that the low temperature T 0 activated process does not dominate. Although it is difficult to make some more deterministic physics we have to remark in Table 1 the large values of B and their increase with increasing the alloy composition. This clearly indicates (having in mind that τ rad is essentially a function of material parameters like the effective masses, the dielectric constant, …; that is to say having in mind that τ rad does not vary so much in these samples), that the characteristics non radiative decay time τ 1 decreases when increasing the alloy composition. This is consistent with: i) our X-ray diffraction patterns which are slightly worse for the alloy than it is for the binary compound as an indication of an increasing number of topological defects with increasing the alloy composition and:

ii) the behaviour of the absolute intensities of the photoluminescence features taken at 4K which decrease with increasing the alloy composition. Table I: Summary of the fitting procedure and correlation function between the fit and the data

IV. CONCLUSION

  time. Further introducing the non radiative part τ nr of the decay time in equation 1 we get:

  photoluminescence intensity with temperature for indium-rich InGaN alloys grown by Metal Organic Vapour Phase Epitaxy (as far as we are aware the materials produced by Molecular beam Epitaxy do not exhibit this robustness). This demonstration indicates that this growth technique has reached a sufficient degree of maturity to provide samples with interesting optical properties in a near future. Further and unavoidable improvements of the crystalline quality of these materials will then open the way towards realization of heterostructures based on InN and related alloys, with remarkable and broad fields of application in the near infrared region.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 12345 Figure 1: theta-2theta X-ray diffraction features

TABLE I

 I

		A	To (K)	B	T 1 (K)	Correlation
	InN	4.9 ± 0.5	77 ± 3	93 ± 14	304 ± 16	0.99958
	In 0.82 Ga 0.18 N 2.4 ± 0.3	77 ± 5	66 ± 15	400 ± 30	0.99834
	In 0.72 Ga 0.28 N 0.8 ± 0.1	78 ± 7	156 ± 31	638± 30	0.99801


Lotfi Alloui
email: lotfi.alloui@lgep.supelec.fr

Frédéric Bouillault

Souri Mohamed Mimoune

Numerical Study of the Influence of Flux Creep and of Thermal Effect on Dynamic Behaviour of Magnetic Levitation Systems with a High-Tc Superconductor Using Control Volume Method

Keywords: 47, 11, Df Finite volume methods -84, 71, Ba Superconducting magnets; magnetic levitation devices

This paper displays some simulation results of dynamic responses of the high-Tc superconductors (HTSC)-Permanent magnet (PM) levitation systems taking into account the influence of the flux creep phenomena and of the thermal effect. We focus on the establishment of a three-dimensional numerical code to solve the nonlinear and coupled equations. A new control volume method is proposed for the resolution of the partial derivative equations of the treaded physical phenomena. The influence is comprehensively displayed by comparing the predictions of dynamic responses of such systems in which the thermal effect in the superconductor is and is not taken into account. The electromagnetic and thermal coupling is ensured by an alternate algorithm. The thermal effect highlights the influence of the temperature on the value of the magnetic levitation force, levitation stabilization time and shows that the vibration center of levitated body had drifted downward.

Introduction

When the stable levitation phenomenon of the permanent magnet (PM) over the high temperature superconductor (HTSC) was found, several projects began to explore the potentials of the HTSC applications such as magnetic bearings [START_REF] Masaie | Numerical evaluation of rotational speed degradation in the superconducting magnetic bearing for various superconducting bulk shapes[END_REF], flywheel energy storage systems [START_REF] Hull | Superconducting bearings[END_REF] and linear transportation system [START_REF] Wang | The first man-loading high temperature superconducting maglev test vehicle in the world[END_REF]- [START_REF] Ren | A hybrid maglev vehicle using permanent magnets and high temperature superconductor bulks[END_REF], etc.

The magnetic levitation transportation systems have the advantage of contactless support which offers high speed with low maintenance and few energy consumption. In the conception of these systems it becomes very important to evaluate precisely the electromagnetic force. Generally, these electromagnetic forces are influenced by the heat losses caused by induced current inside the superconductor. In this paper, the results show that heat losses lead to some important advantages and disadvantages.

In the levitation phenomena, the stability is due to the diamagnetic response of the High Tc superconductors and is greatly enhanced by the additional flux pining phenomenon: the magnetic field penetrates the HTSC and creates fluxoid (shielding current density). The magnetic force results from the interaction between the magnetic field and the fluxoid. The motion of these fluxoid leads to flux flow and flux creep phenomena that make the magnetic force time-dependent since the flux influences the dynamic magnetic force's characteristics of the superconductor. Flux creep is observed both in the conventional low Tc superconductors and in the high Tc oxides. But flux creep is mainly present in high superconductors since usually high temperatures and low pinning potentials are involved. This is demonstrated, for instance, by the 'giant' flux creep and by the broadening of the resistive transition in magnetic fields which is a direct consequence of thermally activated flux creep.

Generally, the flux creep becomes an important phenomenon to consider. Yoshida et al. [START_REF] Yoshida | Magnetic field and force analysis of high Tc superconductor with flux flow and creep[END_REF] who measured the dynamic magnetic force of HTSC show that good agreement between numerical results and measurement data are obtained only if flux creep is considered. In the literature, some macroscopic mathematical models were proposed to reflect the flux creep and the flux flow phenomena of superconductivity such as the model proposed by or the model of power law proposed by Rhyner [START_REF] Rhyner | Magnetic properties and ac losses of superconductors with power law current-voltage characteristics[END_REF]. These models lead to good agreements with experimental data [START_REF] Kasal | Simulation of dynamic levitation force taking flux creep into account[END_REF].

To represent the physical phenomena of the magnetic levitation, several works with different numerical methods and different magnetic formulations have been used to evaluate precisely the magnetic force which is the important parameter to consider [START_REF] Yoshida | Magnetic field and force analysis of high Tc superconductor with flux flow and creep[END_REF]- [START_REF] Kasal | Simulation of dynamic levitation force taking flux creep into account[END_REF]. However, these works do not take into consideration the thermal effect which can play an important role. The induced current in the superconductor creates a power density leading to heat transfer appearance and thus to a temperature increase. In practice, the temperature influences the dynamic magnetic force's characteristics of the superconductor. For instance, T. Suzuki and al [START_REF] Jiang | The magnetic levitation performance of YBaCuO bulk at different temperature[END_REF]- [START_REF] Suzuki | Temperature dependency of levitation force and its relaxation in HTS[END_REF], studied experimentally the effects of the temperature on the levitation force and the time relaxation. They found that the levitation force and the time relaxation of the levitation force decreased with increasing temperature. So it is interesting to develop a numerical code which allows dynamics behaviour modeling of magnetic levitation systems to take into account the thermal effect and which can explain the observed phenomena.

So, this paper presents 3D simulation results of dynamic responses of the high-Tc superconductor-PM levitation systems.

After the description of the considered system, the mathematical models based on the A-V formulation and the numerical code which uses control volume method CVM for solving the nonlinear and coupled problem are examined.

Basic formulation

In this section we introduce the main equations of the superconductor HTSC -permanent magnet PM levitation system. In such system, a cylindrical PM moves (vertically or laterally) over a rectangular high Tc superconductor as shown in Fig. 1.

The effect of flux creep in the superconductor is simulated by using a logarithmic dependence between the critical current and the thermodynamic activation energy, together with the Arrhenius law, as suggested in [START_REF] Brandt | Superconductor disks and cylinders in an axial magnetic field. I. Flux penetration and magnetization curves[END_REF]. The resulting relation between the electrical field and the current density inside the superconductor can be written in the form of power law as:

J (T) J E T T n c c J J ) E(J, ) ( = ( 1 )
Jc is the critical current density which may depend of temperature T and Ec is the critical field of the superconductor. On the basis of Anderson's excitation creep model, the exponent n can be identified with Uc/(k B T), where Uc is the critical activation energy and k B is Boltzmann's constant [START_REF] Sugiura | Numerical electromagnetic field analysis of type-II superconductors[END_REF]. The limits cases of n are the ohmic conductors law for (n=1) and the Bean's model for (n→∞). Generally, n varies from 5 to 25 to describe the flux creep behaviour [START_REF] Kasal | Simulation of dynamic levitation force taking flux creep into account[END_REF]. For the used material in our work (YBaCuo), the good value of n is taken equal to 20 at T = 77 K [START_REF] Berger | Etude des phénomènes couples magnétothermiques dans les supraconducteurs à haute température[END_REF]. In this paper we have adopted a linear model to describe the dependence of the critical current on temperature [START_REF] Braeck | Superconducting trapped-field magnets: temperature and field distributions during pulsed-field activation[END_REF]:

0 B c 0 0 0 0 c0 c T k U n with T T n n(T) Tc T 1 Tc T 1 J (T) J = = - - = ( 2 )
n 0 , J c0 are respectively the n-exponent and the critical current density at T=T 0 (T 0 is the liquid nitrogen bath temperature).

In addition to the gravity force, the PM is submitted to an electromagnetic force resulting from its interaction with the superconductor. This force is the opposite of the magnetic one exerted by the PM on the superconductor, ie:

∫ × ∇ × = Ω dΩ ) A ( J F s ( 3 )
J is the shielding current density, Ω represents the volume of the superconductors and A s = A m +A stands for the total vector potential. A is the induced part of A s due the shielding current in the superconductor, and A m is the excitation part of A s due to the presence of PM. For the calculation of A m , a permanent magnet with a magnetic moment m has been modelled by succession of loops with surface current density j = m∧n [START_REF] Fou | Numerical Modelling of the association of magnet and HTS superconductors[END_REF]. By using the Biot-Savart law, A m is calculated in each point of the whole domain.

The A-V formulation (V: is the electric scalar potential) is used to calculate the induced parts of A s given by:

0 V) t T)( , ( σ V) t T)( , ( σ ) (ν 0 = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∇ + ∂ ∂ - ∇ = ∇ + ∂ ∂ + × ∇ × ∇ A E J A E A s ( 4 )
Where ν 0 is the magnetic reluctivity of vacuum used for superconducting and nonconducting parts, and J s is the source current density defined by:

t T σ ∂ ∂ - = m s A ) (E, J ( 5 )
In (4) the magnetic vector potential A can not be unique unless a gauge constraint is used. The most useful one is the Coulomb gauge (∇.A =0) introduced in (4) by adding a penalty term [START_REF] Coulomb | Finite element three dimensional magnetic field computation[END_REF]- [START_REF] Biro | On the use of the magnetic vector potential in the finite element analysis of 3d eddy currents[END_REF]. Then the left hand size of the first equation in (4) becomes a Laplacien:

0 ) A )( (E, . J ) A )( (E, A . A) ( s = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∇ + ∂ ∂ - ∇ = ∇ + ∂ ∂ + ∇ ∇ - × ∇ × ∇ V t T σ V t T σ ν ν 0 0 ( 6 )
σ is the nonlinear equivalent conductivity of the superconductor and is defined using (1) by:

1 1 - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ) ( E ) (E, T n Ec Ec Jc(T) T σ ( 7 )
So in the superconducting region it is necessary to evaluate

V t t ∇ - ∂ ∂ - ∂ ∂ - = A A E m ( 8 )
In order to obtain the distribution of the temperature in the superconductor, we should resolve the thermal equation in the superconductor:

W T) T) .(κ t T (T) ρ(T)C

p = ∇ ∇ - ∂ ∂ ( ( 9 )
where κ(T) is the thermal conductivity, ρ(T)Cp(T) is the heat capacity per volume, and W = σ(E,T).E 2 is the power magnetic density generated by the electromagnetic phenomenon in the superconductor's sample. Convective boundary conditions are used:

) ( n . ) ( o T T h T T - = ∇ κ - ( 1 0 )
which expresses the thermal flux exchange between the superconductor and the cryogenic fluid. In [START_REF] Suzuki | Temperature dependency of levitation force and its relaxation in HTS[END_REF], h is the convective heat transfer coefficient.

The superconducting material used in this application is the YBaCuO. The thermal and electrical characteristics of YBaCuO are given in table I extracted from literature [13-17-18-19].

To determine the position of the PM, we have to study the movement of the levitated body in vertical (z) and lateral (x) directions. Using the second Newton's law, the equation of movement of the PM when the air damping is neglected, can be expressed as:

⎪ ⎩ ⎪ ⎨ ⎧ = = = - ⇒ = ∑ y m F x m F z m mg F m y x z
Taking into account initial conditions:

) , , ( ) , , ) , , ( ) , , ( :
0 0 0 0 0 0 ( and 0 z y x z y x z y x z y x t
In the general case, where the magnetic permeability takes different values, this basic formulation is not well appropriated for the electromagnetic equation ( 4) which obeys a conservation principle defined by the rotational law. Then it is necessary to modify this basic formulation to derive the algebraic equations. In order to do this, twelve additional nodes are needed WS, WN, ES, EN, SB, ST, NB, NT, BW, BE, TW and TE, so in total twenty one nodes are required. The system of equations in (6) will be integrated over the control volume v e as follow

∫∫∫ ∫∫∫ ∫∫∫ ∫∫∫ ∫∫∫ = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∇ + ∂ ∂ - ∇ = ∇ + ∂ ∂ + ∇ ∇ - × ∇ × ∇ e e e
v v v v v ν ν ) A ( . J ) A ( A) . . ( A)
. ( s [START_REF] Berger | Etude des phénomènes couples magnétothermiques dans les supraconducteurs à haute température[END_REF] System [START_REF] Berger | Etude des phénomènes couples magnétothermiques dans les supraconducteurs à haute température[END_REF] gives differential terms of second order. The terms of the two identical first order derivative lead to the classical basic control volume scheme which handles with free divergence laws (figure 2a), for example:

∫∫∫ Δ Δ ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ e w n s t b e w y z y x z d y d x d x x y A A ν ν ( 14
)
An approximation of first order is used to evaluate numerically expression [START_REF] Fou | Numerical Modelling of the association of magnet and HTS superconductors[END_REF]:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ - ∂ ∂ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂
W P P E y y y ν ν ν A A A A A A A ( 1 5 )
After arranging, (15) leads to the algebraic relationship written as:

P y W y E y A ) ( A A w e w e c c c c + - + ( 1 6 )
The terms of the second order crossed derivatives lead to the new control volume scheme (figure 2b), for example:

∫∫∫ ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ e w n s t b n s ΔxΔz x ν dxdydz x ν y y y A A ( 1 7 )
Using this new configuration of the control volume, the magnetic potential in the equation ( 14) will be calculated in the following way:

) A A A (A ) A A A (A A A A A A y W y SW y E y SE y W y NW y E y NE y sw y se y nw y ne y n s - - + Δ - - - + Δ = - - - = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ x 4 ν x 4 ν Δx ν Δx ν x ν s n s n ( 1 8 )
After arranging, [START_REF] Lkebe | Anisotropic thermal diffusivity and conductivity of YBCO(123) and YBCO(211) mixed crystals[END_REF] leads to an algebraic relationship written as:

) A A A (A ) A A A (A W y SW y E y SE y W y NW y E y NE y - - + - - - + ys yn c c ( 1 9 )
All integrant terms resulting from system (13) are evaluated. Therefore, all the coefficients c i are obtained and describe the physical and the geometrical properties of the cellule surrounding each node. An algebraic system of equations is constructed with appropriate boundaries conditions in the whole domain and expressed as: The global control volume matrices [M] and [N] are non symmetric and non linear. According to the definition of conductivity given by ( 7), we know that the effective conductivity σ of the superconductor is nonlinearly dependent on both electric field intensity E and shielding current density J. At each step, the Crank-Nicolson method is employed to integrate the algebraic system [START_REF] Murakami | Melt processed high-temperature superconductors[END_REF] which is solved iteratively until convergence is reached.

[ ] [ ] [ ] [ ] [ ] [ ] [ ] ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ 0 J V A N 0 N N V A 0 M 0 M S 22 12 11 2 1 & & ( 2
In the same way the control volume method is used to solve the algebraic equation of the differential thermal equation [START_REF] Kasal | Simulation of dynamic levitation force taking flux creep into account[END_REF]. The obtained algebraic system when taking into account the boundary conditions [START_REF] Suzuki | Temperature dependency of levitation force and its relaxation in HTS[END_REF] is [Q] represents the column matrix related to power density and heat exchange in the superconductor. The system is also solved iteratively using the Crank-Nicolson method to treat the time derivative term.

[ ][ ] [ ][ ] [ ] Q T K T C = + & ( 2 1
The coupling algorithm of the electromagnetic and thermal equations consists of an iterative algorithm based on the bloc diagonal Gauss-Seidel method. It permits to separate the non symmetric and the non linear magnetic equation from the non symmetric and non linear thermal one.

Numerical results and discussion

The developed numerical code is used to display some tests results of considered levitation systems. On the one hand, the first use of the numerical code allows simulating the behaviour of the sample superconductor without taking into account the thermal effect. On the other hand it permits to validate the code. The goals of the simulation are to obtain a qualitative physical understanding of the stability features and to describe the difference between cooling the superconductor close to the magnet and cooling it at long distance. We can also predict the appearance of lateral restoring forces and calculate what happens when arbitrary displacements of the magnet occur around the superconductor. The levitated system consists of a rectangular-shaped superconductor (36 mm long, 36 mm wide and 10 mm thick) and a cylindrical PM (20 mm diameter, 10 mm thick and 250g weight).

Fig. 3(a,b,c) shows the magnetic induction vector field and the current density distribution of the component J y which arise when the permanent magnet polarized of 1.1 T goes down toward a cool superconductor. The starting position is located at 20 mm from the superconductor, then the PM will be brought closer to the superconductor at a constant speed rate and will be stopped at 2mm (Fig. 3.c). We observe that the magnetic field due to PM penetrates the HTS bulk from the upper side, and is excluded in the lower part of the sample. The variation in time of this field induces a current density in the region penetrated by the field and leads to a repulsion force between the PM and the superconductor (see fig. 3.d). Therefore, we conclude that is not necessary to use a superconductor with large height, indeed the superconductor develops a surface current such as the current penetration depth is sufficient to exclude the field created by the magnet. Over this depth, the contribution of the bottom tends progressively to be negligible.

The levitation force in Fig. 3.d, which is calculated when the superconductor is cooling at long distance to the magnet, shows a hysteresis property. A repulsive force is obtained when the permanent magnet is moved toward the HTSC as shown in upper line in Fig. 3.d. If the magnet is next moved away from the HTSC, the repulsive force decreases and a weak attractive force appears, which is called the fishing effect [START_REF] Sugiura | Numerical electromagnetic field analysis of type-II superconductors[END_REF]- [START_REF] Murakami | Melt processed high-temperature superconductors[END_REF].

Fig. 4(a,b,c) shows the magnetic induction vector field and the current density distribution of the component J y when the PM is first close to the surface of the superconductor, then moves up from the superconductor and finally goes down. In this case a new physical phenomenon arises. In the first step when the PM is close to the surface of the superconductor (Fig. 4.a)

this later is cooled. In this case, the superconductor stays passive and does not develop induced current. This can be explained by the fact that the flux structure due to the magnet is frozen within the sample and the condition Js = ∇×H =0 is fulfilled.

Thus no change occurs until the magnet is shifted. In the second step when the magnet is moved away from the superconductor, a current distribution appears due to the variation in the time of the magnetic field (Fig. 4.b). We can notice that in this case, this current induces an attractive force (see Fig. 4.d). If the movement of the magnet is reversed, the interaction force becomes repulsive because the new induced current flows in the opposite direction to the previous scheme (see Fig. 4.c), there is no hysteresis property. Fig. 4.d shows the actual behaviour of the levitation force for a specific process. a,b,c) shows the magnetic induction vector field and the current density distribution of the component J y that appear when the magnet is laterally displaced. Notice that these studies can be obtained only by using a 3D approach contrary to the previous ones where 2D modeling can lead to approximated results. An attractive X horizontal force appears when the PM is shifted to the right hand side of the superconductor (x direction). This force increases until the edge of the PM becomes close to the edge of the superconductor, after that the force decreases and tends to zero (Fig. 5.d).

The numerical code is examined a second time in order to predict the experimental phenomena which arise when the temperature is increased and thus the force is reduced [START_REF] Jiang | The magnetic levitation performance of YBaCuO bulk at different temperature[END_REF]- [START_REF] Suzuki | Temperature dependency of levitation force and its relaxation in HTS[END_REF]. This effect can be observed for example in the dynamic response of the levitation behaviour [START_REF] Jiang | The magnetic levitation performance of YBaCuO bulk at different temperature[END_REF]- [START_REF] Suzuki | Temperature dependency of levitation force and its relaxation in HTS[END_REF]. First, some basic simulations are made to show if the influence of the superconductor's temperature on the levitation behaviour is notable. The results of theses simulations point out the difference between the displacement responses and the magnetic force of the levitation system when the thermal effect is and is not considered. Fig. 6 and fig. 7 display the dynamic response of the displacement and the magnetic force of a free levitated PM, where the initial conditions are taken as z 0 = 8 mm and 0 z & =0, respectively.

In fig. 6, it is found that there is a noticeable downward drift of the vibration centre in the response. The thermal effect does not appear initially. But from the instant t = 0.22S a shift appears between the dynamics obtained from the model which takes into account the thermal effect and the one which does not. We can conclude that the thermal effect can not be neglected. By comparing the dynamics of the levitation system obtained we can see that the thermal effect accelerates the stability of the levitation body. Fig. 7 displays the magnetic force varying with time. From it, one sees that when the thermal effect is taken into account, the magnetic force decays faster than when the thermal effect is not considered. This can be explained by the fact that the thermal effect in the high Tc superconductor due to an increasing temperature leads to a decreasing current density.

Figure 8 shows the evolution of the temperature as a function of time along Z axis including the center of the superconductor. Figure 9 shows the temperature profile in superconductor at different instants t = 0.5s, t = 1 and t = 1.575 s, respectively. From these figures, we can see that the temperature increases from the top (near to the power source) to the bottom of the superconductor.

Fig. 10 shows the dependence of the final stable position (static equilibrium position) of the levitation system with the initial position z 0 . This result is very important in the levitation systems application. The value of the final stable position of the PM depends strongly on the z 0 value. The final stable position becomes higher when the initial position is elevated.

Therefore the initial position z 0 plays an important role in the levitated body.

Let z ths and z s be the final stable positions of the levitated PM when the thermal effect is and is not considered in the numerical simulation, respectively. Denote their absolute error by ε = ⎢z ths -z s ⎢. It is evident that the error ε characterizes the influence of the thermal effect in the superconductors on the superconducting levitation systems. The prediction exhibits that the absolute error is mainly dependent on the initial position z 0 for a free levitated PM. In this case the characteristic curve of the absolute error varying with the initial position z 0 is plotted in figure 11. From it, we see that the absolute error increases where the initial position z 0 of the PM increases. When the initial position is fewer than 4 mm the absolute error becomes small and the influence of the thermal effect can be neglected.

Conclusions

In this paper, the shielding current and the electromagnetic field in the superconductor are determined by solving the A-V magnetic equation coupled to the thermal one in order to show the influence of the temperature on the dynamic behaviour of the levitation systems. The flux creep phenomenon is described by using the power law relationship. A three-dimensional numerical code is employed to simulate the dynamic responses of the high-Tc superconductor-PM levitation system. The code is based on a new control volume method for the spatial scheme and a Cranck-Nickolson algorithm for the time one. An iterative solver is used for the resolution of the final algebraic coupled equations. The numerical results exhibit that the drift of the vibration center of the levitation systems is mainly dependent upon the thermal effect when the initial position of the PM above a superconductor is large enough. Hence, the influence of the thermal effect in the superconductor model which leads to larger drift of vibration center, has to be considered in the dynamic simulation of a high superconductor PM levitation system in general case.

Figure captions

  Figure (2.b) presents a projection of a new control volume according to plan X-Y.

  of the magnetic vector potential and the electrical scalar potential; matrix of differentiation of A and V with respect to time variable and matrix of the permanent magnet excitation related to the last term of (6).

  ) Here, [T] denotes the column matrix of the nodal values of the temperature; [ ] T & the column matrix of differentiation of T with respect to time variable and [C] is a nonlinear matrix related to the column matrix [ ] T & and [K] the linear matrix related to [T].

Fig. 5 .

 5 Fig.5.(a,b,c) shows the magnetic induction vector field and the current density distribution of the component J y that

Fig. 1 .Fig. 3 .

 13 Fig. 1. A high Tc-superconducting-PM levitation system Fig. 2. a) Basic control volume scheme, b) New control volume scheme Fig. 3. Current density distribution of the component J y in an HTSC puck and magnetic field induction of an HTSC-magnet system. The figures correspond to different gaps as the magnet is brought close to the HTSC and then brought away. The arrows indicate the direction of the movement of the magnet; a) gap = 20mm; b) gap = 2 mm, c) gap = 20 mm; d) forcedisplacement plot. Minimum gap = 2 mm.

Fig. 4 .

 4 Fig.4. Current density distribution of the component J y in an HTSC puck and magnetic field induction of an HTSC-magnet system. The figures correspond to different gaps as the magnet is brought away to the HTSC and then brought close. The arrows indicate the direction of the movement of the magnet; a) gap = 2mm; b) gap = 21 mm; c) gap = 4 mm; d) forcedisplacement plot. A positive sign corresponds to repulsion, and a negative sign means attraction.

Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .

 5678 Fig.5. Current density distribution of the component J y in an HTSC puck and magnetic field induction of an HTSC-magnet system correspond to different lateral displacement. The arrows indicate the direction of the movement of the magnet. (a) lateral displacement = 2.1mm; (b) lateral displacement = 5.1 mm, (c) lateral displacement = 7 mm, and (d) represent horizontal force-displacement plot. A negative sign corresponds to attraction.Fig.6. Simulation of response of the levitated PM of the levitation system when the thermal effect is and is not taken into account. Fig.7. Simulation of magnetic force excreted on the levitated PM of the levitation system when the superconductor is dealt with and without thermal effect Fig.8. Evolution of the superconductor's temperature.

Fig. 9 .

 9 Fig.9. Temperature distribution in superconductor at different instants respectively at t = 0.25, t=0.75s and at t= 1.47s.

Fig. 10 .Fig. 11 .Fig

 1011 Fig.10. the dependence of the final stable position (static-equilibrium position) of the freely levitated PM from the superconductor with the initial position Fig.11. absolute error of prediction of response displacement between two cases of the superconductor with and without consideration the thermal effect

Table 1

 1

	Alloui et al
	EPJ Applied Physics

TABLE 1 .

 1 SAMPLE PARAMETERSρC p (T)max thermal capacitance [77K ; 300K] -0.66+0.03T-8.3*10 -5 *T 2 +7.8*10 -8 *T 3 [J.cm -3 .K -1 ][17]-[START_REF] Lkebe | Anisotropic thermal diffusivity and conductivity of YBCO(123) and YBCO(211) mixed crystals[END_REF] ρC p (T)min thermal capacitance [77K ; 300K] -0.35+0.019T-4.1*10 -5 *T 2 +2.6*10 -8 *T 3 [J.cm -3 .K -1 ] [17]-[START_REF] Lkebe | Anisotropic thermal diffusivity and conductivity of YBCO(123) and YBCO(211) mixed crystals[END_REF]

	Symbol	Quantity	value or expression
	κ	thermal conductivity	6 [W.m -1 .K -1 ]	[19]
	Tc	critical temperature	92 [K]	[19]
	Ec	critical electrical field	1 [µV/cm]	[13]
	n 0	n-exponent at 77K	20	[13]
	J c0	critical current density at 77K	440 [A/mm 2 ]	[ 1 9 ]
	h	convective heat transfer coefficient	400 [W/(m 2 .K)]	[13]

& & & & γ F (1 1 )

& & & & & &

Here F z , F x and F y are the component of the electromagnetic force F given by equation [START_REF] Wang | The first man-loading high temperature superconducting maglev test vehicle in the world[END_REF], m is the mass of the permanent magnet and g is the gravitational acceleration. Using equation [START_REF] Brandt | Superconductor disks and cylinders in an axial magnetic field. I. Flux penetration and magnetization curves[END_REF], it's possible to determine the new position of PM due to the variation of the field.

Control Volume Method

The numerical method used in this paper is the control volume method (CVM). This method can be regarded as a special version of the weighted residuals method. It consists in dividing the domain into a number of nonoverlapping subdomains or control volumes such that there is one control volume surrounding each node (Fig. 2.a). The weighting function is setting to be unity over one subdomain at a time and zero everywhere else. This variant of the method of weighted residuals is called the subdomain method or CVM.

The basic formulation of CVM deals with differential equations which seem to obey a conservation principle defined by the divergence law. In this case, six nodes surrounding each principle node 'P' are necessary to discretize the algebraic equation : West 'W', East 'E', South 'S', North 'N', bottom 'B' and Top 'T' (Fig. 2.a). This basic formulation with seven nodes in total is suitable for heat and mass transfer, fluid flow, chemical reaction and others related processes.