Abstract-This paper presents a general method for studying the mechanical properties of a ferrofluid seal by using a three-dimensional analytical approach based on the coulombian model of a magnet. The fundamental Maxwell's equations lead us to define the concept of magnetic energy of the ferrofluid seal by using only the three-dimensional equations of the magnetic field created by ring permanent magnets radially magnetized.

Our study corresponds to the specific case when the ferrofluid is submitted to a very high magnetic field. Under these conditions, the mechanical properties of the ferrofluid depend only on the magnetic field created by the permanent magnets. Throughout this paper, the remanence polarization J of the magnets used is higher than 1T. The magnetic field we use in order to align the magnetic particles is very intense, greater than 400 kAm -1 . Consequently, the magnetic particles are assumed to be saturated and the magnetic field they create can be omitted.

In this paper, a cylindrical structure consisting of two outer ring permanent magnets and an inner non-magnetic cylinder is considered. In addition, a ferrofluid seal is placed between them. The calculation of the magnetic pressure of the ferrofluid seal has been analytically established in three dimensions in order to determine its shape. Moreover, the geometrical evolution of the ferrofluid seal shape is presented when the inner non-magnetic cylinder crushes the ferrofluid seal. The radial stiffness of the ferrofluid seal is determined in three dimensions when the inner cylinder is decentered. Furthermore, a way of obtaining the ferrofluid seal static capacity is discussed.

It must be emphasized here that the three-dimensional method presented in this paper differs from the ones commonly used in the literature because the only term we consider in the equation of equilibrium describing the state of the ferrofluid seal is determined accurately. Such a method has been used for the design of ironless loudspeakers that have demonstrated some promising acoustical tests.

INTRODUCTION

Analytical calculations of the magnetic field produced by ring permanent magnets are useful for the design of devices using permanent magnets and ferrofluid seals. In this paper, we show how these analytical calculations can be used for studying and optimizing ferrofluid seals in ironless structures. To our knowledge, this way of studying a ferrofluid seal appears for the first time in the literature and we think that it can be useful for the design of magnetic bearings using ferrofluid seals or ironless loudspeakers. Indeed, some recent papers dealing with ironless loudspeakers have clearly taken the same assumptions as the ones used in this paper for studying audio speakers containing ferrofluid [START_REF] Rosensweig | Study of audio speakers containing ferrofluid[END_REF]. The first studies dealing with ferrofluids were certainly carried out by Rosensweig [START_REF] Rosensweig | Ferrohydrodynamics[END_REF] and many authors have been interested in the promising applications using ferrofluids [START_REF] Raj | Advances in ferrofluid in ferrofluid technology[END_REF].

On the other hand, some progress in the analytical calculations of the magnetic field produced by ring permanent magnets [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF]- [START_REF] Selvaggi | Calculating the external magnetic field from permanent magnets in permanent-magnet motors -an alternative method[END_REF] have allowed us to determine accurately the magnetic energy of the ferrofluid by using the fundamental's Maxwell equations.

We can say that studies on ferrofluids deal generally with either their chemical and physical properties or their mechanical behavior [START_REF] Choi | Simulation of hydrostatical equilibrium of ferrofluid subject to magneto-static field[END_REF]. Ferrofluids have various engineering applications [START_REF] Shah | Ferrofluid squeeze film in a long bearing[END_REF], such as damping [START_REF] Bajkowski | A model for a magnetorheological damper[END_REF], for medical applications [START_REF] Park | New design of the magnetic fluid linear pump to reduce the discontinuities of the pumping forces[END_REF] or for ironless loudspeakers [START_REF] Lemarquand | Ironless loudspeakers[END_REF]- [START_REF] Berkouk | Analytical calculation of ironless loudspeaker motors[END_REF]. These last applications are the ones we are currently studying.

Moreover, ferrofluids are also often used as squeeze films [START_REF] Tarapov | Movement of a magnetizable fluid in lubricating layer of a cylindrical bearing[END_REF] in bearings or as rotating shaft seals. Since then, numerous studies in the field of ferrofluid dynamic bearings have been carried out. Both static and dynamic characteristics of these bearings have been studied theoretically [START_REF] Walker | Ferrohydrodynamics thrust bearings[END_REF]- [START_REF] Cunha | A new boundary integral formulation to describe three-dimensional motions of interfaces between magnetic fluids[END_REF]. Moreover, recent trends in the ferrofluid lubrication applications are described and discussed, taking into account various phenomena including cavitation [START_REF] Chen | Some design and prototyping issues on a 20 krpm hdd spindle motor with a ferro-fluid bearing system[END_REF]- [START_REF] Meng | An analysis on the magnetic fluid seal capacity[END_REF]. This paper presents a general method for studying the mechanical properties of a ferrofluid seal by using a three-dimensional analytical approach based on the coulombian model of a magnet. The structure considered is composed of two parts that are either made out of permanent magnets (the stator) or non magnetic materials (the moving part). Consequently, no reluctance effect is observed. The moving part has a translation movement, but this study does not take its effects into account. Thus, we only deal with the static properties of the ferrofluid seal. The magnetic field created by the ring permanent magnets is very intense. Therefore, the ferrofluid is largely saturated. Although numerical approaches are commonly used to study ferrofluid seals, analytical approaches with exact formulations allow a great flexibility in optimizing the dimensions of the configurations studied. Consequently, in our specific use of the ferrofluid seal, it is more interesting to use an exact analytical term rather than a lot of simplified terms in the Bernoulli equation for describing the state of the ferrolfuid seal.

It is noted that this study is entirely based on the magnetic pressure calculation which is analytical throughout this paper. Indeed, such an analytical approach allows a precise study of the ferrofluid seal shape when it is crushed, which is often an occurring state during the functioning. The seal capacity is also determined and discussed and a analytical method for obtaining the radial stiffness of the ferrofluid seal is also presented.

We think that this method paves the way for new and promosing applications using both permanent magnets and ferrofluid seals. All the analytical expressions used in this paper have been published in previous papers and are a good alternative to the finite element method for the design of such structures.

STRUCTURE DESCRIPTION

Geometry

The studied structure consists of two outer ring permanent magnets and an inner non-magnetic cylinder. Such a structure is considered in this paper because it allows us to have important gradient values. Moreover, this structure is not very complicated to build since it requires only two ring permanent magnets. As a result, the magnets are stacked and their magnetizations are both radial, but in opposed directions (Fig. 1). The non-magnetic cylinder is centered with the two ring permanent magnets. Moreover, both outer ring permanent

h Z 0 R r in r out Figure 1
. Geometry : two outer ring permanent magnets and an inner non-magnetic cylinder with a ferrofluid seal between them; the ring inner radius is r in , the ring outer radius is r out , the height of a ring permanent magnet is h. magnets have the same dimensions and the same magnetization in magnitude. The ring inner radius is r in , the ring outer radius is r out and the ring permanent magnet height is h. The z axis is a symmetry axis.

Basic equations

The magnetic field produced by the ring permanent magnets can be determined by using a fully analytical approach. The fundamental equations of Maxwell lead us to determine the magnetic energy in the ferrofluid. Indeed, let us consider the four fundamental Maxwell's equations:

∇. B = 0 (1)

∇ ∧ H = j (2) ∇. D = ρ (3) ∇ ∧ E = - ∂ B ∂t ( 4
)
where B is the magnetic induction field, H is the magnetic field, j is the volume current density, D is the electric flux density, E is the electrostatic field and ρ is the electrical charge. In our structure, no currents are considered as the magnetic field is created only by the permanent magnets. The vector fields B and H are defined for all points in space with the following relation:

B = µ 0 H + J (5)
where µ 0 is the permeability of the vacuum and J is the polarization vector of the magnet. When the magnetic field is determined outside of the magnet, we have J = 0. The analogy between the Mawxell's equations leads us to write that:

∇. H = - ∇. J µ 0 = σ * µ 0 (6)
where σ * corresponds to a fictitious magnetic pole density. On the other hand, the magnetic field H verifies the following equation:

∇ ∧ H = 0 (7)
Thus, H can de deducted from a scalar potential φ( r) by

H = -∇ (φ( r)) (8)
Let us consider a structure owing several ring permanent magnets, that is, an assembly of magnetic pole surface densities and magnetic pole volume densities; with ( 6) and ( 7), we can write that

φ( r) = 1 4πµ 0   i S i J k .d S i r -r + j V j -∇. J k r -r dV j   (9)
where J k is the magnetic polarization vector owing to the ring permanent magnet k and r -r is the distance between the observation point and a magnetic charge contribution. We deduct immedialtely that the magnetic field created by the ring permanent magnets is determined as follows: As there is a cylindrical symmetry in our configuration, the magnetic field depends only on two variables r and z. This relation is used in the next section for showing the magnetic field in our configuration.

H = -∇.   1 4πµ 0   i S i J k .d S i r -r + j V j -∇. J k r -r dV j     (10)

Magnetic field created by the two ring permanent magnets

As shown in Fig 1, our geometry is comoosed of two ring permanent magnets. The magnetic field created by the two outer ring permanent magnets depends on both the radius, r, and the altitude, z. The radial component of the magnetic field created by the two outer ring permanent magnets is denoted H r (r, z) and the axial component is denoted H z (r, z). It is noted that the expressions of these components have been established in previous papers [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF]- [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF]. The magnetic field H is represented in Fig. 2. The ring inner radius r in equals 25 mm, the ring outer radius equals 28 mm and the permanent magnet height, h, equals 3 mm. The magnetic pole surface density +σ * equals 1 T. Figure 2 shows that the magnetic field H is the strongest near the permanent magnets, especially where the magnetic field gradient is the most important.

MAGNETIC PRESSURE ON THE FERROFLUID SEAL

Assumptions

It is emphasized that the structure considered in this paper contain no iron. Moreover, the further calculations are presented for magnets with 1 T remanent magnetization in order to normalize the results.

In fact, the magnets used in the prototypes are Neodymium Iron Boron ones for which the remanent magnetization can reach 1.5 T. Furthermore, the interesting regions of space are the ones where the ferrofluid goes. Besides, in the field gradient due to the external magnetic field created by the two outer ring permanent magnets, the magnetic particles are attracted towards the higher intensity regions of the magnetic field. With these 1 T normalized magnets and for the proposed configurations, the magnetic field intensity there is greater than 400 kA/m. Of course, all the field values are proportional to remanent magnetization value. Then, we use commercial ferrofluids either from the company Ferrotec or Ferrolabs. Such ferrofluids have a saturation magnetization, M s , smaller than 32 kA/m and a particle concentration below 5.5 %. It is to be noted that for bearing or loudspeaker applications, a great bearing effect is seeked which requires high saturation magnetizations. Therefore, the magnetic field, H, created by the permanent magnets is far higher than the ferrofluid critical field [START_REF] Matthies | Numerical simulation of normalfield instability in the static and dynamic case[END_REF]. So, the ferrofluid is totally saturated and its magnetization is denoted M s . Another consequence is that the magnetic relative permeability of the ferrofluid equals one. Thereby, the presence of ferrofluid does not modify the magnetic field created by the permanent magnets and the field created by the ferrofluid itself is omitted. Moreover, the ferromagnetic particles are assumed to be small saturated spheres which can be freely oriented in all the directions of space. Besides, this study is a static one. Thus, all the particles of the saturated ferrofluid are aligned with the permanent magnet orienting field. Consequently, the ferrofluid magnetization has the same direction as the orienting field. Furthermore, the aggregation in chains of the ferrofluid particles is omitted [START_REF] Ivanov | Ferrofluid aggregation in chains under the influence of a magnetic field[END_REF]. It is noted that, when the device is at laboratory temperature and at rest, the aggregation phenomenon is observed.

Moreover, the thermal energy E T (E T = kT where k is Boltzmann's constant and T is the absolute temperature in degrees Kelvin) and the gravitational energy E G (E G = ∆ρV gL where V is the volume for a spherical particle, L is the elevation in the gravitational field, g is the standard gravity, ∆ρ is the difference between the ferrofluid density and the outer fluid) are neglected.

Eventually, this paper deals with the ferrofluid free boundary surface. Its shape depends on the result of the force competition at this boundary. On one hand, the magnetic pressure, p m , is exerted. It is given by:

p m (r, z) = µ 0 M s . H(r, z) = µ 0 M s H r (r, z) 2 + H z (r, z) 2 (11)
where both magnetic field components H r (r, z) and H z (r, z) are analytically calculated [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF]. It is to be noted that this expression is given below in a fundamental form. However, this expression can be written in a form more suitable for numerical treatments [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF]. The magnetic pressure p m (r, z) is given by [START_REF] Selvaggi | Calculating the external magnetic field from permanent magnets in permanent-magnet motors -an alternative method[END_REF].

p m (r, z) = µ 0 M s H r (r, z) 2 + H z (r, z) 2 (12)
where M s is the intensity of magnetization of a magnetic particle. The radial component H r (r, z) of the magnetic field created by the two permanent magnets is defined by [START_REF] Choi | Simulation of hydrostatical equilibrium of ferrofluid subject to magneto-static field[END_REF].

H r (r, z) = σ * πµ 0 i(1 + u) (ς(u 1 ) -ς(u 2 )) (13)
with

ς(u) = ξ 1 (-(a 1 d + b 1 (c + e 1 )))F * i sinh -1 [ √ -c+d-e 1 √ c+e 1 +du ], c+d+e 1 c-d+e 1 d √ -c + d -e 1 e 1 d(1+u) c+e 1 +du √ 1 -u 2 + ξ 1 (b 1 c -a 1 d)Π * e 1 c-d+e 1 , i sinh -1 [ √ -c+d+e 1 c+e 1 +du ], c+d+e 1 c-d+e 1 d √ -c + d -e 1 e 1 d(1+u) c+e 1 +du √ 1 -u 2 + ξ 2 (-(a 2 d + b 2 (c + e 2 ))) F * i sinh -1 [ √ -c+d-e 2 √ c+e 2 +du ], c+d+e 2 c-d+e 2 d √ -c + d -e 2 e 2 d(1+u) c+e 2 +du √ 1 -u 2 + ξ 2 (b 2 c -a 2 d)Π * e 2 c-d+e 2 , i sinh -1 [ √ -c+d+e 2 c+e 2 +du ], c+d+e 2 c-d+e 2 d √ -c + d -e 2 e 2 d(1+u) c+e 2 +du √ 1 -u 2 - η 3 ((a 3 d -b 3 e 3 )) F * i sinh -1 [ √ -d-e 3 √ e 3 +du ], -d-e 3 d+e 3 d √ -d -e 3 (-c + e 3 ) d(1+u) e 3 +du √ 1 -u 2 - η 3 (b 3 c -a 3 d)Π * -c+e 3 d+e 3 , i sinh -1 [ √ -d+e 3 e 3 +du ], -d+e 3 d+e 3 d √ -d -e 3 (-c + e 3 ) d(1+u) e 3 +du √ 1 -u 2 - η 4 (a 4 d -b 4 e 4 ) F * i sinh -1 [ √ -d-e 4 √ e 4 +du ], -d+e 4 d+e 4 d √ -d -e 4 (c + e 4 ) d(1+u) e 4 +du √ 1 -u 2 - η 4 (b 4 c -a 4 d)Π * -c+e 4 d+e 4 , i sinh -1 [ √ -d-e 4
The axial component of the magnetic field created by the two ring permanent magnets is given by [START_REF] Bajkowski | A model for a magnetorheological damper[END_REF].

H z (r, z) = σ * πµ 0   -r in K * -4rr in (r-r in ) 2 +z 2 (r -r in ) 2 + z 2   + σ * πµ 0   r in K * -4rr in (r-r in ) 2 +(z-h) 2 (r -r in ) 2 + (z -h) 2   - σ * πµ 0   r in K * -4rr in (r-r in ) 2 +z 2 (r -r in ) 2 + z 2   + σ * πµ 0   r in K * -4rr in (r-r in ) 2 +(z+h) 2 (r -r in ) 2 + (z + h) 2   ( 15
)
ξ i = d(-1 + u) c + e i + du ( 16
)
η i = d(-1 + u) e i + du ( 17
)
where K * [m] is given in terms of the incomplete elliptic integral of the first kind by ( 18)

K * [m] = F * [ π 2 , m] (18)
F * [φ, m] is given in terms of the elliptic integral of the first kind by [START_REF] Berkouk | Analytical calculation of ironless loudspeaker motors[END_REF]:

F * [φ, m] = θ=φ θ=0 1 1 -m sin(θ) 2 dθ ( 19
)
Parameters a 1 r in rz b 1 -r 2 in z c r 2 + r 2 in d -2rr in e 1 z 2 a 2 -r in r(z -h) b 2 r 2 in (z -h) e 2 (z -h) 2 a 3 r in rz b 3 -r 2 in z e 3 r 2 + r 2 in + z 2 a 4 r in r(-z -h) b 4 -r 2 in (-z -h) e 4 r 2 + r 2 in + (z + h) 2
Table 1. Definition of the parameters used in ( 14)

Π * [n, φ, m] is given in terms of the incomplete elliptic integral of the third kind by ( 20)

Π * [n, φ, m] = φ 0 1 (1 -n sin(θ) 2 ) 1 -m sin(θ) 2 dθ ( 20
)
The parameters used in ( 14) are defined in Table 1. Moreover, when we input [START_REF] Choi | Simulation of hydrostatical equilibrium of ferrofluid subject to magneto-static field[END_REF] in Mathematica, we have to take the real part of H r (r, z) because of the noise calculus. What's more, the parameter i used in ( 14) is the imaginary number (i 2 = -1).

As the ferrofluid considered in this paper has a saturation magnetization of 32 kA/m and the magnetic field is greater than 400 kA/m, the value of this pressure is higher than 12800 N/m 2 . On the other hand, the surface tension exists. But when the values of both the surface tension coefficient, A, (A equals 0.0256 kg/s 2 for the used ferrofluids) and the radius of curvature are considered, the effect of the surface tension can be omitted: this latter does not deform the free boundary surface.

Finally, for hydrodynamic pressures which equal zero or have low values, the seal free boundary surface is a magnetic iso-pressure surface. Then, its mechanical properties, such as its capacity or its stiffness can be studied. Further dynamic studies could also be carried out with a perturbation method from this rest position [START_REF] Kuzhir | Free boundary of lubricant film in ferrofluid journal bearings[END_REF]. The magnetic pressure p m (r, z) has been determined with [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF]. Figure 3 shows that the magnetic pressure is greater near the two ring permanent magnets, especially where the magnetic field gradient is the strongest. It can be noted that such a figure is useful to see quickly where the magnetic pressure is the most important. Moreover, Fig.

(3) gives also an indication about how the magnetic pressure varies according to r and z; we see that this variation is very important either to the permanent magnets. Furthermore, this representation is interesting since it shows that the potential energy is concentrated in a very small ferrofluid volume. As a consequence, it gives an indication about what quantity of ferrofluid should be used to design a ferrofluid seal. Indeed, some viscous effects can appear when the inner cylinder moves between the two outer ring permanent magnets. If a large quantity of ferrofluid is used, the ferrofluid seal is thick, the potential energy increases but the viscous effects become an actual drawback according to the dynamic movement of the inner cylinder. If a too small amount of ferrofluid is used, the viscous effects disappear but all the interesting properties of the ferrofluid seal (damping, stability, linearity,...) disappear as well. In short, to a given geometry ( here two ring permanent magnets with an inner non-magnetic cylinder) corresponds to an adequate quantity of ferrofluid which has some interesting physical properties with very little viscous effects.

Potential energy of the ferrofluid seal

The potential energy of the ferrofluid seal is defined by ( 21):

E m = - (Ω)
p m (r, z)dV [START_REF] Walker | Ferrohydrodynamics thrust bearings[END_REF] where (Ω) is the ferrofluid seal volume. Indeed, this energy allows the calculation of the seal mechanical properties. As a remark, the magnetic pressure is given in N/m 2 and the potential energy in J.

STUDY OF THE FERROFLUID SEAL SHAPE

Ferrofluid seal shape under a magnetic field

We first analyze the ferrofluid shape without taking into account the inner non-magnetic cylinder. As we assume each magnetic particle to be a small sphere which can be freely oriented in all the directions of space, each magnetic particle is oriented with the direction of the orienting field in order to minimize its potential magnetic energy. More exactly, the magnetic particles go where the magnetic field is the most important by always being aligned with it. Figure 4 is a two dimensional contour representation of the potential energy of the ferrofluid seal. Such a figure is useful to see the points in space which have the same potential energy. This is interesting because this element of information allows us to study the deformation of the seal when the inner cylinder is radially decentered. The aim of this figure is to show where the ferrofluid goes at first: the warmest areas are the first occupied ones. The shape of the ferrofluid seal between r = 24.6 mm and r = 25 mm can be approximated as half an ellipse. Consequently, when the ferrofluid seal thickness is smaller than 0.4 mm, its contour can be written in terms of an equation of an ellipse [START_REF] Tiperi | Overall characteristics of bearings lubricated ferrofluids[END_REF].

(r -r i ) 2 a 2 i + z 2 b 2 i = 1 (22)
To illustrate [START_REF] Tiperi | Overall characteristics of bearings lubricated ferrofluids[END_REF], Table [START_REF] Rosensweig | Ferrohydrodynamics[END_REF] gives the parameters of the equations of ellipse describing the ferrofluid seal in Fig. 4 for r between 24.6 mm and 25 mm. The approximation error between the equations of ellipse and the real contour shape of the ferrofluid seal is also given. When some ferrofluid is added, the shape of the ferrofluid seal changes and for Table 2. Parameters (in mm) of the equations of an ellipse describing small ferrofluid seals r < 24.6 mm, as represented in Fig. 5, the ferrofluid seal contour cannot be represented in terms of an equation of ellipse. Furthermore, we can point out that if a large amount of ferrofluid is added, the ferrofluid seal size becomes very large because of the magnet edge effects. In our configuration, the thickness of the ferrofluid seal for which the shape changes is 0.6 mm. This element of information is very important since it gives clearly what quantity of ferrofluid sould be used to make a small seal. We have also represented the final three-dimmensional form of the ferrofluid seal with the same dimensions (r in = 25 mm. J = 1 T, h = 3 mm) in Fig 6. . Representation of growing ferrofluid volumes under a magnetic field from the warmest areas to the coolest ones: the warmest colors correspond to the region in space where the magnetic energy is the strongest

Shape of the crushed ferrofluid seal

We take now the non-magnetic cylinder into account. The device has to be dimensioned in such a way so that the cylinder crushes the seal, in order to reach the intended watertightness. Figure 7 shows where the ferrofluid goes when the inner non-magnetic cylinder crushes the ferrofluid seal. The phenomenon is axisymmetrical, as the inner cylinder is supposed here to be centered. As seen in Fig. 4, for small quantities of ferrofluid, the ferrofluid seal contour can be described in terms of an equation of ellipse. Therefore, when a small ferrofluid seal is crushed, the quantity of ferrofluid goes symmetrically under and above the initial ferrofluid seal by describing an equation of ellipse. The change of the ferrofluid seal shape generates an energy reduction. Table 3 presents the energy reduction when the cylinder, being first at a distance in length which equals 0.3 mm from r in , crushes radially the ferrofluid seal. The energy reduction is very important (68%) when the radial displacement of the inner cylinder is 0.2 mm.

CALCULATION OF THE STATIC CAPACITY OF THE FERROFLUID SEAL

The aim of this section is to determine the ferrofluid seal static capacity. Indeed, the pressure underneath the non-magnetic cylinder can be Consequently, a pressure gradient appears and the ferrofluid seal is thus deformed. Therefore, the calculation of the ferrofluid seal capacity is required. For this purpose, a second configuration is considered, which corresponds to the case when a cylindrical air gap appears in the seal along the cylinder because of an applied pressure on one side of the seal. This second configuration is represented in Fig. 8. Two steps are necessary to determine the capacity of the ferrofluid seal. The first step is the determination of the potential energy of the ferrofluid seal without any air gap in the seal. To do so, a numerical integration of ( 21) leads to a first value of the potential energy which is denoted E m (1). The used numerical method is the GaussKronrod method.

A three-dimensional representation of the ferrofluid seal perfored in shown in Fig 9 . The second step is the determination of the potential energy of the ferrofluid seal with the air gap. Again, a numerical integration leads to a second value of the potential energy denoted E m [START_REF] Rosensweig | Ferrohydrodynamics[END_REF]. The energy difference is denoted ∆E m and verifies [START_REF] Raikher | Orientational dynamics in magnetic fluids under strong coupling of external and internal relaxations[END_REF]:

∆E m = E m (1) -E m (2) (23)
The energy difference corresponds to the pressure work δW (P ) and satisfies [START_REF] Miyake | Sliding bearing lubricated with ferromagnetic fluid[END_REF]: ∆E m = δW (P ) = P Sd [START_REF] Miyake | Sliding bearing lubricated with ferromagnetic fluid[END_REF] where S is the surface area of the air gap and d is the thickness of the perforation. Consequently, the capacity P lim verifies [START_REF] Chang | A theoretical and experimental study of ferrofluid lubricated four-pocket journal bearing[END_REF]:

P lim = δW (P ) Sd ( 25
)
A numerical application has been done with several ring inner radii (Fig. 10). We define the ferrofluid thickness as the thickness of the 4. Volume and Magnetic field corresponding to a given ferrofluid thickness ferrofluid between the outer ring permanent magnets and the inner cylinder. Figure 10 shows that the thinner the ferrofluid seal is, the higher the pressure gradient to which the ferrofluid seal can resist is. Table 4 gives both the volume corresponding to each ferrofluid thickness and the smallest magnetic field Hlim in the concerned volume.

RADIAL STIFFNESS OF THE FERROFLUID SEAL

Determination of the radial stiffness

This section presents a way of obtaining the radial stiffness of the ferrofluid seal. For this purpose, a third configuration is considered (Fig. 11), which corresponds to a decentered inner cylinder. Figure 12 illustrates the dependence of the radius with θ. The ferrofluid seal is thus more crushed when the cylinder comes closer to the magnets. The problem here is no longer axisymmetrical, but fully three-dimensional. Indeed, where the radial room for the ferrofluid is narrower, the seal becomes larger along the Oz axis, but the ferrofluid also moves towards the region of increasing room and so, the angle θ interferes. The way of obtaining the radial stiffness is done in two steps. The potential energy E m (1) is first calculated when the nonmagnetic cylinder is centered (Fig. 1) with [START_REF] Walker | Ferrohydrodynamics thrust bearings[END_REF]. Then, the potential energy E m (2) is calculated when the non-magnetic cylinder is not centered (Fig. 11). In the second configuration, the limits of the integrals which allow us to determine the potential energy depend on the angle θ. Consequently, the potential energy E m (2) is determined with [START_REF] Zhang | Static characteristics of magnetized journal bearing lubricated with ferrofluids[END_REF].

E m (2) = 2π 0 r in r 0 +e cos(θ) z 0 -z 0 e m (r, z)rdrdθdz ( 26
)
where r 0 and e are determined by the equation of the circle which is not centered given by [START_REF] Osman | Static and dynamic characteristics of magnetized journal bearings lubricated with ferrofluid[END_REF].

r(θ) = r 0 + e cos(θ) (27)
where r 0 = 24.7 mm and e = 0.1 mm. Figure 12 shows how the inner cylinder is decentered. Consequently, the ferrofluid shape is not axisymmetrical since the inner cylinder is radially decentered.

As the ferrofluid shape is not axisymmetrical, the limits of the integrals in Eq. 26 depend on the angle θ. Thanks to the determination of E m (1) and E m (2), we can calculate the radial force F r with [START_REF] Shah | Anisotropic permeable porous facing abd slip velocity squeeze film in axially undefined journal bearing with ferrofluid lubricant[END_REF].

F r = E m (1) -E m (2) 2∆r (28)
where ∆r = 0.001 m is the radial distance in length for which the inner cylinder is moved away. The radial stiffness k r is determined by [START_REF] Cunha | A new boundary integral formulation to describe three-dimensional motions of interfaces between magnetic fluids[END_REF].

k r = F r ∆r ( 29
)
The numerical application gives k r = 5.6 N/mm. As a conclusion, for a ferrofluid seal whose thickness equals 0.3 mm, the radial stiffness is k r = 5.6 N/mm.

Decentering supported by the ferrofluid seal

The determination of the maximum possible decentering depends on the quantity of ferrofluid used in the seal. A numerical application is given here to illustrate that. The ferrofluid quantity we use corresponds to the case when its thickness equals 0.4 mm. With the inner cylinder, the air gap equals 0.3 mm. It is noted that the inner cylinder is radially moved away. As some ferrofluid is moved off from the place where the ferrofluid seal is crushed, the ferrofluid seal thickness increases at the opposite from the place where the ferrofluid seal is crushed. The moving off centre for which a perforation can appear in the seal corresponds to 0.14 mm. As a conclusion, in this configuration with r in = 25 mm, r out = 28 mm and h = 3 mm, the moving off centre supported by the ferrofluid seal equals 0.14 mm.

CONCLUSION

This paper has presented an analytical way of studying analytically the mechanical properties of a ferrofluid seal. This static study is based on the potential energy criterion. Moreover, the magnetic field created by the ring permanent magnets being very intense, the magnetic field created by the magnetic particles is omitted. The expression of the magnetic pressure has been determined analytically and has been plotted in order to see where the potential energy is the most important. The shape of the ferrofluid seal depends on the quantity of ferrofluid which is used and we have shown that for small quantities, the ferrofluid seal contour can be described in terms of an equation of ellipse. Moreover, if the thickness of the ferrofluid seal is greater than 6 mm, the ferrofluid seal spreads over the whole inner faces of the ring permanent magnets. It can be noted that the shape of the ferrofluid seal changes when the inner cylinder is taken into account. This change comes with an energy reduction which depends on the radial displacement of the cylinder, so on the seal crushing. The ferrofluid seal capacity has been determined and we have shown that it depends on the thickness of the ferrofluid seal. At last, a way of obtaining the radial stiffness of a ferrofluid seal is presented and we have found a radial stiffness which equals 5.6 N/mm. Furthermore, the moving off centre supported by the ferrofluid seal has been studied and we have found that it depends on the quantity of ferrofluid used in the seal. These results can be useful to design many engineering and industrial applications which use ferrofluid seals.

  MECHANICAL PROPERTIES OF A FERROFLUID SEAL: THREE-DIMENSIONAL ANALYTICAL STUDY BASED ON THE COULOMBIAN MODEL R. Ravaud, G. Lemarquand Laboratoire d'Acoustique de l'Universite du Maine, UMR CNRS 6613 Avenue Olivier Messiaen, 72085 Le Mans, France guy.lemarquand@univ-lemans.fr

Figure 2 .

 2 Figure 2. Magnetic field created by the two outer ring permanent magnets; r in = 25 mm, h = 3 mm, σ * = 1 T

e 4

 4 +du ], -d+e 4 d+e 4 d √ -d -e 4 (-c + e 4 ) d(1+u) e 4 +du √ 1 -u 2

3 Figure 3 .Figure 3

 333 Figure 3. Three-dimensional representation of the magnetic pressure of the ferrofluid seal

Figure 4 .

 4 Figure 4. Representation of growing ferrofluid volumes under a magnetic field from the warmest areas to the coolest ones; the warmest colors correspond to the region in space where the magnetic energy is the strongest

Figure 5

 5 Figure 5. Representation of growing ferrofluid volumes under a magnetic field from the warmest areas to the coolest ones: the warmest colors correspond to the region in space where the magnetic energy is the strongest

Figure 6 .Table 3 .

 63 Figure 6. Tree-dimensional representation of the ferrofluid seal volume; the ring inner radius equals 25 mm. J = 1 T, h = 3 mm

Figure 7 .

 7 Figure 7. Ferrofluid seal which is crushed

Figure 8 .Figure 9 .Figure 10 .

 8910 Figure 8. Geometry: two ring permanent magnets and an inner nonmagnetic cylinder; air gap in the seal

Figure 11 .Figure 12 .

 1112 Figure 11. Radially moved away cylinder, crushed ferrofluid


G Masson

S Cogan

N Bouhaddi

J P Lombard

J Bonini

Parameterized reduced models for efficient optimization of structural dynamic behavior

Substructure synthesis is a method of model order reduction which is generally more efficient, computationally speaking, than analyzing the complete structural system. However, these methods are not necessarily well adapted for use within an optimization process since they do not preserve the fidelity of the reduced models when structural modifications within the reduced substructure are introduced. As a result, a costly model reduction must be performed at each iteration step.

This paper presents a new method to improve standard reduction methods by taking into account a priori knowledge of the potential structural modifications. Indeed, this information proves to be salutary in creating a single enriched model reduction transformation that preserves the precision of the reduced substructure model throughout the optimization process. The proposed approach consists in extending the standard transformation matrix by a set of static residual vectors which are optimized with respect to the design variables to be modified. The proposed method can be used with a variety component mode synthesis approaches with any type of substructure natural modes: free-free, cantilever or hybrid modes. The proposed methodology is illustrated on the basis of a simulated test case taken from the aero-engine industry.

Introduction

The ever increasing demand for faster engineering analysis in the design process has resulted in a substantial amount of research and development on faster and more accurate approximate reanalysis method. The difficulty of any model reduction procedure lies in how to complete the representation basis in order to reduce truncation effects. Two main classes of methodologies can be found in the literature. The first class seeks to generate a set of Ritz vectors capable of representing with precision the structural behavior under a wide variety of structural modifications. For example, Balmès 1 10 studied the possibility of using a constant basis of Ritz vectors to create parametric families of reduced models whereas Bouazzouni A., et al 2 9 developed a method for optimally constructing additional vectors by using the dynamic behavior of the structure before modification combined with the a priori knowledge of the design variables. Both of these approaches have already been used effectively in an industrial context. The second class of reduction methods are based on a high order polynomial expansions of model responses about a nominal point in parameter space. This approach proves to be particularly interesting for small number of design variables and can be used for topological optimization. In this paper, we will extend the approach developed in [START_REF] Bouazzouni | Selecting a ritz basis for reanalysis of frequency response functions of modified structures[END_REF] for use with substructure synthesis techniques. The latter represent a economic means for evaluating the structural behavior of complex mechanical assemblies and are known collectively as component mode synthesis methods. This approach represents a structure as an assembly of individually reduced substructures called superelements. For example, the Craig-Bampton (CB) method is one of many techniques of component mode synthesis used intensively in the aerospace industry 3 4 . However, these methods are not always well adapted to the industrial problem. This is especially true when parametric studies are to be performed with respect to design variables contained within individual superelements. The designer has the choice of either re-using the nominal model reduction transformation or performing a new superelement analysis for each modified component. The first option generally leads to inaccurate results while the latter is often impracticable due to cost considerations.

We propose in this paper a new method to improve the standard superelement reduction methods by taking into account an a priori knowledge of the potentially modifiable design variables with the objective of constructing a unique reduction transformation matrix which will preserve the predictiveness of the superelement for a wide range of structural modifications. The approach consists in generalizing the procedure proposed in 2 to individual structural components or substructures by completing the standard reduction transformation of the component by a set of optimized static residual vectors which depend on the set of design parameters retained for the component in question. The global parameterized model is thus constituted of an assembly of individually reduced and parameterized component models. The potential of the proposed methodology is demonstrated on the basis of a simulated numerical test case taken from the aero-engine industry.

Model reduction problem

The basic philosophy behind the creation of a superelement from the global stiffness and mass matrices of a structural component is briefly described below in the context of a linear conservative elastodynamic substructure. The homogeneous equation of motion for such a substructure is give by:

( ) 0 = - υ υ λ y M K (1)
Where M;K represent the discrete mass and stiffness matrices of the system and υ υ λ ; y are a given pair of eigensolutions. The generic substructure model reduction is performed via a transformation matrix T 0 such that:

υ υ φ q T 0 = (2)
Where q ν is a vector of generalized coordinates representing the contribution of each column of T 0 . Note that some of these coordinates may correspond to physical degrees of freedom. Substituting (2) into (1) and premultiplying by the transpose of T 0 yields:

( ) 0 ˆ= - υ υ λ q M K (3)
where K ˆ;

n n R M , ˆ∈
are the reduced order stiffness and mass matrices. For example, in the case of a Craig-Bampton reduction [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF] , the substructure degrees of freedom (dof) are partitioned into junction and interior dof sets. The junction dof correspond to those dof used to assemble the substructure to the so-called residual structure. The partitioned mass and stiffness matrices of the substructure thus have the following form:

      =       = ii ij ji jj ii ij ji jj M M M M M K K K K K (4)
the transformation matrix T 0 has the following form:

      Φ = ij G I T 0 0 ( 5
)
Where

ij ij K K G ii 1 - - =
and Φ are a subset of eigenvectors of the slave eigenvalue problem

[ ] 0 = - υ υ φ σ ii ii M K .
The corresponding superelement matrices of the reduced substructure are given by:

        Φ Φ Φ + Φ Φ + Φ + + =       Φ Φ + = ii T ij T ij ii T ji ii T ij jj ij T ij ij ii T ij ii T ij ji jj M M G M M M G M M G G M G M K G K K K ˆ0 0 ˆ(6)
The Craig-Bampton superelement matrices represent the dynamic behavior of the substructure as seen through the junction dof. The static behavior is exact for all junction loads. The question we are concerned with here is: what happens if a design variable (e.g. a plate thickness or beam property) of a superelement is modified? In order to save the cost of calculating a new T 0 , we might be consider using the same transformation (evaluated on the basis of a nominal model) on the modified global matrices of the substructure. Experience shows that the precision of this solution degenerates rapidly with the amplitude of the perturbation. We propose an alternative solution that allows the transformation matrix to completed with a new basis of vectors which are optimally chosen with respect to the design variables to be modified

Optimal static residual vectors

The strategy described in [START_REF] Bouazzouni | Selecting a ritz basis for reanalysis of frequency response functions of modified structures[END_REF] for improving the robustness of the reduction transformation with respect to structural modifications will be briefly reviewed in this section. The novelty of this approach lies in its use of an a priori knowledge of the nature and localization of the potential design modifications while retaining the amplitude of the modification as a variable. The objective is to construct an optimal set of Ritz vectors, using the nominal model, to complete the standard reduction transformation.

Candidate design modifications

We define a parametric correction of the impedance Z of the substructure model as follows :

M K Z ∆ - ∆ = ∆ 2 ω (7)
By introducing notion of zones or groups of finite elements, equation ( 7) can be written as follows :

[ ] ∑ = ∆ - ∆ = ∆ np i i zone i i zone i p M p K Z 1 2 ) ( ) ( ω (8)
where ∆p i is the variation of the np design parameter p i .

In general, the correction matrix is a nonlinear function of the design parameters p i and, for example, the stiffness modification of a zone can be expressed as :

zone i i i i i zone i K p p p p K α α α α ) ( ) ( ∑ ∆ = ∆ ∆ (9)
The homogeneous equations of motion of the modified substructure is given by :

[ ] 0 ) ( ) ( ) ( = ∆ + ω ω ω y Z Z (10)

Force basis associated to a structural modification

By introducing the notion of force

) (ω ∆ f associated with structural modifications ) (ω Z ∆
, equation (10) can be rewritten:

) ( ) ( ) ( ω ω ω ∆ = f y Z (11)
where:

) ( ) ( ) ( ω ω ω y Z f ∆ - = ∆ (12)
The vector ) (ω ∆ f represents the forces acting on the nominal structure due to the structural modifications. Given that the response vector y is unknown, the force vector cannot be determined exactly either. The essential idea of the methodology proposed in 2 is expressed in the following steps: 1. use (12) to generate a basis of force vectors that, if it does not contain the exact force vector with respect to a specific design modification, will at least represent a space containing these vectors. This is accomplished by injecting known nominal model response vectors into (12). 2. The resulting force basis is then used to generate static response vectors, once again on the basis of the nominal model. 3. The first two steps are repeated for each candidate design parameter. In practice, many different types of model response vectors may be injected into (12), including global substructure modes or sensitivity vectors, slave system modes, static vectors, etc. These forces are obviously frequency dependent and are generally evaluated at the frequency corresponding to the response vector in question. Let B∈R N,m denote the basis of response vectors used for this purpose, hence:

) ( ) ( ω ω Bc y = (13)
For a modification of the design parameter p i , equation ( 12) becomes:

[ ] m m i i i b Z b Z F )( ) ( 1 1 ω ω ∆ - ∆ - = ∆ L (14)
The final force basis, representative of the set of global set of candidate structural modifications, is then obtained by concatenation of the series of bases generated for each zone:

[ ] np F F F F ∆ ∆ ∆ ∆ = ... 2 1
(15)

Enriched transformation matrix

Once we have built our force basis, we can compute a series of response vectors which will be used to complete the standard model reduction transformation matrix. These responses are obviously frequency dependent and in some cases it is worthwhile Once again, these responses are generally frequency dependent and in some cases it is worthwhile to calculate them explicitly as in [START_REF] Grihon | Une méthode de réduction de modèle pour la réanalyse modale et aéroélastique de structure d'avion[END_REF] . However, in order to avoid re-decomposing the global stiffness matrix of the substructure, the static solution is often evaluated instead:

∆ - ∆ = FK S 1 (16)
Practically speaking, the basis S ∆ is rarely of full rank and must be reconditioned before use. The form this reconditioning takes depends on both the contents of both S ∆ and T 0 as will be discussed in the next section. Letting Q∈R N,q be the reconditioned matrix, the improved transformation matrix T is given by: [ ]

Q T T 0 = (17)
where T 0 is a classical transformation matrix. In the case where T 0 is the Craig-Bampton transformation, the enriched superelement has the form:

          Φ Φ Φ Φ Φ + Φ Φ + Φ + + =           Φ Φ Φ Φ + = Q M Q M Q Q M M M G M M M G M M G G M G M Q K Q K Q Q K K G K K K ii T ii T ii T ii T ij T ij ii T ji ii T ij jj ij T ij ij ii T ij ii T ii T ii T ii T ij ji jj 0 0 ˆ0 0 0 0 ˆ(18)
In the case of the Craig-Bampton reduction, the complementary basis Q is calculated with respect to the slave system. Hence, the displacements on the junction dof are zero.

Remarks

Several practical remarks are in order concerning the practical application of equation ( 17):

• The static responses S ∆ are often linearly dependent and it is important to perform an adequate reconditioning of this matrix before use. The form of this pre-processing depends on the contents of T 0. For example, if part or all of T 0 the modal matrix of the substructure, then it is advantageous to replace S ∆ by the corresponding static residual matrix formed by removing the static contribution of the retained modes.

• When possible, it is important not to combine the contributions of several zones to construct a single (sub)basis. This has the effect of coupling otherwise independent parameters.

Industrial application

Description of test case

The proposed methodology will be illustrated on a subassembly of an aero-engine model (figure 1) comprising four Craig-Bampton superelements and a residual structure containing, by definition, the ensemble of junction dof for all superelements as well as physical dof corresponding to unreduced structural elements. The numerical model contains a total of about 40000 dof and is under free-free boundary conditions. Substructure 1, a casing, contains 9000 dof with 300 interface dof while substructures 2 and 3 are both stators containing respectively 4200 dof with 600 interface degrees of freedom and 4400 dof with 600 interface dof. Substructure 4 contains 6000 dof with 60 interface dof and is also a casing. A analysis range is defined which contains 6 rigid body modes and 12 flexible modes in the initial configuration. Superelements for all four components have been created using the Craig-Bampton technique.

In order to illustrate the proposed methodology, we have perturbed the nominal substructure models as follows:

-concerning the first substructure, we define two zones where the Young's modulus are modified by a factor 2.

-concerning the second and the third substructures which are two tubes we choose to introduce local modifications. We define two zones for each substructure where the Young's modulus are modified respectively by a factor 2 and 0.5.

-and for the last substructure, we define a zone where the thickness is the modification parameter. The nominal values of the shell element thickness are divided by 2. These structural perturbations lead to maximum frequency shifts on the order of 20%. We will compare the precision between two types of calculation: 1. Re-using the standard Craig-Bampton transformation of each the nominal substructure models on the corresponding perturbed models. 2. Enrichment of the standard Craig-Bampton transformation by a set of optimally selected Ritz vectors. In order to demonstrate the improved convergence properties optimal basis vectors, the dimensions of the two transformation matrices are chosen to be identical, that is to say, the number of slave eigenmodes included in Φ are increased for the standard Craig-Bampton reduction. The dimension of the reduction bases for both methods are reported in table 1 and2

Comparison criteria

The relative precision of the two analyses will be compared on the basis of three criteria: 1. Relative frequency error:

100 * exact exact R f f f f - = ε
where:

• f exact is an eigenfrequency of the assembled perturbed system system calculated using an exact update of the modified superelements. • f R is an eigenfrequency of the assembled system system calculated using either the standard or enriched Craig-Bampton transformation on the perturbed model. where:

• y exact ∈R N,1 an eigenvector of the assembled perturbed system calculated using an exact update of the modified superelements. • y R ∈R N,1 an eigenvector of the assembly system calculated using either the standard or enriched Craig-Bampton transformation on the perturbed model.

Transfer function

This criteria is used in a qualitative way only and allows the impact the predictive errors to be interpreted in terms of amplitude envelopes.

Results

Table 3 reports the relative frequency error and MAC criteria for the two analyses. We note a general improvement in the precision of all eigenfrequencies and eigenvectors. The remaining modes are local to individual substructures. Figures 2 and3 show the exact (nominal and perturbed) and approximate transfer functions for two different collocated observation-excitation points. In both cases, the modal damping coefficient was defined to be 2%. Note that axis labels have been purposely omitted. The dotted and dashed lines represent the exact transfer function of the perturbed and nominal models, respectively, while the dot-dashed and continuous lines represent the standard and enriched Craig-Bampton results, respectively. We can observe that the standard Craig-Bampton transformation largely underestimates the magnitude of the frequency shifts of the perturbed structure while the proposed methodology yields a significantly more accurate prediction.

Conclusions

In this paper, we have extended an existing methodology for the approximate reanalysis of linear elastodynamic behavior to superelement technology. This approach allows design parameters of a superelement to be modified, for example in the context of an optimization algorithm, without having to perform a complete superelement analysis at each point in parameter space. The method is based on the enrichment of the standard superelement reduction transformation by a basis of representation which is optimized with respect to the design parameters to be modified. The proposed methodology can be integrated in a variety of component mode synthesis techniques and we have illustrated it use in the context of the Craig-Bampton superelement. An industrial test case derived from the aero-engine industry was used to illustrate the efficiency of the optimized basis in representing the modified behavior of a perturbed structure. The proposed approach is particularly effective when dealing with very large models having a large number of subassemblies. The improved model reduction transformation can be prepared in advance for each substructure and then used directly during the optimization procedure, thus avoiding the exact recalculation of the modified superelements. Current work is in progress to take into account local non-linearities at the interfaces of superelements.

  Figure 1

2 .

 2 Modal assurance criteria (MAC):

  dotted line represent the exact transfer function of the nominal model dashed line represent the exact transfer function of the perturbed model continuous line represent the transfer function of the enriched Craig-Bampton results dot-dashed line represent the transfer function of the standard Craig-Bampton results dotted line represent the exact transfer function of the nominal model dashed line represent the exact transfer function of the perturbed model continuous line represent the transfer function of the enriched Craig-Bampton results dot-dashed line represent the transfer function of the standard Craig-Bampton results

Figure 2 :

 2 Figure 2 : Frequency Response Function

Figure 3 :

 3 Figure 3 : Frequency Response Function

Table 1 :

 1 . Composition of the standard Craig-

			Substructure
		1	2	3	4
	Constraint modes 300	600	600	60
	Normal modes 20	40	40	30
	Optimal vectors 0	0	0	0
	Reduced model	320	640	640	90
	size
	Bampton reduction
			Substructure
		1	2	3	4
	Constraint modes 300	600	600	60
	Normal modes 15	20	20	20
	Optimal vectors 5	20	20	10
	Reduced model	320	640	640	90
	size

Table 2 :

 2 Composition of the enriched Craig-Bampton reduction

Acknowledgements

The authors would like to thank the SNECMA, Villaroche for their support in developing and testing the technology presented in this paper.

Standard Craig-Bampton


Abstract-This paper presents some improved analytical expressions of the magnetic field produced by arc-shaped permanent magnets whose polarization is radial with the amperian current model. First, we show that the radial component of the magnetic field produced by a ring permanent magnet whose polarization is radial can be expressed entirely in terms of elliptic integrals. Such an expression is useful for optimization purposes. We also present a semi-analytical expression of the axial component produced by the same configuration. For this component, we discuss the terms that are difficult to integrate analytically and compare our expression with the one established by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]. In the second part of this paper, we use the amperian current model for calculating the magnetic field produced by a tile permanent magnet radially magnetized. This method was in fact still employed by Furlani for calculating the magnetic field produced by radially polarized cylinders. We show that it is possible to obtain a fully analytical expression of the radial component based on elliptic integrals. In addition, we show that the amperian current model allows us to obtain a fully analytical expression of the azimuthal component. All the expressions determined in this paper are compared with the ones established by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF] or in previous works carried out by the authors.

INTRODUCTION

More and more analytical approaches were proposed by many authors [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]- [START_REF] Furlani | Field analysis and optimization of ndfeb axial field permanent magnet motors[END_REF] for calculating the magnetic field created by arc-shaped permanent magnets whose polarizations can be radial [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF]- [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF] or axial [START_REF] Furlani | A three-dimensional field solution for permanent-magnet axial-field motors[END_REF]- [START_REF] Perigo | General expressions for the magnetic flux density produced by axially magnetized toroidal permanent magnets[END_REF]. Strictly speaking, there are mainly four analytical methods that are used for calculating the three magnetic field components produced by tile permanent magnets or ring permanent magnets radially magnetized. These four methods derived directly from the Maxwell's equations and the properties of the magnetic field shape.

Consequently, some authors used the coulombian model or the magnetic potential [START_REF] Babic | Calculation improvement of 3d linear magnetostatic field based on fictitious magnetic surface charge[END_REF]- [START_REF] Ravaud | Analytical expression of the magnetic field created by tile permanent magnets tangentially magnetized and radial currents in massive disks[END_REF] while other authors used the amperian current model and the Biot-Savart Law [START_REF] Azzerboni | Magnetic field evaluation for disk conductors[END_REF]- [START_REF] Babic | New expression for calculating magnetic fields due to current-carrying solid conductors[END_REF] for calculating the magnetic field created by currents or permanent magnets. However, these analytical methods are not suitable for the study of iron-core structures. This is why some alternative methods, either based on series expansions [START_REF] Jian | Analytical calculation of magnetic field distribution in coaxial magnetic gears[END_REF] or finite-element methods [START_REF] Chau | Transient analysis of coaxial magnetic gears using finite element comodeling[END_REF] are used for the design of permanent-magnet topologies [START_REF] Chau | Design of a magnetic-geared outer-rotor permanent-magnet brushless motor for electric vehicles[END_REF].

The coulombian model implies the existence of magnetic charges located on the faces of a magnet and inside it. The amperian current model implies the determination of the vector potential produced by fictitious currents flowing around or inside a magnet. The two approaches are in fact equivalent because the magnetic field produced by a solenoid or an axially magnetized cylindrical magnet has the same shape but generally not the same value. It is well known that magnetic fields produced by permanent magnets are generally higher than the ones produced by currents.

The problem is thus to guess what model is the most appropriate for calculating the three components of the magnetic field produced by permanent magnets. It does not seem to be more difficult to use the amperian current model rather than the coulombian model for calculating the magnetic field created by parallelepiped magnets. Indeed, the coulombian model applied to a parallelepiped magnet whose polarization is directed along the normal direction implies two times the calculation of two surface integrals whereas the amperian current model implies four times the calculation of two surface integrals.

However, the two analytical methods lead to a fully analytical expression of the magnetic field [START_REF] Akoun | 3d analytical calculation of the forces exerted between two cuboidal magnets[END_REF]. This is why this configuration is often used by authors for optimization purposes: it has a very low computational cost [START_REF] Yonnet | Permanent magnet bearings and couplings[END_REF]- [START_REF] Charpentier | Calculation of ironless permanent magnet coupling using semi-numerical magnetic pole theory method[END_REF]. However, as shown in [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF], this method can only be used in the near-field if parallelepiped magnets are used for representing arc-shaped permanent magnets. It is emphasized here that numerous approaches based on series expansions have been proposed [START_REF] Selvaggi | Computation of the three-dimensional magnetic field from solid permanent-magnet bipolar cylinders by employing toroidal harmonics[END_REF]- [START_REF] Zhilichev | Calculation of magnetic field of tubular permanent magnet assemblies in cylindrical bipolar coordinates[END_REF] while semi-analytical approaches have led to compute the magnetic field in all points in space [START_REF] Rakotoarison | Using coulombian approach for modeling scalar potential and magnetic field of a permanent magnet with radial polarization[END_REF]. Moreover, other analytical methods allow the determination of the force between coil currents [START_REF] Conway | Noncoaxial inductance calculations without the vector potential for axisymmetric coils and planar coils[END_REF] [START_REF] Kim | Mutual inductance of noncoaxial circular coils with constant current density[END_REF].

For arc-shaped permanent magnets, it seems to be more difficult to guess what model is the most appropriate. For arc-shaped permanent magnets whose polarization is radial, some new analytical expressions are derived from the vector potential, not only for ring but also for tile permanent magnets. These expressions are compared with the ones established by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF], Babic [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF], Rakotoarison [33] and Ravaud [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. These expressions have been found with the help of Mathematica and Maple. Moreover, we discuss the validity and the computational cost of all these expressions. A Mathematica file is available online for showing the validity of the expressions obtained and the equivalence between the amperian current and couloumbian models [36].

The first part of this paper deals with the calculation of the magnetic field created by ring permanent magnets radially magnetized and the second part of this paper deals with the magnetic field produced by tile permanent magnets radially magnetized. By using the amperian current model, we show how some expressions can be significantly improved and why some expressions cannot still be expressed entirely in terms of special functions.

RING PERMANENT MAGNETS WHOSE POLARIZATION IS RADIAL

Notation and geometry

In this first part, we study the magnetic field produced by ring permanent magnets whose polarization is radial. For this purpose, let us consider the geometry shown in Fig 1 . Figure 1-A represents a ring permanent magnet radially magnetized according to the coulombian model. The ring inner radius is r 1 and its outer one is r 2 . Its height is h = z 2 -z 1 . Its magnetic polarization J is directed towards 0. Some fictitious charge surface densities are located on the inner and outer faces of the ring and a charge volume density appears inside the ring. Consequently, as mentioned in [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF], up to three integrals must be determined for obtaining an analytical expression of the three magnetic field components created by this ring permanent magnet. The contribution of the charge surface densities can be expressed in terms of elliptic integrals of the first, second and third kind [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF] but the contribution of the charge volume density can be expressed in an analytical part based on elliptic functions and a numerical part [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF]. Consequently, the coulombian model applied to this configuration does not seem to be the best model for obtaining a fully analytical expression of the magnetic field.

In this paper, we use the amperian current model for calculating the magnetic field produced by a ring permanent magnet radially The magnetic polarization J is transformed into a surface density currents k flowing on the upper and lower faces of the ring. Such a model has been employed for example by Lang [START_REF] Lang | Fast calculation method for the forces and stiffnesses of permanent-magnet bearings[END_REF] in the case of axially magnetized ring permanent magnets.

The ratio between J and k is expressed as follows:

k = J µ 0 (1)
By using the amperian current model in the case of radially magnetized ring permanent magnets, the expression of the radial component can be improved significantly because it can be expressed only in terms of elliptic integrals and no further numerical integrations are required. Consequently, the time necessary to calculate this component is lower with the amperian current model than the coulombian model. However, the two models are equivalent and give the same results as it is shown in a next section. We present now the new magnetic field expressions determined with the amperian current model. The analytical expression of the radial component can be determined by using either the vector potential or directly the Biot-Savart Law. Indeed, by using the analogy between the coulombian model and the amperian current model, the expression of the radial component can be expressed in terms of elliptic integrals. By denoting r the observation point and r a point located on the charge distribution on the magnet, it is useful to define the reciprocal distance r -r as follows: 1

r -r = 1 r 2 + r2 -2rr cos( θ) + (z -z) 2 (2)
By using the amperian current model, and as stated in [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF], the potential vector created by a ring permanent magnet whose polarization is radial can be expressed as follows:

A( r) = µ 0 4π   Sup k( r ) r -r d Sup - S down k( r ) r -r d Sdown   (3)
where S up is the upper surface of the ring (z = z 2 ) and S down is the lower surface of the ring (z = z 1 ). Moreover, it is to be noted that in this analogy, the current volume density k v is 0 because the polarization is supposed to be perfectly radial.

k v = ∇ ∧ J µ 0 = 0 (4)
Thus, the magnetic field can be obtained by the following expressions:

H( r) = 1 µ 0 ∇ ∧ A( r) ( 5
)
As there is a symmetry according θ in our configuration, the azimuthal component H θ (r, z) equals zero and the magnetic field depends only on r and z, that is :

H( r) = H(r, z) = H r (r, z) u r + H z (r, z) u z (6)
For the rest of this paper, we use the operator 2 i=1 2 k=1 (-1) i+k [•] for the case of ring permanent magnets whose polarization is radial. However, according to the convention chosen for calculating the three magnetic components, this operator can be replaced by

2 i=1 2 k=1 (-1) i+k+1 [•]
. The modulus of the magnetic field created is always the same but the sign of the magnetic components depend directly on the way we define the axes. Such a way of writting the magnetic field expressions has also been adopted by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF] or Babic [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF]. The only influence of the choice of (-1) i+k [•], (-1) i+k+1 [•] or (-1) i+j+k [•] and (-1) i+j+k+1 [•] is the change in sign of the magnetic field component expressions.

Radial component determined with the amperian current model

The radial component H r (r, z) can be determined by calculating the projection of H( r) along u r :

H r (r, z) = H( r). u r = 1 µ 0 ∇ ∧ A( r) . u r (7)
The interest of using the amperian current model for this configuration lies in the fact that only two integrals must be calculated. Consequently, we obtain an analytical expression of the radial field component based on elliptic integrals. By using Mathematica, the arguments of the elliptic integrals used are defined as follows:

φ +,- 1 = (b + 2e)x bx ± √ 2 xe 2 (x -c) φ 2 = i sinh -1   -1 b + 2e b -2e cos( θ)   φ 3 = b + 2e b -2e (8)
It is noted that we use mainly the elliptic integrals of the second and third kind for the calculation of the radial component. These special functions have also been used in [5][6]. However, no further numerical integrations are required here. In short, the radial component H r (r, z) can be expressed as follows:

H r (r, z) = J 4πµ 0 2 i=1 2 k=1 (-1) i+k (g(i, k, 2π) -g(i, k, 0)) (9)
with

g(i, k, θ) = 2(z-z k )f (r 2 +(z-z k ) 2 , r 2 +r 2 i +(z-z k ) 2 , rr i , -r 2 -2(z-z k ) 2 , θ) (10) where f (a, b, e, c, x, θ) = η 2ξ 1 (2ce 2 + ξ 2 ) F * [φ 2 , φ 3 ] +η -e 2 (c -x)(bx √ 2 + 2ξ 1 ) Π * φ + 1 , φ 2 , φ 3 +η e 2 (c -x)(bx √ 2 -2ξ 1 ) Π * φ - 1 , φ 2 , φ 3 -2ηax xe 2 -ce 2 √ 2 + bξ 1 Π * φ + 1 , φ 2 , φ 3 -2ηax -xe 2 + ce 2 √ 2 + bξ 1 Π * φ - 1 , φ 2 , φ 3 (11)
where F * [x, y] and Π * [x, y, z] are the elliptic integrals of the second and third kind that have been used in previous papers [5][6].

In addition, the parameters ξ 1 , ξ 2 , η are defined as follows:

ξ 1 = e 2 x(x -c) ξ 2 = x(b 2 -2e 2 ) η = i -e 2 sin( θ) 2 (b-2e) 2 csc( θ) 2 -1 b+2e xξ 1 (2ce 2 + ξ 2 ) (12)
Therefore, by using the amperian current model applied to a ring permanent magnet radially magnetized, it is possible to obtain a fully analytical expression of the radial field component. This result is of great importance because it can be used for optimizing the radial field in ironless loudspeakers [START_REF] Lemarquand | Ironless loudspeakers[END_REF][38].

Validation of the radial component expression

A comparison has been performed for verifying the accuracy of the radial component expression with the expressions determined in [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF] and [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. It is noted that this expression can also be compared to the expression determined by Furlani by using θ s (2) -θ s (1) = 2π and by neglecting the current densities flowing in the axial direction (see [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]). We take the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1T , z = 0.0035 m. We represent in Fig 2 the radial component H r (r, z = 0.0035 m) versus the axial displacement z. Figure 2 clearly shows that the analytical expression determined in this paper is accurate and can be used instead of the one established in [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF] because it has a lower computational cost. In fact, the choice of the model used is very important. Indeed, it is more difficult to obtain a fully analytical expression of the radial component H r (r, z) in the coulombian approach rather than in the amperian approach. On the other hand, this does not imply that the amperian current model is always more appropriate than the coulombian model, as shown in [START_REF] Ravaud | Analytical expression of the magnetic field created by tile permanent magnets tangentially magnetized and radial currents in massive disks[END_REF] in which radial currents in massive disks generate a magnetic field that can be deducted more easily from the coulombian model. In fact, more generally, for cylindrical geometries, the authors feel that the analyze of the geometry studied must receive a peculiar attention according to the choice of the model used.

Axial Component H z (r, z)

The axial component H z (r, z) can be determined by calculating the projection of H( r) along u z :

H z (r, z) = H( r). u z = 1 µ 0 ∇ ∧ A( r) . u z ( 13
)
This axial component H z (r, z) is thus given by :

H z (r, z) = J 4πµ 0 2 i=1 2 k=1 (-1) i+k (I 1 (r, z) + I 2 (r, z)) (14)
with

I 1 (r, z) = -4r i (r -r i ) 2 + (z -z k ) 2 K * -4rr i (r -r i ) 2 + (z -z k ) 2 (15)
where K * [x] is the complete elliptic integral of the first kind [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF].

I 2 (r, z) = 2π 0 ln r i -r cos( θ) + r 2 + r 2 i + (z -z k ) 2 -2rr i cos( θ) d θ (16)
It is emphasized here that this expression is given in a more compact form than the one determined in a previous paper by the authors [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. Consequently, it seems to be more judicious to use the amperian current model for calculating the axial component of the magnetic field produced by a ring permanent magnet radially magnetized rather than the coulombian model. However, there is still a part of this expression that can be determined by using a numerical integration.

Validation of the expression of the axial component H z (r, z)

A comparison has been performed for verifying the accuracy of the axial component expression with the expressions determined in [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF] and [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. It is noted that this expression can also be compared to the expression determined by Furlani by using θ s (2) -θ s (1) = 2π. We take the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1T , r = 0.024 m. We represent in Fig 3 the axial component H z (r = 0.024, z) versus the axial displacement z. Figure 3 clearly shows that the analytical expression determined in this paper is accurate and can be used instead of the one established in [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF] because it has a lower computational cost and it is given in a more compact form.

TILE PERMANENT MAGNETS WHOSE POLARIZATION IS RADIAL

Notation and geometry

In this second part, we study the magnetic field produced by tile permanent magnets whose polarization is radial. For this purpose, let us consider the geometry shown in Fig 4 . Figure 4-A represents a tile permanent magnet radially magnetized according to the coulombian model. The ring inner radius is r 1 and its outer one is r 2 . Its height is h = z 2 -z 1 . Its magnetic polarization J is directed towards 0. Its angular width is θ 2 -θ 1 . Some fictitious charge surface densities are located on the inner and outer faces of the ring and a charge volume density appears inside the ring. As stated in the first section of this paper, we consider the current surface densities flowing around the tile permanent magnet.

In addition, we define the following operator for the rest of this paper:

℘ i,j,k [•] = J 4πµ 0 2 i=1 2 j=1 2 k=1 (-1) i+j+k [•] (17) 3.2. Radial Component H r (r, θ, z)
We use the same equation as Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF] for calculating the radial component of the magnetic field produced by a tile permanent magnet whose polarization is radial.

H r (r, θ, z) = J 4πµ 0 z 2 z 1 r 2 r 1 r sin(θ -θ 1 ) r -r (θ 1 , z) 3 drdz - J 4πµ 0 z 2 z 1 r 2 r 1 r sin(θ -θ 2 ) r -r (θ 2 , z) 3 drdz + J 4πµ 0 r 2 r 1 θ 2 θ 1 (z -z 2 ) cos(θ -θ) r -r ( θ, z 2 ) 3 rd θdr - J 4πµ 0 r 2 r 1 θ 2 θ 1 (z -z 1 ) cos(θ -θ) r -r ( θ, z 1 )
3 rd θdr [START_REF] Babic | New expressions for calculating the magnetic field created by radial current in massive disks[END_REF] where 1 r -r (α, β)

= 1 r 2 + r2 -2rr cos(θ -α) + (z -β) 2 (19)
It is noted that Eq (18) can be decomposed in two parts. The contribution of the charges located on the curved planes (inner and outer faces of the ring) can be expressed in terms of elliptic integrals by using (9) in which

2 i=1 2 k=1 (-1) i+k (g(i, k, 2π) -g(i, k, 0)) is transformed into 2 i=1 2 k=1 (-1) i+k (g(i, k, θ -θ 2 ) -g(i, k, θ -θ 1 )
). Thus, the radial field component H r (r, θ, z) can be decomposed as follows:

H r (r, θ, z) = J 4πµ 0 2 i=1 2 k=1 (-1) i+k (g(i, k, θ -θ 2 ) -g(i, k, θ -θ 1 )) + J 4πµ 0 z 2 z 1 r 2 r 1 r sin(θ -θ 1 ) r -r 3 drdz - J 4πµ 0 z 2 z 1 r 2 r 1 r sin(θ -θ 2 ) r -r 3 drdz (20)
The last r and z integrations can be performed analytically and we obtain the following expression for the radial field component:

H r (r, θ, z) = ℘ i,j,k [g(i, k, θ -θ j ) + k i,j,k + l i,j,k ] (21)
k i,j,k = sin(θ -θ j ) ln [Y i,j,k ] (22)
l i,j,k = -cos(θ -θ j ) arctan [X i,j,k ] (23)
Y i,j,k = (z -z k ) + r 2 + r 2 i + (z -z k ) 2 -2rr i cos(θ -θ j ) X i,j,k = (z -z k )(r i -r cos(θ -θ j )) r sin(θ -θ j ) r 2 + r 2 i + (z -z k ) 2 -2rr i cos(θ -θ j ) (24)
As a conclusion, it is possible to obtain an analytical expression of the radial field component produced by a tile permanent magnet radially magnetized. This result is interesting for the optimization of radial fields in electrical machines.

Validation of the expression of the radial field component

The expression of the radial field component produced by a tile permanent magnet radially magnetized is compared with the expression obtained by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]. We take the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1T , r = 0.02 m. We represent in Fig. 5 H r (r = 0.024, z) versus the angular displacement θ. Let us now consider the azimuthal component H θ (r, θ, z) of the magnetic field produced by a tile permanent magnet radially magnetized. We use the same equation as Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF] for calculating this component. Therefore, by using the amperian current model, the azimuthal component is expressed as follows :

H θ (r, θ, z) = J 4πµ 0 2 j=1 (-1) j    r 2 r 1 θ 2 θ 1 (z -z) sin(θ -θ) r -r j (r, θ) 3    rdrd θ + J 4πµ 0 2 j=1 (-1) j    r 2 r 1 z 2 z 1 r -r cos(θ -θ j ) r -r (r, z) 3    drdz (25) where 1 r -r j (α, β) = 1 r 2 + α 2 -2rα cos(θ -β) + (z -z j ) 2 1 r -r (α, β) = 1 r 2 + α 2 -2rα cos(θ -θ j ) + (z -β) 2 (26)
We find:

H θ (r, θ, z) = ℘ i,j,k [-cos(θ -θ j ) ln [Y i,j,k ] -sin(θ -θ j ) arctan [X i,j,k ]] +℘ i,j,k (z -z k ) r (r i -r cos(θ -θ j ) + ξ i,j,k ) (27)
where

ξ i,j,k = r 2 + r 2 i + (z -z k ) 2 -2rr i cos(θ -θ j ) (28)
Therefore, it is possible to obtain a fully analytical expression of the azimuthal field component by using the amperian current model.

Validation of the expression of the azimuthal field component

The expression of the azimuthal field component produced by a tile permanent magnet radially magnetized is compared with the expression obtained by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]. We take the following dimensions:

r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1T , r = 0.024 m, θ 2 -θ 1 = π 4 .
We represent in Fig. 7 the axial component H z (r = 0.024, θ, z = 0.002 m) versus the angular displacement θ.

Figure 6 shows that the fully analytical expression obtained in this work for azimuthal component gives the same azimuthal field as the one established by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]. However, the expression we give in this paper does not use any special functions or numerical integrals. Consequently, such an expression has a lower computational cost and is useful for optimization purposes.

Axial Component H z (r, θ, z)

The axial component H z (r, θ, z) is also determined with the same integral formulation of Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF].

H z (r, θ, z) = J 4πµ 0 2 j=1 (-1) j r 2 r 1 θ 2 θ 1 rr cos(θ -θ) r -r (α, γ) 3 drd θ - J 4πµ 0 2 j=1 (-1) j r 2 r 1 θ 2 θ 1 r2 r -r (α, γ)
r 2 + γ 2 -2rγ cos(θ -α) + (z -z j ) 2 (30)
After the two integrations, we obtain a compact form of the axial component:

H z (r, θ, z) = ℘ i,j,k -2r i (r -r i ) 2 + (z -z k ) 2 F * θ -θ j 2 , - 4rr i (r -r i ) 2 + (z -z k ) 2 + J 4πµ 0 θ 2 θ 1 ln r i -r cos(θ -θ) + ξ θ d θ (31)
with

ξ θ = r 2 + r 2 i + (z -z k ) 2 -2rr i cos(θ -θ) (32)
where F * [φ, m] gives the elliptic integral of the first kind [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF]. This expression has also been improved compared to the one determined in [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF] because it is presented in a more compact form.

Validation of the expression of the axial field component

The expression of the axial field component produced by a tile permanent magnet radially magnetized is compared with the R. Ravaud (thick line = this work ), ( Points = Furlani's work [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]) expression obtained by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]. We take the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1T , r = 0.015 m, θ 2 -θ 1 = π 4 . We represent in Fig. 7 the axial component H z (r = 0.015 m, z = 0.002 m) versus the angular displacement θ.

Figure 7 shows that the reduced expression obtained in this work for the axial component gives the same axial field as the one established by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF].

CONCLUSION

This paper has presented some improved analytical expressions of the magnetic field produced by ring and tile permanent magnets radially magnetized. The radial field is given in terms of elliptic integrals and no further numerical integrations are required. Consequently, the time necessary for calculating this radial field is very low and allows us to carry out easily parametric studies. In the case of tile permanent magnets radially magnetized, the amperian current model leads to a fully analytical expression of the azimuthal component, as in the case of the coulombian model. However, only two integrals must be determined with the amperian current model and we think that this model is more appropriate than the coulombian model for the study of this azimuthal field. The axial field produced by a tile permanent magnet or a ring permanent magnet can be expressed is in an analytical part based on elliptic integrals and a numerical part. This expression has been improved compared to the one established by Furlani but we have not succeeded in integrating analytically the last term. Nevertheless, we have presented the expression of the axial field in a compact form and this expression seems to be easier to use than the one determined by the authors in a previous paper with the coulombian model [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. In short, the amperian current model seems to be more appropriate than the coulombian model for calculating the magnetic field produced by radially magnetized arc-shaped permanent magnets.

  COMPARISON OF THE COULOMBIAN AND AMPERIAN CURRENT MODELS FOR CALCULATING THE MAGNETIC FIELD PRODUCED BY RADIALLY MAGNETIZED ARC-SHAPED PERMANENT MAGNETS R. Ravaud, G. Lemarquand Laboratoire d'Acoustique de l'Universite du Maine, UMR CNRS 6613 Avenue Olivier Messiaen, 72085 Le Mans, France

Figure 1 .

 1 Figure 1. Representation of a ring permanent magnet radially magnetized : A: coulombian model, B: amperian current model

Figure 2 .

 2 Figure 2. Representation of the radial component H r (r, z = 0.0035 m) versus the radial displacement r with the following values: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1T ; (thick line = amperian current model), (points = coulombian model)

Figure 4 -Figure 3 .Figure 4 .

 434 Figure 3. Representation of the axial component H z (r = 0.0249, z) versus the axial displacement z with the following values: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1T ; (thick line = amperian current model), (points = coulombian model)

Figure 5 .

 5 Figure 5. Representation of the radial component H r (r = 0.02, θ, z = 0.0035 m)) versus the angular displacement θ with the following dimensions: 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1T , r = 0.02 m; (thick line = this work ), ( Points = Furlani's work[1])

Figure 6 . 1 r

 61 Figure 6. Representation of the azimuthal component H θ (r = 0.024, θ, z = 0.002 m)) versus the angular displacement θ with the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z -z 1 = 0.003 m, J = 1T , r = 0.024 m, θ 2 -θ 1 = π 4 rad; (thick line = this work ), ( Points = Furlani's work[1])

Figure 7 .

 7 Figure 7. Representation of the axial component H z (r = 0.015, θ, z = 0.002 m)) versus the angular displacement θ with the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1T , r = 0.015 m, θ 2 -θ 1 = π 4 ; (thick line = this work ), ( Points = Furlani's work [1])


François Musy

Artem Napov

Yvan Notay

Ronan Perrussel
email: ronan.perrussel@ec-lyon.fr

Riccardo Scorretti

Franc ¸ois Musy

Krylov-based algebraic multigrid for edge elements

come

I. INTRODUCTION

In finite element method libraries, the linear system solvers play a key role in terms of performances in the computing time and in the memory consumption. Multigrid methods are among the most efficient iterative linear system solvers for elliptic problems [START_REF] Hackbusch | Multi-Grid Methods and Applications[END_REF]. In this work, we focus on the discretization with the lowest order edge element of a curl-curl equation [START_REF] Hackbusch | Multi-Grid Methods and Applications[END_REF] which gives rise to a linear system Ax = f . Several Algebraic MultiGrid (AMG) algorithms have been proposed for "efficiently" solving this system [START_REF] Beck | Algebraic multigrid by components splitting for edge elements on simplicial triangulations[END_REF]- [START_REF] Musy | Compatible coarse nodal and edge elements through energy functionals[END_REF]. We use the main ideas of the algorithm proposed by Reitzinger and Schöberl (RS), who were the first to propose an edge prolongation matrix satisfying a commutativity property [START_REF] Reitzinger | An algebraic multigrid method for finite element discretizations with edge elements[END_REF]. Our motivation comes from the fact that this algorithm has the fastest setup time and gives the sparsest prolongation and coarse grid matrices among the available methods. However, it also gives the poorest rate of convergence, leading, in the literature, to a non-optimal multigrid solver. We combine the RS algorithm ideas with a Krylov-based multigrid cycle in order to recover classical multigrid performance. Numerical experiments are performed on 2D problems.

curl δ curl U + γU = f on Ω ⊂ R d (d = 2 or 3),

II. COMPONENTS OF THE ALGORITHM

A. Recursive Krylov-based multigrid cycle

The multigrid preconditioning algorithm on grid k (denoted by MGp) is given by Algorithm 1 where matrix A k represents the discrete problem on grid k (grid 0 is the coarsest) and P k is the prolongation matrix from grid k -1 to grid k.

Algorithm 1: INPUT r k , OUTPUT z k = MGp(r k ) 1) Relax using smoother M k : z k ← M -1 k r k . 2) Restrict residual: r k-1 ← P T k (r k -A k z k ).
3) Compute an approximate solution θk-1 to:

A k-1 θ k-1 = r k-1 .
(2) 4) Prolongate coarse-grid correction:

z k ← z k + P k θk-1 . 5) Relax using smoother M k : z k ← z k +M -1 k (r k -A k z k ).
If k-1 is zero, the solution of ( 2) is exact. If not, a common strategy is to perform once MGp(r k-1 ) (V-cycle) or twice (Wcycle) to approximately solve [START_REF] Beck | Algebraic multigrid by components splitting for edge elements on simplicial triangulations[END_REF]. Following [START_REF] Notay | Recursive krylov-based multigrid cycles[END_REF], we propose to consider a K 2 -cycle which is a W-cycle completed by the two first iterations of a flexible Krylov subspace method. An example with the flexible conjugate gradient is given by Algorithm 2.

Algorithm 2:

INPUT r k-1 , OUTPUT θk-1 = CS(r k-1 ) 1) First iteration: d k-1 ← MGp(r k-1 ); α k-1 ← r T k-1 d k-1 d T k-1 A k-1 d k-1 ; θk-1 ← α k-1 d k-1 ; r k-1 ← r k-1 -α k-1 A k-1 d k-1 .
2) Second iteration:

c k-1 ← MGp(r k-1 ); d k-1 ← d k-1 - c T k-1 A k-1 d k-1 d T k-1 A k-1 d k-1 c k-1 ; θk-1 ← θk-1 + r T k-1 d k-1 d T k-1 A k-1 d k-1 d k-1 .
The iteration cost in time and memory of both K 2 -and W-cycle is roughly equivalent. Moreover, their convergence rate are theoretically similar [START_REF] Notay | Recursive krylov-based multigrid cycles[END_REF]. Nevertheless, in practice the K 2 -cycle has a better convergence rate than the W-cycle [START_REF] Notay | Recursive krylov-based multigrid cycles[END_REF].

B. Prolongation matrix and smoother

Reitzinger and Schöberl proposed to construct an edge prolongation matrix P edg satisfying a commutativity property:

P edg G H = G h P nod . (3)
In this equality, P nod is a nodal prolongation matrix obtained from a nodal auxiliary matrix and G h and G H are respectively fine and coarse edge-node incidence matrix. At the finest level, G h is given by the relation between vertices and edges on the finite element mesh. At the same level, the nodal auxiliary matrix B, following the proposition in [START_REF] Kaltenbacher | Algebraic multigrid methods for nodal and edge based discretizations of maxwell's equations[END_REF], contains information about the edge-node incidence, the lengths of the edges and coefficient δ from (1). The construction of P nod is then performed by the double pairwiseaggregation algorithm proposed in [START_REF] Notay | An aggregation-based algebraic multigrid method[END_REF]. Once P nod is known, the definition of the coarse edge incidence matrix G H and the edge prolongator P edg are straightforward. To apply the construction recursively, coarse grid matrices are obtained by Galerkin product:

A H = (P edg ) T AP edg , B H = (P nod ) T BP nod . (4)
The smoother is a symmetric version of the smoother proposed in [START_REF] Hiptmair | Multigrid method for Maxwell's equations[END_REF] in a geometric multigrid context. It uses on each grid the corresponding edge-node incidence matrix.

III. NUMERICAL RESULTS

The behavior of the method with an increasing size of the problem and several kinds of parameters δ and γ is studied.

A. Structured mesh

The examples are taken from [START_REF] Jones | A multigrid method for variable coefficient maxwell's equations[END_REF]. The domain is a unit square and Dirichlet boundary conditions are enforced. The mesh with triangles is structured but this fact is not used by the solver. For the parameters, three situations are considered:

1) Homogeneous parameters: δ = γ = 1.

2) Oscillating with discontinuities for δ = f (x, y) and γ = 1. Function f have the following definition:

f = C(2 + sin(40πx)) 2 (2 + cos(40πy)) 2 with C =          10 in ]0, 0.5[×]0, 0.5[, 10 4 in ]0.5, 1[×]0, 0.5[, 10 -1 in ]0, 0.5[×]0.5, 1[, 10 2 in ]0.5, 1[×]0.5, 1[.
3) Oscillating with discontinuities for δ = f (x, y) and for γ = f (y, x).

The behavior of the iterative method is evaluated by computing the average convergence rate σ est in energy norm:

σ est = er t k Aer k er t 0 Aer 0 1/(2k f )
with er k the error at the k-th iteration and k f the iteration where the stopping criterion is reached.

An examination of the convergence rate of the two-grid algorithm provides information to predict the behavior of the multigrid cycle. In Table I, the two-grid convergence rate is quasi-independent of the size of the problem and is bounded away from 1. Case 3 is the worst situation probably because the aggregation does not take into account variation of γ; it may explain the bad convergence rate for one particular mesh in this case (in bold). The number of unknowns is divided roughly by 4 between fine and coarse grids which is the best trade-off between coarsening and the overall arithmetic complexity.

With such convergence rates for the two-grid solver, it is illustrated in Table II that the convergence rates of the W-cycle (and of course of the V-cycle) deteriorates as the number of grids increases. On the contrary, the convergence rate of the K 2 -cycle has a remarkable stability. The change of cycling proposed for the RS algorithm seems to make this algorithm robust and optimal in 2D. The proposed algorithm can be straightforwardly extended to 3D problems but our actual Matlab implementation does not enable us to propose results on challenging problems and also to discuss computational time and memory requirements. We are working on an implementation in Fortran90 in order to propose valuable comparisons on realistic problems.

TABLE II COMPARISON

 II BETWEEN A W-CYCLE AND A K 2 -CYCLE SOLVER. CASE 3 ON THE MESH WITH 48896 D.O.F.. Coefficients δ and γ are those defined in Case 3. The domain is still a unit square but an unstructured mesh, refined at the center of the square, is considered. Table III confirms the results obtained for the structured meshes.

	# grids	3	4	5	6	7
	d.o.f. coarsest grid 3279	844	216	55	11
	K 2 -cycle, σest	0.68	0.68	0.68	0.68	0.68
	W-cycle, σest	0.83	0.89	0.93	0.95	0.97
	B. Unstructured mesh