Thierry Lubin
email: thierry.lubin@univ-lorraine.fr

Amir Abbas Vahaj

Akbar Rahideh

Cross section of the axial-flux reluctance magnetic coupling Design Optimization of an Axial-Flux Reluctance Magnetic Coupling based on a Two-Dimensional Semi-Analytical Model

In this paper, design optimization of an axial flux reluctance magnetic coupling is presented. The optimal design procedure is based on a 2-D semi-analytical model defined at the mean radius combined with a multi-objective genetic algorithm (NSGA-II). In order to take into account the end-effects in the radial direction, a correction factor is defined to improve the torque and the axial force determination. The obtained results are compared with those of 3-D nonlinear finite element simulations and experimental results. It is shown that the proposed semi-analytical model is very accurate and requires very little computing time.

Introduction

Magnetic couplings can be used to transmit rotational movement between two shafts through a sealed wall in order to avoid any contact with the product to be treated for safety or hygienic reasons. The torque transmission between the two shafts is carried out without contact by the interaction between magnetic fields in the air-gap.

There are two main types of magnetic couplings: synchronous and asynchronous (eddy-current), with radial and axial flux topologies. The synchronous magnetic coupling [START_REF] Hornreich | Optimal design of synchronous torque couplers[END_REF][START_REF] Charpentier | Optimal design of cylindrical air-gap synchronous permanent magnet couplings[END_REF] consists of rare-earth permanent magnets (PMs) placed on both rotors that move at the same speed [START_REF] Fontchastagner | Some coaxial magnetic couplings designed using an analytical model and an exact global optimization code[END_REF]. This kind of coupling has maximum transmittable torque, which is the highest value of torque that can be transmitted before the stall and called pull-out torque [START_REF] Lubin | Simple analytical expressions for the force and torque of axial magnetic couplings[END_REF][START_REF] Li | Ideal radial permanent magnet coupling torque density analysis[END_REF]. For eddycurrent magnetic couplings [START_REF] Mohammadi | Double-Sided permanent-magnet radial-flux eddy-current couplers: three-dimensional analytical modeling, static and transient study, and sensitivity analysis[END_REF][START_REF] Wang | A simple method for performance prediction of permanent magnet eddy current couplings using a new magnetic equivalent circuit model[END_REF], one of the two rotors consists of a conductive plate, which usually is copper. The torque transmission is due to the induced currents which depend on the relative speed between the two rotors (slip). This results in Joule losses [START_REF] Li | Analytical Analysis of an Adjustable-Speed Permanent Magnet Eddy-Current Coupling With a Non-Rotary Mechanical Flux Adjuster[END_REF][START_REF] Aberoomand | Design optimization of double-sided permanentmagnet axial eddy-current couplers for use in dynamic application[END_REF], and adequate cooling may be necessary to avoid overheating [START_REF] Lubin | 3-D analytical model for axialflux eddy current couplings and brakes under steadystate conditions[END_REF].

Among classical magnetic couplings technologies, reluctance torque magnetic couplings have received very little attention in the literature [START_REF] Nagrial | Analysis and performance of variable reluctance (V. R.) torque coupler[END_REF]. As shown in Fig. 1, an axial-flux reluctance magnetic coupling consists of two ferromagnetic discs facing each other and separated by a small air-gap. One side is equipped with sector-shaped rareearth PMs magnetized in the axial direction with alternate polarity. The other side is made with mild steel salient poles with the same number as PMs. If one of the rotors is rotated, a change of stored energy will happen in the magnetic circuit. This change is resisted by a restoring torque toward the aligned position. The angular shift between the two discs is directly related to the transmitted torque value. In steadystate, the two rotors rotate at the same angular speed.

The main advantage of the reluctance magnetic couplings compared to more classical synchronous magnetic couplings with PMs on both sides is their durability in hightemperature environments, i.e. the side with salient-poles can be used in high-temperature environments and the other side with magnets can be used in room temperature. Besides, since fewer magnets have been used, the manufacturing price will reduce. On the other hand, the main drawback is lower torque density [START_REF] Nagrial | Analysis and performance of variable reluctance (V. R.) torque coupler[END_REF].

This paper is part of recently developed studies on the analytical modeling of magnetic couplings for their design [START_REF] Ravaud | Permanent magnet couplings: field and torque three-dimensional expressions based on the Coulombian model[END_REF][START_REF] Dolisy | A new analytical torque formula for axial field permanent magnets coupling[END_REF][START_REF] Wang | A general analytical model of permanent magnet eddy current couplings[END_REF][START_REF] Dai | Analytical Modeling of Axial-Flux Permanent Magnet Eddy Current Couplings with a Slotted Conductor Topology[END_REF]. Analytical models that can be found in the literature focus on more conventional topologies, such as those with permanent magnets on both sides [START_REF] Ravaud | Permanent magnet couplings: field and torque three-dimensional expressions based on the Coulombian model[END_REF][START_REF] Dolisy | A new analytical torque formula for axial field permanent magnets coupling[END_REF] or those with a copper disc on one side based on eddy currents Fig. 2. 2-D model at the mean radius for one pole-pitch, y-z plane. [START_REF] Wang | A general analytical model of permanent magnet eddy current couplings[END_REF][START_REF] Dai | Analytical Modeling of Axial-Flux Permanent Magnet Eddy Current Couplings with a Slotted Conductor Topology[END_REF]. For these structures where all the regions (magnet, air-gap, or copper) have the same geometry and therefore the same periodicity for the analytical solution, the analytical models are relatively simple. Here, the analytical model is more complicated due to the slotting effect (Fig. 2) which the operation of this magnetic coupling is based on it. Models already have been developed for the study of reluctance magnetic couplings but they are not suitable for use in an optimization procedure because they are too simple to predict the torque correctly [START_REF] Nagrial | Analysis and performance of variable reluctance (V. R.) torque coupler[END_REF] or they require too much computation time, i.e. 3-D finite element simulations [START_REF] Wu | Analysis and design optimisation of magnetic couplings using 3D finite element modelling[END_REF].

In this paper, we propose a 2-D semi-analytical model based on the subdomains method (SDM). SDM is a semi-analytical method that consists of dividing the problem into the physical regions named subdomains. Then Maxwell's equations are solved analytically in each subdomain [START_REF] Dubas | Analytical solution of the magnetic field in permanent-magnet motors taking into account slotting effect: no-load vector potential and flux density calculation[END_REF][START_REF] Zhu | An accurate subdomain model for magnetic field computation in slotted surface-mounted permanent magnet machines[END_REF][START_REF] Rahideh | Analytical model of slotless brushless PM linear motors considering different magnetization patterns[END_REF][START_REF] Lubin | Exact analytical method for magnetic field computation in the air-gap of cylindrical electrical machines considering slotting effects[END_REF]. This method gives accurate results when the slotting effect must be taken into account. This method has been applied for the study of permanent magnet machines [START_REF] Dubas | Analytical solution of the magnetic field in permanent-magnet motors taking into account slotting effect: no-load vector potential and flux density calculation[END_REF][START_REF] Zhu | An accurate subdomain model for magnetic field computation in slotted surface-mounted permanent magnet machines[END_REF], PM linear motors considering different magnetization patterns [START_REF] Rahideh | Analytical model of slotless brushless PM linear motors considering different magnetization patterns[END_REF] and study of the reluctance motors [START_REF] Lubin | Exact analytical method for magnetic field computation in the air-gap of cylindrical electrical machines considering slotting effects[END_REF].

The semi-analytical model proposed in this paper is developed at the mean radius in Cartesian coordinate. A new correction factor is defined to take into account the 3-D effects in the radial direction. Then the model is used for multi-objective design optimization by helping the genetic algorithm. The results are compared with 3-D non-linear finite element simulations and experimental tests.

2-D Semi-Analytical Model

To simplify the analytical modeling, we will consider a 2-D model defined at the mean radius R m . As a result, the curvature effects and 3-D effects in the radial direction will be neglected. It has been shown that curvature effects can be neglected without great errors for axial-field actuators [START_REF] Dolisy | A new analytical torque formula for axial field permanent magnets coupling[END_REF], [START_REF] De La Barrière | 3-D formal resolution of Maxwell equations for the computation of the no-load flux in an axial flux permanent-magnet synchronous machine[END_REF]. On the other hand, the decay of the magnetic field on the radial edges of the magnetic coupling must be taken into account otherwise the pull-out torque and the axial force will overestimate. Hence, a correction factor will be introduced in the formula of the torque and axial force in order to take this effect into account.

2-D Semi-analytical Model in the y-z

The 2-D model of the reluctance magnetic coupling at the mean radius is shown in Fig. 2. By considering this simplification, the cylindrical structure of Fig. 1 is made equivalent to a linear reluctance magnetic coupling. The model is then defined in a Cartesian coordinate system where the y-coordinate represents the circumferential direction and the z-coordinate the axial direction. An infinite dimension is considered in the x-direction. In addition, to obtain the simplest model, the following assumptions have been made:

-The iron parts (salient-pole and back-iron) have infinite magnetic permeability. -The permanent magnets have a relative recoil permeability μ r = 1.

The first hypothesis means that the magnetic saturation on the performance of the magnetic coupling has been neglected. This is the most important assumption and the principal limitation of the proposed model. However, in the next section, it will be shown that neglecting the magnetic saturation is not an issue, because the optimal shape of the salient poles corresponds to an unsaturated condition for the studied example. The second hypothesis implies using rare-earth permanent magnets (NdFeB or SmCo), but this is necessary to obtain good performances in terms of torque density [START_REF] Nagrial | Analysis and performance of variable reluctance (V. R.) torque coupler[END_REF].

As shown in Fig. 2, the entire domain is divided into 3 regions: Region I corresponds to the PMs, Region II is for the air-gap, and Region III shows the slot between two salient-pole. To compute the torque for different relative positions between the two rotors, the magnets (region I) can be shifted by an angle δ (torque angle) from the center of a salient-pole. As the problem presents a periodicity in the ydirection, the whole domain will be limited by 0 ≤ y ≤ τ, where τ is the pole-pitch. Also, because the whole domain is current-free, a magnetic scalar potential formulation is used to solve the problem. From the Maxwell equations, we obtain the following partial differential equation for each region (i = I, II, III).

22 22 0 ii yz       for ,, i I II III  (1)
where Φ i is the magnetic scalar potential in region i. It is note worthy to say that the magnetization vector of the magnets presents only one independent component in the zdirection, hence the Laplace equation ( 1) is valid in the PMs region. According to the periodicity of the magnetic field distribution along the y-directions in the air-gap and PMs regions (i = I, II), solution of (1) is as follow:

       
where A i , B i , C i , and D i are the integration constants which can be determined from the interface conditions between the regions, and  k =k/, k is a positive odd integer The solution in the slot region is a bit more difficult to determine [START_REF] Lubin | Exact analytical method for magnetic field computation in the air-gap of cylindrical electrical machines considering slotting effects[END_REF]. For this region, we have to consider the following boundary conditions:

    1 , 0 and 1 , 0 22 III III zz                     (3)
Considering ( 1) and ( 3), the magnetic scalar potential in region III is as follows:

    1,2,3 ( , ) sin 1 2 nn N zz III III III n n y z A e B e y                   (4)
where A III and B III are the integration constants in region III and  n =n/β, n is a positive integer.

As indicated previously, the 10 integration constants are obtained from the interface conditions between the regions, which lead to a linear system of 10 equations given in the appendix. The boundary conditions at z = 0 (region I) and z = c (region III), and the interface condition between regions I and II at z = b are easy to handle.

On the other hand, the interface condition between region II and region III at z = b given by ( 5) is more difficult and corresponds to the main difficulty for this type of semianalytical model. The method has been developed with great detail in [START_REF] Dubas | Analytical solution of the magnetic field in permanent-magnet motors taking into account slotting effect: no-load vector potential and flux density calculation[END_REF]- [START_REF] Lubin | Exact analytical method for magnetic field computation in the air-gap of cylindrical electrical machines considering slotting effects[END_REF] for cylindrical or rectangular problems and will therefore not be repeated here.

        , 1 , 1 , 22 0
elsewhere

III II y b y yb                 (5)
The next step is to determine the torque and axialforce expression from the Maxwell stress tensor. A line placed in the middle of the air-gap (region II) is taken as the integration path so the electromagnetic torque T 2D and the axial force F 2D can be expressed as follows:

2 0 0 2 ( ) ( ) m D yII zII RH T p B y B y dy     (6)   22 2 0 0 ( ) ( ) zII yII D H F p B y B y dy     (7)
where B yII and B zII are the y-and z-component of the flux density in the air-gap which are derived from the magnetic scalar potential (2) as:

0 II yII B y      0 II zII B z      ( 8
)
where  0 is the vacuum permeability. Substituting ( 8) into ( 6) and ( 7), the analytical expressions for the electromagnetic torque and for the axial force become:

  22 20 1,3,5 2 K II II II II D k T p H k A D B C     (9)   22 20 1,3,5 2 K II II II II D m k H F p k A B C D R     (10)
At this point, it is worth noting that it is not possible to give a direct relation between the torque and the geometrical parameters, as it was the case for more simple magnetic couplings geometries without slotting effect [START_REF] Lubin | Simple analytical expressions for the force and torque of axial magnetic couplings[END_REF][START_REF] Lubin | 3-D analytical model for axialflux eddy current couplings and brakes under steadystate conditions[END_REF]. When the slotting effect is taken into account by using the sub-domain method [START_REF] Dubas | Analytical solution of the magnetic field in permanent-magnet motors taking into account slotting effect: no-load vector potential and flux density calculation[END_REF]- [START_REF] Lubin | Exact analytical method for magnetic field computation in the air-gap of cylindrical electrical machines considering slotting effects[END_REF], it is necessary to numerically invert a matrix to solve the system of 10 linear equations between the integration constants given in the appendix:

     M X S  (11)
   

t I I I I II II II II III III X A B C D A B C D A B  (12)
where [M] is the topological matrix, (X) is a vector of the integration constants which should be determined, and (S) is the source term which depends on the magnet distribution.

A numerical solution of ( 11) can be obtained by using mathematical software (Matlab). Therefore, the proposed model is known as semi-analytical.

Because infinite dimension has been considered in the x-direction (radial direction), the torque and the axial force obtained from ( 9), [START_REF] Lubin | 3-D analytical model for axialflux eddy current couplings and brakes under steadystate conditions[END_REF] will have overestimation compared to the real value. This is due to the decay of the flux density on the radial edge of the magnets which is not taken into account in the torque and axial force expressions. To account for this 3-D effect, the following sub-section is dedicated to the definition of a correction factor.

Radial end effects, correction factor

In [START_REF] Lubin | Simple analytical expressions for the force and torque of axial magnetic couplings[END_REF] that focused on axial flux magnetic coupling with PMs arranged on both discs, the radial dependence of the axial flux density in the air-gap has been measured. The results are recalled in Fig. 3. As it can be observed, the axial flux density drops down rapidly near the edges of the magnet whereas it was considered constant and equal to its value at the mean radius in the previous section, as shown by the thick line curve in Fig. 3. This assumption leads to an In order to take into account this 3-D effect in a simple way and have a better estimation for the torque and the axial force, a correction factor K e is introduced by carrying-out a supplementary 2-D analytical model in the xz plane as shown in Fig. 4. The objective is to obtain the axial flux density distribution along the x-direction. We have to solve a problem with two regions having the same dimension, the correction factor is estimated under a salientpole at y = 0. The domain is limited in the x-direction by imposing a perfect magnetic boundary condition at x = ± L/2. In order to have no impact of this artificial boundary on the magnetic field determination, we impose L>>H. . The axial component of the flux density defined at z = a +b is given by:

  1,3,5 sinh ( ) cos sinh 4 with cos 2 K zII k k r k k H k L x A x k L ab L B k AH kL                       (13)
Fig. 5 shows the flux density distribution along the xdirection obtained from ( 13) with H = 30 mm, a = 10 mm and b = 5 mm. It can be observed that the waveform is close to the one obtained by measurement given in Fig. 3, but with a symmetrical aspect around x = 0.

Because the magnetic force varies with the square of the flux density, the correction factor K e will be defined as: ()

H avg zII H B B x dx H    (14)
where B max is defined at x = 0 as shown in Fig. 5.

Then the torque and axial force with the 3-D effects will be computed by combining ( 9) and ( 10) with ( 14) as:

32 D e D T K T  and 32 D e D F K F  ( 15
)
Variation of the correction factor as a function of airgap length is given in Fig. 6 for H = 30 mm and H = 100 mm and considering a fixed value for the PMs thickness (a = 10 mm). These results show that the length of the airgap and the radial length of the magnets have a significant impact on the value of the correction factor.

Comparison with 3-D Finite Element Simulation

In order to show the effectiveness of the proposed semi-analytical model, the results are compared with those obtained from 3-D finite-element simulations by considering the actual geometry of the reluctance magnetic coupling (Fig. 1). For the 3-D finite element simulations, we have used COMSOL Multiphysics® software with a scalar potential formulation (magnetostatic). The non-linear magnetic property (B-H curve) for the salient-pole has been considered in the FE analysis. The B-H characteristic is shown in Fig. 7. The mesh in the different regions has been refined until convergent results are obtained. The meshes in the magnet and in the iron part are shown in Fig. 8. Only 1 pole of the magnetic coupling has been considered in the FE analysis with anti-periodic boundary conditions in the circumferential direction. The air surrounding the magnetic coupling is considered in the FE simulation (an infinite box surround the system).

The geometrical parameters are given in Table 1. These values correspond to the ones obtained during an optimization procedure that will be presented in the next part.

The contours of magnetic flux density in the iron part under no load condition ( = 0) obtained with the 3-D FE simulation is shown in Fig. 9. The flux density inside the teeth is around 1.4T, which is slightly below the knee point of the B-H curve shown in Fig. 7.

Flux Density Distribution in the Air-Gap

Fig. 10 shows the flux density distribution along the y-direction in the middle of the air-gap, at the mean radius, under no-load (δ = 0°) and full-load condition (δ = 9°). Due to the presence of the salient-pole, the flux density in the air gap has large variations, especially for the y-component. We can observe that very good agreements are obtained between the analytical results and the 3-D FE simulations

Pull-Out Torque

For a reluctance magnetic coupling [START_REF] Nagrial | Analysis and performance of variable reluctance (V. R.) torque coupler[END_REF], the torque is a function of sin(2pδ), where δ is the torque angle, whereas it is a function of sin(pδ) for a magnetic coupling with PMs on the two discs [START_REF] Lubin | Simple analytical expressions for the force and torque of axial magnetic couplings[END_REF]. This means that the pull-out torque T max is reached for δ = /4p that corresponds to δ = 9° for the studied coupling with p = 5.

Fig. 14. Torque versus the slot-opening to pole pitch ratio.

The air-gap length has a large influence on the performance of a magnetic coupling. Its value is imposed by the thickness of the sealing barrier between the two mechanical shafts. In order to show the effectiveness of the correction factor K e on the torque prediction, Fig. 11 shows the pull-out torque as a function of air-gap length obtained without considering the correction factor [START_REF] Aberoomand | Design optimization of double-sided permanentmagnet axial eddy-current couplers for use in dynamic application[END_REF] and taking it into account [START_REF] Dai | Analytical Modeling of Axial-Flux Permanent Magnet Eddy Current Couplings with a Slotted Conductor Topology[END_REF].

One can observe that the pull-out torque is well predicted by using the proposed semi-analytical model, even with changing the air-gap length. This result clearly shows the obtained gain in precision thanks to the correction factor. The error exceeds 30% without using the correction factor (9) and it is reduced to less than 10% when the correction factor is used [START_REF] Dai | Analytical Modeling of Axial-Flux Permanent Magnet Eddy Current Couplings with a Slotted Conductor Topology[END_REF]. Another important result is a significant decrease in magnetic coupling performance as the air-gap increases. The pull-out torque goes from a value of 20 Nm for an air-gap of 2 mm to about 7 Nm for an air-gap of 6 mm, which corresponds to a decrease of around 65%.

Axial Force

The main disadvantage of an axial-flux magnetic coupling concerns the axial force between the two discs which affect the bearing losses of the drive [START_REF] Lubin | Simple analytical expressions for the force and torque of axial magnetic couplings[END_REF]. It is therefore important to predict it accurately. The maximum value for the axial force corresponds to the no-load condition (i.e. δ = 0°). Fig. 12 shows the axial force as a function of the air-gap length with or without the correction factor. For the studied coupling, the axial force reaches a value of about 1000 N for an air-gap of 3 mm. Again, the correction factor greatly improves the predictions compared to 3-D nonlinear FE simulations.

Influence of Geometrical Parameters

In this section, the impact of some geometrical parameters on the performance of the magnetic coupling is studied. The main objective is to verify the correctness of the semi-analytical model in predicting the trend, compare to the 3-D non-linear FE results. The geometrical parameters are those given in Table 1. Only one parameter varies and other dimensions have been fixed.

It is well-known that the number of pole-pairs has a great impact on magnetic couplings performances. For a given air-gap length, there is always an optimum value. Fig. 13 shows that the analytical model is able to predict this optimal value which is equal to p = 4 for an air-gap b = 3 mm. In general, the optimal value for the number of polepairs decreases when the air-gap increases [START_REF] Lubin | Simple analytical expressions for the force and torque of axial magnetic couplings[END_REF].

It can be seen in Fig. 13 that the torque prediction obtained with the semi-analytical model is underestimated for large number of pole-pairs and slightly overestimated for Fig. 15. Torque versus the salient-pole thickness. Fig. 16. Flowchart of optimum design with GA p < 4. This is mainly due to the curvature effect that has been neglected in the proposed 2-D analytical model, while 3-D FE simulations take this effect into account [START_REF] Lubin | Improved 3-D analytical model for axial-flux eddy-current couplings with curvature effects[END_REF].

For a reluctance magnetic coupling, another important parameter is the slot-opening to pole pitch ratio β as shown in Fig. 2. Fig. 14 shows that the maximum torque occurs for a large value of β  0.8, which corresponds to a small opening angle for the salient-pole (20%) compared to the PMs pole-arc to pole pitch ratio (83%). This result is in agreement with a previously published paper about reluctance magnetic couplings [START_REF] Nagrial | Analysis and performance of variable reluctance (V. R.) torque coupler[END_REF]. When the salient-pole opening becomes very small (i.e. β > 0.9), magnetic saturation occurs and the torque decreases more rapidly than the one predicted by the semi-analytical model. This loss of accuracy for the analytical model will not be a problem for the optimization procedure because this area (very low value for the salient-pole opening angle) never corresponds to an optimal value. Fig. 15 depicts the torque versus the salient-pole thickness. It can be observed that the torque increases with the salient-pole thickness and becomes constant for c ≥ 15 mm. This result indicates that it is important to choose this value carefully to obtain the best performance (not too small and not too big). Once again, very good prediction is obtained with the semi-analytical model.

Design optimization and experimental result

Design optimization procedure with GA

As the semi-analytical model is very fast in terms of computational time and very efficient in terms of precision, it will be used in a design optimization procedure. It is worth noting that 50 ms is needed to compute the torque with the analytical model whereas it takes over 40 s with 3-D FE simulations, which corresponds to a time saving of more than 800. It is well-known that Genetic Algorithm (GA) is a robust and effective tool in optimization problems. Here the multi-objective genetic algorithm NSGA-II, which is available in Matlab has been used. This algorithm will be coupled with the proposed semi-analytical model.

The flowchart for the optimization procedure is given in Fig. 16. Seven parameters (R 1 /R 2 , R 2 , p, a, c, , and ) need to optimize. Two objective functions that correspond to minimize the magnet volume (V m ) and the external radius (R 2 ) have been defined. Also, one constraint is imposed which corresponds to the desired pull-out torque (16 Nm) while the air-gap length is fixed to b = 3 mm. The optimization problem can be summarized as given in [START_REF] Wu | Analysis and design optimisation of magnetic couplings using 3D finite element modelling[END_REF]. It should be noted that an additional constraint could be added in [START_REF] Wu | Analysis and design optimisation of magnetic couplings using 3D finite element modelling[END_REF] to limit the maximum flux density in the salient-pole to avoid magnetic saturation effects (B teeth ≤ 1.6T for example). Although the flux density in the salientpole is not calculated with the proposed semi-analytical model (infinite permeability), it can be estimated from the zcomponent of the flux density in the air-gap (8) by applying the magnetic flux conservation. This issue has not been addressed in this paper.

     
The bounds of the variables in vector x (16) were chosen as follows: R 1 /R 2 = 0.1-0.9; R 2 = 20-100 mm; p = 2-10; a = 2-20 mm; c = 2-25 mm;  = 0.1-0.9; β = 0.1-0.9.

The optimization procedure uses 100 individuals evolving for 50 generations. It takes less than 30 s to obtain the Pareto front given in Fig. 17 It can be observed that the optimal value for the volume of the magnets and therefore the price of the magnetic coupling increases when the outer radius of the magnets decreases and vice versa. The solution noted A in Fig. 17 corresponds to the optimal parameters given in Table 2. These geometrical parameters are those corresponding to the prototype built for the experimental validations and presented in the following section.

Design optimization procedure with PSO

Particle Swarm Optimization (PSO) is a stochastic optimization method for non-linear functions. It is based on the reproduction of a social behavior. The origin of this method comes from observations made during computer simulations of flights in a group of birds and schools of fish [START_REF] Kennedy | Particle swarm optimization[END_REF]. These simulations highlighted the ability of individuals in a moving group to maintain an optimal distance from each other and to follow a global movement in relation to local movements in their neighborhood. This social behavior based on the analysis of the environment and the neighborhood is then a method of optimum search by observing the trends of neighboring individuals. Each individual tries to optimize his chances by following a tendency that he moderates by his own experiences. PSO has been largely used for the design optimization of electrical machines and actuators.

A multi-objective PSO technique is used for the design optimization of the studied magnetic coupling. The objective functions and the constraint are the same as the one defined in [START_REF] Wu | Analysis and design optimisation of magnetic couplings using 3D finite element modelling[END_REF]. The resulting Pareto front is given in Fig. 18. The obtained results from this method are close to those obtained with GA (Fig. 17). The magnets volume for an external radius R 2 = 6 cm is slightly higher (73 cm 3 ) than that obtained with the GA method (70 cm 3 ).

Static torque measurement

To validate the semi-analytical model from an experimental point of view and to check the performance of the magnetic coupling, we have manufactured the prototype shown in Fig. 19. Fig. 19(a) shows the salient-pole side which is composed of ten ferromagnetic teeth. The thickness of the back iron has been chosen to avoid magnetic saturation. We used low carbon steel (grade AISI-1010) with high relative permeability in the linear region of the B-H curve and a high level of saturation to improve the performance. The second disc is made using sector-type NdFeB glued on the iron yoke as shown in Fig. 19(b). The geometrical parameters are given in Table 1.

Fig. 20 shows the reluctance magnetic coupling placed on the test bench. The axial coupling is inserted between two electrical machines (DC motors, 3 kW, 1500 rpm). In fig. 20, the air-gap value is b = 3 mm. The air-gap length has been set by inserting non-magnetic plates of known thickness between the two discs. Figs. 20(a) and 20(b) respectively show a photograph and a block-scheme representation of the test bench arrangement for the static torque measurement. As shown in Fig. 20(b), the static torque was measured thanks to weights (250 g, 500 g, 1 kg) suspended to a rod (l = 1 m) locked to one rotor, the other being fixed. The relative angular position δ was measured using an incremental encoder with a resolution of 4096 pulses/revolution (precision of 0.088 degrees) and the data was transferred into a computer.

The measured static torque as a function of load angle δ is shown in Fig. 21 and is compared with the ones obtained with the semi-analytical model.

Three values for the air-gap were considered (b = 3 mm; 5 mm, and 7 mm). As shown in Fig. 21, good agreement between the measured and predicted torque values have been obtained, which clearly show the validity of the proposed analytical model. The error on the pull-out torque prediction is never greater than 10%. The small differences between the model and the tests can be explained by the test procedure to measure the load torque, which is relatively simple (weight suspended at the end of a long bar) and may lead to measurement errors.

Limitations of the semi-analytical model

In the proposed semi-analytical model, the curvature effect and the magnetic saturation of the iron parts have been neglected. Although we have obtained accurate results in sections 3 and 4 for the studied magnetic coupling, some lack of precision may occur for particular geometrical configurations of the magnetic coupling such as very low air-gap and tooth opening values (saturation effect), or for very low or very large values of the pole-pairs number where the curvature effect is more pronounced.

Curvature effect

As shown in [START_REF] Dolisy | A new analytical torque formula for axial field permanent magnets coupling[END_REF] and [START_REF] De La Barrière | 3-D formal resolution of Maxwell equations for the computation of the no-load flux in an axial flux permanent-magnet synchronous machine[END_REF], the curvature effects can relative error increases but remains below 10%, which is acceptable. This result confirms the ones obtained in [START_REF] Dolisy | A new analytical torque formula for axial field permanent magnets coupling[END_REF] and [START_REF] De La Barrière | 3-D formal resolution of Maxwell equations for the computation of the no-load flux in an axial flux permanent-magnet synchronous machine[END_REF] where it was shown that curvature effect can be neglected without great errors for axial-field actuators.

Magnetic saturation effect

As indicated previously, magnetic saturation is the most important assumption of the proposed model. Magnetic saturation occurs in the teeth and is highly dependent on the air-gap and tooth opening value, the smaller the air-gap or tooth opening, resulting in the greater the magnetic saturation effect.

In order to highlight the magnetic saturation effect on the pull-out torque, it has been computed for two values of the air-gap length (1 mm and 3 mm) and by varying the slot opening , as shown in Fig. 23. Fig. 23 shows that magnetic saturation has a great impact on the pull-out torque for very low air-gap and tooth opening dimensions. It can be seen that the pull-out torque is well predicted until  < 0.8 for an air-gap of 3 mm whereas this value drops to  < 0.7 for an air-gap of 1 mm. Fig. 23 shows that the optimal dimensions for the salient-pole correspond to the limit of the magnetic saturation (knee-point of the BH curve). It should be noted that an additional constraint could be added in the optimization procedure (Fig. 16) to limit the maximum flux density in the salient-pole to avoid magnetic saturation effects (B teeth ≤ 1.6T), then infinite permeability assumption will be valid.

Conclusion

Reluctance magnetic coupling can be a solution for torque transmission between two shafts without contact. Compared to more classical synchronous magnetic couplings with permanent magnets on the two discs, it can be used in high temperature or corrosive environments. In this paper, we have developed a new analytical tool for the design optimization of a reluctance magnetic coupling. A correction factor has been proposed to take into account the 3-D effects in the radial direction. Comparisons with 3-D finite element simulations and measurements have shown the ability of the semi-analytical model in predicting the pull-out torque and the axial force with precision. Moreover, as the semi-analytical model is very effective in terms of computational time, it has been used in a design optimization procedure with multi-objective genetic algorithm. Finally, we have shown the limits of the model in terms of magnetic saturation.

Nomenclature

R 1

Inner radius of the magnets.

Appendix

The 10 integration constants are obtained by solving a linear system of 10 equations by using Matlab:

   

Fig. 3 .

 3 Fig. 3. Radial dependence of the flux density in the air-gap [4].

Fig. 4 . 6 .Fig. 7 .Fig. 5 .

 4675 Fig. 4. 2-D model in the x-z plane at y = 0 (under a salient-pole) Fig. 6. Impact of the air-gap length on the correction factor.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Flux line distribution under no-load condition obtained with 3-D FE simulation

Fig. 8 .

 8 Fig. 8. 3-D FE simulations: mesh distribution in the magnet and in the iron parts.

Fig. 11 .

 11 Fig. 11. Pull-out torque versus the air-gap length with or without the correction factor.

Fig. 12 .

 12 Fig. 12. Axial force versus the air-gap length with or without correction factor.

Fig. 13 .

 13 Fig. 13. Torque versus the number of pole-pairs.

  .

Fig. 17 .Fig. 19 .

 1719 Fig. 17. Pareto front obtained with GA: Magnets volume vs. Outer radius R 2 for T max = 16 Nm.

Table 2 :Fig. 18 .

 218 Fig. 18. Pareto front obtained with PSO: Magnets volume vs. Outer radius R 2 for T max = 16 Nm.

Fig. 21 .Fig. 20 .

 2120 Fig. 21. Measured and predicted static torque for three air-gap values.

  be analyzed by considering a dimensionless number  defines as the ratio of the radial excursion of the magnets H to the pole-pitch  : for  means that the curvature is pronounced. In order to change the value of  given in[START_REF] Dai | Analytical Modeling of Axial-Flux Permanent Magnet Eddy Current Couplings with a Slotted Conductor Topology[END_REF], the pole-pitch value is varied by changing the pole-pairs number from p = 1 to p = 10, which gives 0.42 <  < 4.24. The other geometrical parameters are kept constant and are those given in TableIwith an air-gap value fixed to 3 mm. Fig.22gives the relative error on the pull-out torque between the semi-analytical model and the 3-D FE simulations. It can be seen that the pull-out torque is predicted without error if the curvature coefficient  ≃ 1.8.For lower or greater values of the curvature coefficient, the

Fig. 22 .

 22 Fig. 22. Curvature effect: relative error on the pull-out torque between 3-D FEM and the semi-analytical model.

Fig. 23 .

 23 Fig. 23. Magnetic saturation effect: pull-out torque versus the slot opening for two air-gap values: 1 mm and 3 mm.

R 2

 2 Outer radius of the magnets. R m Mean radius of the magnets (R m = (R 1 +R 2 )/2). H Radial length of the magnets (H = R 2 -R 1 ). a Magnets thickness. b Air-gap length. c Salient-poles thickness. α Permanent magnets pole-arc to pole-pitch ratio β Slot-opening to pole-pitch ratio p Pole-pairs number. B r Remanence of the magnets. τ Pole pitch (τ = R m π/p).

Table 1 :

 1 Geometrical parameters of the studied reluctance magnetic coupling

	Symbol	Quantity	value
	R1	Inner radius of the magnets	30 mm
	R2	Outer radius of the magnets	60 mm
	a	Magnets thickness	10 mm
	b	Air-gap length	variable
	c	Salient-pole thickness	15 mm
	α	PMs pole-arc to pole-pitch ratio	0.83
	β	Slot-opening to pole-pitch ratio	0.77
	p	Pole-pairs number	5
	Br	Remanence of the magnets	1.25 T


G Boudebs

V Besse

H Leblond

Saturation of the Kerr effect in carbon disulfide

Accurate measurements of the nonlinear index in carbon disulfide shows an unusual saturation law of the Kerr effect which can be used to built a model of light propagation describing the filamentation in this liquid.

Introduction

Filamentation of a high intensity optical pulse was observed even before the concept of optical soliton was introduced [1]. The high intensities reached in filaments allow observing many nonlinear phenomena, such as high harmonic or supercontinuum generation [START_REF] Couairon | Femtosecond filamentation in transparent media[END_REF]… Recently two-dimensional spatial soliton propagation essentially stabilized due to χ (3) χ (5) balance in carbon disulfide (CS 2 ) at λ = 920 nm was reported [3], naturally paving the way to light filamentation in this liquid. However, at both wavelengths 1064 and 532 nm the standard model for the description of nonlinear dynamics of a light beam in this medium, which is based on the cubic nonlinear Schrodinger (NLS) equation, eventually refined by adding several terms to account for higher order nonlinear coefficients, does not remain valid as an accurate description of the saturation of the Kerr effect. We report in this paper on accurate measurements of effective n 2 really encountered by the material versus increasing intensities. These measurements lead to an alternative model in filamentation when combined to NLS equation [START_REF] Besse | Filamentation of light in carbon disulfide[END_REF].

Experiments: observation of the filamentation and nonlinear coefficient measurements

First we investigate the propagation of a laser beam in a cell filled with CS 2 , with thickness L = 10mm. The experimental setup is composed of a focusing lens (f = 20cm) into the entry of the cell and an imaging system (x10 magnification) that can be translated in order to image the entry and the output of the cell. The light source is an Nd: YAG laser, emitting at wavelength λ = 532 nm, 12 ps long pulses with Gaussian time profile, linearly polarized, operating at a repetition rate of 10 Hz. The output beam is acquired using a CCD sensor triggered with the laser. Beam shaping is achieved using a spatial filter, so that the spatial profile of the collimated beam is close to Gaussian at the input of the setup. By adjusting the beam waist at the focus ( 0 R z 13µm ω = λ π = ), we fixed z R , the Rayleigh distance to be one tenth of L. Above an input energy per pulse of 1.4 µJ, we observed the formation of filaments at the output of the cell. Changing the input intensity allows to increase the number of filaments (Fig 1-a). To give a theoretical account of these observations, an accurate measurement of the nonlinear parameters of CS 2 is required. We performed it using the well-known Z-scan method [5]. The standard model for the evolution of the light intensity I and the beam phase ϕ during propagation in a nonlinear medium, in the thin sample approximation can be set as:

2 3 dI αI βI γI dz = -- -, (1)
( )

2 2,eff 2 4 d kn I k n I n I dz = = + ϕ , ( 2
)
where z is the propagation direction, k the wave vector. α, β and γ are the linear, two-photon, and three-photon absorption coefficients, n 2 and n 4 the third-order and fifth-order nonlinear refraction indices, respectively. Measurements using a nonlinear fitting of the curve yielded by numerical solution of the differential evolution equation and analytical solutions of equations ( 1) and ( 2), with the necessary approximations, have been given [6]. Eqs. ( 1) and ( 2) can be solved analytically in special cases where some of the nonlinear coefficients vanish [7]. In the general case, an analytical expression of z versus I is found and inverted numerically [8]. It is known that CS 2 has neither linear nor two-photon absorption (α = 0, β = 0). Moreover, one of the absorption bands of CS 2 is centred at about 200 nm (see figure 1(a) in [3]), which is approximately the third harmonic λ/3 of the wavelength λ = 532 nm. This closeness to resonance increases the fifth order optical susceptibility, ( ) 5 χ ∝ γ . We used the recent D4σ-4f-Z-scan method [9], together with the numerical inversion related to the analytic solution of (1) and (2) (see [8] for details). The change in transmittance of the material is observed via a cooled CCD camera (-30 • C). For this, the used setup is a 4f system. A top-hat beam object is located in the front focal plane focus of the first lens, its Fourier transform is obtained in its rear focal plane, where is placed the sample, here a 1 mm thick cell filled with CS 2 . The same point is also the object plan of the second lens, which performs another Fourier transform into the CCD. A recent improvement of the technique allows to avoid the division of the two normalized Z-scan transmittances (closed and open aperture ones) by recording the beam waist relative variations according to the second order momentum. We measured the effective third-order refraction index n 2,eff with increasing intensity, for two wavelengths λ = 532 nm and λ = 1064 nm (see Fig. 1-b). We observe a steep increase of n 2 as the input central peak intensity I 0 increases. Then, after having reached a maximum value, n 2 decreases rapidly again. We observed the same evolution for both wavelengths, except that, at λ = 1064 nm, the variations are slower and the intensity range is wider. These results call into question the use of CS 2 as a reference material for calibrating nonlinearities measurement systems since it is not a constant even in the weakly nonlinear regime. The definition of n 4 is a fortiori problematic.

Results and conclusions

The measurements show that the usual formula 2,eff In other words, n 4 is positive at the considered wavelengths, and not negative as in [3]. Indeed, measurements performed at moderate intensities have provided the values n 4 = (1.2 ± 0.3)×10 -32 m 4 /W 2 at λ = 532 nm, n 4 = (2.2 ± 0.4)×10 -33 m 4 /W 2 at λ = 1.064 µm [8], in contrast with n 4 = -2×10 -35 m 4 /W 2 at 800 nm [10], or n 4 = -5.2×10 -35 m 4 /W 2 at 920 nm [3]. The change in the sign of n 4 with wavelength is related to the location of the third harmonic with respect to the absorption bands of the material. The main point is that, in the present situation, the saturation of the Kerr effect cannot be accounted for by the n 4 term. One may consider the option of retaining more terms in the expansion of the response function of the material in a power series of the electric field (higher order nonlinear susceptibilities). However, a difficulty would arise then. From theoretical point of view, the expansion cannot remain valid when the higher order corrections get the same order of magnitude as the leading one. In [3], it was possible to have the n 4 term comparable to the n 2 one, because the third harmonic resonance increased considerably the former with respect to the latter, which was not resonant. Hence we decided to provide an empirical model which describes the experimentally observed variations of n 2 with intensity, by fitting the experimental curve with an adequate Ansatz. After several trials, we found that

( ) ( ) 2 2
I aI 1 b I = +

2,eff n allowed to correctly reproduce the data at 532 nm and 1064 nm. The best fit, for data at 532 nm, is obtained with the set of parameters a = 3.4×10 -32 m 4 /W 2 , b = 5.8×10 -15 m 2 /W; while at 1064 nm, a = 6.3×10 -33 m 4 /W 2 and b = 2.3×10 -15 m 2 /W as illustrated in Fig. 1-b. Numerical computations based on a NLS-type model using the above expression of the nonlinear index allowed to reproduce the experimental observations of the filamentation [START_REF] Besse | Filamentation of light in carbon disulfide[END_REF]. In conclusion accurate measurements of the nonlinear optical parameters have shown unexpected variations of the effective nonlinear index with intensity in CS 2 . The saturation of the Kerr effect, which is responsible for the formation of filaments, cannot be using the standard n 2 -n 4 model leading to an alternative model in filamentation when combined to NLS equation.

Fig. 1 :

 1 Fig. 1: (a) Filamentation observed at the output of the CS2 cell, for value of the pulse energy 14.2 µJ, (b) 2,eff n versus 0 I at 532 nm (solid red line) and at 1064 nm (solid blue line). The dashed lines correspond to fitting taking into account the uncertainties.

  be used any more. The situation strongly differs from the experiments of[3] because the nonlinear index 2,eff n is here increasing with moderate intensity instead of being decreasing.


G Boudebs

V Besse

C Cassagne

H Leblond

Cid B De Araújo

D4σ σ σ σ nonlinear measurement inside a 4f-Z-scan system

The D4σ method using a 4f-Z-scan system is presented in order to increase the sensitivity and the resolution in Z-scan measurements. The numerical calculations and the experimental results validate our approach.

The inspection of linear and nonlinear (NL) optical parameters such as refractive index and absorption coefficient is essential for most of the applications in optics (optical limiting, lasers, optical amplifiers, etc.). Measurements of NL refraction and absorption using the Z-scan method [1] are widely performed nowadays. A variant of this method called Baryscan [2] reporting an increased phase distortion resolution (PDR) ( 50 000 λ

) has been published using stable CW laser input. More recently van Stryland's group reported a dual-arm Z-scan technique [3] improving the PDR up to 1 000 λ using pulsed lasers. These recent studies show that improvement of the PDR is still relevant. Measuring the diffraction efficiency allowed us to compare the sensitivity of different techniques inside a Z-scan 4f-system [4]. But still remain open questions to fully understand the physical phenomenon contributing for sensitivity and PDR improvements. We demonstrate here that the sensitivity of Baryscan is twice lower than that of Z-scan and we show that the use of hard-(physical) apertures (as a razor blade) with pulsed laser presenting pointing fluctuations could be a severe limitation leading to both lower sensitivity and small S/N ratio when compared to the direct measurement of the output beam waist variation using a CCD sensor.

The NL image formation inside the 4f -system shown in fig. 1 (a) is described using a model based on Fourier optics. The general scheme of beam propagation inside the 4f-system is described in details in [4]. We assume that the electric field at the object plane is Gaussian. We propagate the beam using the Helmoltz wave equation and the phase transformation related to lens thickness variation up to the image plane, taking into account the NL response of the material in the focus. We assume cubic nonlinearity and a thin NL medium of thickness L exhibiting (i) linear absorption defined by α (m -1 ), (ii) two-photon absorption defined by β (m/W) and (iii) NL refraction defined by

eff 0 0 0 0 I I log 1 q q = +
, the effective intensity. The beam waist measurement is performed using the ISO standard definition [START_REF]ISO Standard 11146, Test methods for laser beam widths, divergence angles and beam propagation ratios[END_REF]. Based on the second moment of ( ) I x, y the D4σ method gives 4 times the standard deviation of the intensity distribution. The beam radius in the x direction is:

( )( ) ( ) 2 x 2 I x, y x x dxdy I x, y dxdy +∞ +∞ +∞ +∞ -∞ -∞ -∞ -∞ ω = - ∫ ∫ ∫ ∫
where ( ) ( )

x I x, y xdxdy I x, y dxdy

+∞ +∞ +∞ +∞ -∞ -∞ -∞ -∞ = ∫ ∫ ∫ ∫
is the centroid of the beam. For Gaussian beams, the D4σ The detailed description of our experimental setup is given in [4]. Excitation is provided by a Nd:YAG laser linearly polarized 17 ps single pulses at 1064 λ = µm. In the image plane we use a CCD camera performing Zscan, Baryscan and D4σ profiles out from the same acquired images at each z-position by numerically changing the soft aperture. When using the D4σ method one have to calculate the centroid of the beam inside the frames which is equivalent to follow the pointing instability of the pulsed laser and thus reduces the noise that could appear with hard apertures as the ones used with Z-scan and Baryscan.

Fig. 2 (a) shows the experimental results for As2Se3, a chalcogenide glass with relatively high NL absorption. Notice that the Baryscan signal is totally buried inside the noise due to laser pointing instability while the beam waist measurement method gives a higher signal to noise ratio. The same behaviour could be found for Z-scan using hard aperture. The calculated profiles (solid line for D4σ and dashed line for Baryscan) are shown in the same figure. The agreement is very good. Moreover the D4σ method is unambiguously twice more sensitive than the Baryscan method because physically, we consider the totality of the pixels in order to measure the beam waist. So the quantity of information is twice that obtained with the razor blade in Baryscan. This is illustrated in Fig. 2 (a) by the difference between the peaks and the valleys inside the two scans processed for the same acquisitions showing approximately a factor 2 in the sensitivity. Moreover the PDR in Baryscan using pulsed laser with hard aperture in our experiment is approximately λ/180 which is a factor 5 lower than the PDR obtained with the D4σ method.

To obtain the comparison between Z-scan and D4σ and to improve the signal-to-noise ratio in Z-scan the centroid calculation is used for both profiles in Fig. 2(b) to correct for pointing laser fluctuations. The processing have been made using the same acquired images showing unambiguously the same sensitivity. The radius of the soft circular closed aperture in Z-scan is calculated to maximize the optical diffraction efficiency [4]. However, the advantage in the D4σ method when compared to Z-scan is that there is no need to divide two different Z-scan traces in order to obtain the NL refractive response.

In summary, the D4σ method is insensitive to pointing instability of the pulsed laser because no hard aperture is employed as in the Z-scan or Baryscan methods. Numerical calculations allow to obtain simple relations that can be used for the measurements simplifying the procedure especially for NL absorbing material.

Fig. 1 .

 1 Fig. 1. (a) Schematic of the 4f-system. The sample (M) is moved along the focal region. The labels refer to: lenses (L 1 -L 2 ). (b) Relative variation of the beam waist versus z. The other parameters are eff NL0 0.8 ϕ = and 0 q 0.6 = ; (c) Calculated ∆ω pv , versus the effective phase shift.

2 n

 2 (m 2 /W). We define the on axis NL parameters at z = 0,

  . The calculated (ωNL-ωL)/ωL, the relative variation of the beam waist versus z, the position of the sample in the focus, is shown in Fig.1 (b), where ωNL and ωL denote respectively the mean values (measured along x and y) of the beam waist associated to the NL and linear profiles at the output of the 4f-system. The difference between the peak and the valley ∆ωpv shown in Fig.°1 (c), is a function of the effective phase shift at the focus for a Gaussian input beam showing a linear relation ( . Moreover the linearity remains valid in presence of relatively high NL absorption.

Fig. 2 .

 2 Fig. 2. (a) Comparison between D4σ (squares) with that of Baryscan (circles) profiles. The solid and dashed lines show the calculations. (b) Comparison of the beam waist relative variation (squares; vertically shifted to 1) and the usual normalized Z-scan transmittance (stars) for a highly NL absorbing material. The inset shows the open aperture Z-scan transmittance. The NL parameters obtained were: q0 = 0.6 (see the inset in Fig. 2b) and according to Fig. 1, eff NL0 = 0.2 ϕ