Antoine Bruyère

Thomas Henneron

Fabrice Locment

Eric Semail
email: eric.semail@lille.ensam.fr

Alain Bouscayrol

J M Dubus

J C Mipo

Identification of a 7-phase claw-pole starter-alternator for a micro-hybrid automotive application

This paper deals with the identification of a new high power starter-alternator system, using both: a Finite Element Method (FEM) modeling and an experimental vector control. The drive is composed of a synchronous 7-phase claw-pole machine supplied with a low voltage / high current Voltage Source Inverter (VSI). This structure needs specific approaches to plan its electrical and mechanical behaviors and to identify the parameters needed for control purpose. At first, a Finite Element Method (FEM) modeling of the machine is presented. It is used for the predetermination of the electromotive forces and of the torque. Experimental results are in good accordance with numerical results. In a second part, resistive and inductive parameters of the drive are determined by an original experimental approach that takes into account each component of the drive: the battery, the VSI and the machine.

I.

INTRODUCTION

Claw-pole synchronous machines with separate excitation are commonly used to achieve the alternator function in Automotive. This is due to their low cost and their ability to work in a very large speed range. With the studied belt driven starter-alternator system [1], starting the car Internal Combustion Engine (ICE) is made with this claw-pole machine. This kind of system is used to bring the micro-hybrid Stop-Start function, without major modifications of the car powertrain. The new needed starter function adds a new constraint on the electrical machine: the ability to develop a large torque during the start. Starter-alternators already equip small cars, with a small 1.4 ICE. The major interest of this simple system concerns the limited extra cost for the final car [1], when the hybridization operation is achieved.

To take up this challenge for powerful Internal Combustion Engine (ICE) without changing the classical (12 Volts) DC-Bus voltage level [1], a new claw-pole starter-alternator with seven phases and permanent magnets between the claws has been developed [2]-[3] (Fig. 1). The obtained torque density (Nm/m 3 ) allows keeping the economical benefits of using a 12-Volts battery [1].

In order to carry out the starter function, a modeling of the machine is at first necessary to plan the ability of the machine to develop the torque during the short duration of the start. In section II, the predeterminations obtained by a Finite Element Method (FEM) modeling are compared with experimental results for the electromotive forces and for the torque per Ampere.

In section III, resistive and inductive parameters of the whole drive are determined. These parameters are needed for controlling the starter-alternator, in its both functioning modes: in motor mode (for the ICE start) and in alternator mode. As the voltage is low (12V), the resistive and inductive electrical parameter values of the battery and VSI are not negligible in comparison with those of the machine. At first, results obtained with classical determination approach are given and discussed. Secondly, an experimental approach based on an elementary vector control of the machine is developed to find directly the six characteristic time constants needed for controlling the seven-phase drive.

II. DETERMINATION OF EMF AND TORQUE BY FINITE ELEMENT METHOD

This kind of claw-pole structure is known to be sensitive to modeling with a numerical method as FEM [4] because of their 3D characteristics, a thin airgap (about 0.3% of the external diameter) and highly saturated magnetic materials. Moreover, the studied machine is a synchronous 7-phase claw pole machine with permanent magnets between the claws. Sixteen poles and a number of slots per pole and per phase equal to 0.25 imply a modelling of a quarter of the structure at least (Fig. 2). The aim of the modeling is not only to plan the voltage output, in generator mode, but also to determine the maximum available torque in motor mode, during the start of the car ICE. The first step in order to validate the 3D-modeling is to compare the experimental electromotive forces (emf) with the calculated ones. The second step is to obtain, for a defined repartition of the currents among the phases, the maximum possible torque.

The thin airgap leads to a mesh composed of 380000 tetrahedral elements, the air-gap being meshed with two layers of elements. The resolution is made using the magnetostatic assumption and the scalar potential formulation, with CARMEL, a software developed at the L2EP. Fig. 2 gives results of the numerical simulation, related to the magnetic flux density B under no load condition and a 5A current I F supplying the excitation coil. The highest flux density values (in yellow) are obtained under the stator teeth.

Fig. 3 compares the electromotive force waveform named e(t), calculated under no load condition, at a rotation speed N = 1800rpm and I F = 5A, with an experimental measurement in the same conditions. These results validate the numerical model good accuracy, at the first order, for the determination of the voltage outputs.

Fig. 4 gives the calculated torque T as function of the rotor position θ, when the machine currents are imposed to be the same as in an experimental test (Square Wave Voltage supply), for which the torque is maximized. The maximal numerical torque value is T calculated = 62.4 Nm, to be compared with the experimental maximum measured torque: T measured = 66.5 Nm, that shows a difference of about 6 %.

For the determination of the inductive parameters, numerical FEM leads to too large durations of computation in the present case. It is due to the 3D characteristics of the machine and its thin air-gap (which implies a high number of 380000 elements). Moreover, the magnetic materials being highly saturated for this application, the simulation has to take into account the non-linear characteristics of the materials. In Fig. 5, the magnetic flux density is represented for a 2D slice of the whole structure (rotor and stator); z-coordinate of the slice is at the middle of the structure. This figure points out the magnetic saturation areas under no load conditions, with i F = 5 A. These saturated areas evolve with the rotor position and the currents (excitation coil current i F , and stator windings currents). They are at the origin of variable magnetic air-gap effects. As consequence, the inductances depend on the rotor position and on the currents. Finite Element modeling to calculate the inductive parameters for every rotor position, and for every supplying situation, would lead to huge computation and analysis times.

As consequence, experimental methods have been developed for the determination of the machine inductive parameters. For a 7-phase machine, three different dq-reference frames, or "subspaces S1, S2 and S3", can be determined. Consequently, six time constants can be then evaluated, depending on the importance of the reluctance effects. These time constants are needed for the control.

Two methods, each one based on current control, are presented and compared in the next paragraphs. The experimental set-up is described with Appendix A2. It shows the starter-alternator directly mechanically connected to a brushless machine, used to simulate the ICE behavior. The seven-leg VSI is supplied by an electronic programmable source. A programmable electronic load is used in parallel with the voltage source in order to simulate the on-board electrical network energy consumptions. The whole system is managed using a modular dSPACE DS1006 real time control board.

A. Determination of resistance and inductance in the stator frame

The first approach is based on a simplified model. Starting from the inductances definition (1):

[

Flux] = [L SS ] [I s ] (1)
with [Flux] the seven-dimensional vector of flux linked by the seven stator phases, [I s ] the associated current vector and [L SS ] the 7X7 stator inductance matrix.

Several assumptions are made. At first, the magnetic characteristic of the materials is considered to be linear. Secondly, whereas for a seven-phase machine the [L SS ] stator inductance matrix contains 49 values, only 7 different values, one self inductance and six mutual inductances, are considered taking into account standard assumptions of regularity, symmetry and circularity properties for this matrix. These seven inductances lead to six cyclic inductances and one homopolar inductance (or zero-sequence inductance). This last one is omitted since a wye-coupling of the machine implies that the homopolar inductance has no impact on the control.

The procedure of measurement is the following for different values of the rotor position, which is imposed and kept constant by a motion controller:

 a magnetic state is imposed by a choice of a value for the excitation current i F ;  sinusoidal current at chosen frequency is imposed in phase n°1 by using a VSI with a control of the currents;  voltage is measured across each phase successively;  measured voltage is integrated. Flux and corresponding inductance are then deduced using (1). The following measurements are made for several fixed positions. Then, Table 1 gives average values of inductances. Considering average values is in accordance with the circularity assumption used for this first method.

Practically, it has been chosen a 120 Hz frequency with a 17 A RMS sinusoidal current, around a constant 18 A current.

> i F =0A 45 -3 -7 -3 -3 -7 -3 i F =5A 39 -3 -6.5 -3 -3 -6.5 -3
With standard assumptions of symmetry and circularity for the [L SS (θ)] stator inductance matrix, the cyclic inductances in the three d-q reference frames of the seven-phase machine are obtained by (2). In (2), L 0 is the zero-sequence inductance. The six other terms are considered by pairs, in two dimensional subspaces S1, S2 and S3 defined with the generalized Concordia transformation (denoted with the matrix [C] in (2), cf. Appendix A1). Table 2 gives the inductive parameters once the generalized Concordia transformation has been achieved. Considering separately these three subspaces, there is no difference between the inductances along the d-and q-axes. This is due to the circularity assumption of the inductance matrix (variable airgap effects are not taken into account for this first method). Table II gives average cyclic inductance values versus . Finally, with the knowledge of these cyclic inductances it is easy to determine the time constants, dividing the inductance by the resistance. However, using the resistance value to calculate the time constants values introduces a new uncertainty as it is detailed below.

    1 3 3 2 2 1 1 0 ) ( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                                               C C L L L L L L L q S d S q S d S q S d S  SS L (2)
A first approximation of the phase resistances value (R S ), using the stator windings geometrical characteristics, gives R S = 7.6 m Measurement with a digital RLC meter bridge using a low current level (0.5 A) gives a confirmation for this winding resistance value: 8 m at low frequency. However, these values concern only the machine windings. Such low resistance value implies the necessity to take into account all the parasitic resistances usually neglected: resistance of the VSI MOS transistors and resistance at the electrical connections. This uncertainty will be confirmed by the final results of this paper, using the second approach, for which the total resistance value is about 21,7 m.

If this first method is theoretically easy to understand and to use, it must be also noticed that the values of inductances are weak (a few H, especially the mutual terms of [L SS ]). So, the induced voltages are only a few Volts. Consequently, parasitic voltage drops induce uncertainties on the measurements. Moreover, the use of Concordia [C] transformation matrix implies that each determined cyclic inductance is obtained by the sum of six ones: as consequence there is an accumulation of errors on the measurements with this method.

If it will be observed in the following paragraph that the results are acceptable for average values, it must be pointed out that with this approach it is not possible to take into account the effect of all the stator currents on the values of the inductances. Fundamentally the method is based on an assumption of magnetic linearity. Basically it is thus not possible to impose to the machine a magnetic state representative of real working conditions of the drive. Nevertheless, obtained values allow a first estimation of inductive parameters and time constants. With the next approach it is possible to determine directly the cyclic inductances in conditions closed to the working ones.

B. Direct determination of time constants in dq-frames

In the studied case, since the values of resistances and inductances are very low, it is then preferable to consider a direct method of identification that takes into account implicitly the complete drive and gives directly the time constants [10]. The experimental method is based on an elementary vector control for a seven-phase machine. The control is a generalization of classical vector control in dqreference frame developed for three-phase machines. In this case, two L d and L q parameters are obtained for a machine with reluctance effect (leading to 2 time constant values, on d-and q-axes). For a 7-phase machine, three different dq-reference frames, or "subspaces S1, S2 and S3", can be determined. Consequently, six time constants can be then evaluated. The principle of the control is given Fig. 6. In each subspace, there is a control of the dq-currents. The currents controllers are Proportional Integral controllers with compensation of the electromotive force. Constant current references are used except in the channel where the cyclic inductance must be identified. In this channel a simple Proportional controller is used with a square reference added to the constant reference. The PI controllers in the different channels allow ensuring a decoupling with the identified channel. The error and the closed loop time constant relative to the response to the square reference allow obtaining resistance and cyclic inductance. Fig. 7, gives experimental results of the time constant measurements. These parameters depend on the excitation current. For one of the 3 dq-reference frames, the S1-subspace, two different values of time constants, on d-and q-axes, can be clearly identified: the difference can be explained by a variable reluctance effect in the machine. For the two others dq-reference frames, only one inductive parameter appears. These results are compared in table III and table IV with those obtained in table II for two values of I F . For S2 and S3 subspaces the error is satisfying (max 26 %). For S1 subspace, it appears a difference between L d and L q values. Indeed, the variable reluctance effect was not taken into account by the previous method because only average inductances values versus  angle were considered. The new approach confirms that a reluctance effect is quite noticeable for the S1 subspace. If the average value <LS1d> obtained by the first method is compared to the average value (L S1-d + L S1-q )/2 obtained by the last approach, the relative errors are 12% for I F =0A and 27% for I F =0A. The results are coherent.

                                                                                                                                                                                                                      

Fig. 1 .Fig. 3 .

 13 Fig.1. Exploded view of the multiphase claw pole starter-alternator

  : Seven-phase starter-alternator (rear of the picture) mechanically connected to a brushless synchronous machine (front of the picture) (b): Seven-leg Voltage Source Inverter (c): Programmable Electrical source and load

Table I

 I

	gives the average values (in µH) versus  for two
	values of excitation currents.

TABLE I

 I INDUCTANCES VALUES (µH) IN THE STATOR FRAME<L 11 > <L 12 > <L 13 > <L 14 > <L 15 > <L 16 > <L 17


A Bruyere

E Semail

F Locment

A Bouscayrol

J M Dubus

J C Mipo

Identification of sensitive R-L parameters of a multiphase drive by a vector control

This paper focuses on an experimental method to determine the electric parameters of a seven-phase lowvoltage multiphase drive. The drive is a belt driven starteralternator for powerful cars with Hybrid Electrical Vehicles (HEV) functions. The resistive and inductive parameters are necessary to obtain the six characteristic time constants of the control modeling. Classical direct measurements lead to imprecise results because of very low values for the windings electric resistance (a few mΩ) and inductance (a few µH). Effects of the imprecision on the measurements are all the more important that time constants are obtained by a ratio of cyclic inductances by resistance, with cyclic inductances being a linear combination of seven measured inductances. The methodology for identification detailed in this paper is based on a stator current vector control, in a multi-reference frame. This methodology allows us to get directly these time constants. Numerous measurements allow the robustness of the method to be evaluated.

INTRODUCTION

Reducing petrol consumption and CO2 emission is now assumed to be one of the great new societal and economical issues. In this context, automotive suppliers as Valeo, offer new solutions to make cars less costly in terms of energy and less polluting. One of these solutions consists in adding electrical machines running with the classical Internal Combustion Engine (ICE) to make it work more efficiently; it is the principle of Hybrid Electric Vehicles (HEV) [START_REF] Wouk | Hybrids: then and now[END_REF]- [START_REF] Emadi | Vehicular Electric Power Systems[END_REF]- [START_REF] Ehsani | Modern Electric, Hybrid Electric, and Fuel Cell Vehicles[END_REF]. From a small ten years, several hybrid cars structures have appeared on the market. This paper deals with the cheaper one: the starteralternator system [START_REF] Cai | Comparison and review of electric machines for integrated starter alternator applications[END_REF]- [START_REF] Richard | Valeo StARS Technology: A competitive Solution for Hybridization[END_REF]- [START_REF] Cogo | The Valeo belt driven starter alternator reversible systems[END_REF]. This system is composed of a classical car generator, i.e. a claw pole synchronous machine, but with a Voltage Source Inverter (VSI) instead of the diode bridge; the claw-pole machine can then be used as motor. Starter-alternators already equip small cars, with a small ICE, as the Citroën C3, and are used to a unique function: the Stop-start function [START_REF] Richard | Valeo StARS Technology: A competitive Solution for Hybridization[END_REF]- [START_REF] Cogo | The Valeo belt driven starter alternator reversible systems[END_REF]. The major interest of this simple system concerns the limited extra cost for the final car. In order to extend the range of models equipped with a start-alternator or to extend the system capabilities, there is a need to increase its power. Previous works [START_REF] Locment | Comparison of 3-, 5-and 7-leg Voltage Source Inverters for low voltage applications[END_REF] show the interest of multiphase (with more than 3 phases) VSIs to improve the use of the DCbus voltage. Moreover, the use of multiphase machines also enables the increase of torque density by using new windings [START_REF] Dubus | Polyphase stator of a rotating electrical machine with claw-pole rotor and alternator or alternator starter comprising same[END_REF]. It is thus possible to extend the power range, keeping the standard low cost low voltage level of 14V.

With the final objective to explore the capabilities of this kind of system, this paper focuses on the first difficulty: the determination of the control modeling parameters such as resistance, inductances and time constants. In a seven-phase machine with usual voltages, the determination of these parameters is already sensitive [START_REF] Locment | Design and Study of a Multiphase Axial-flux machine[END_REF]. With the studied machine, this sensitivity is amplified because of the low values of resistance (a few mΩ) and inductances (a few μΗ) induced by the low voltage value of the DC-bus. The original methodology of identification presented in this paper lies on a vector control of the multiphase machine currents, in a multi-reference dq frame based on a generalized Concordia transformation [START_REF] Figueroa | Generalized transformations for polyphase phase-Modulation motors[END_REF]- [START_REF] Levi | Multiphase induction motor drives-A technology status review[END_REF]. In the first part of this paper, the 7-phase drive modeling in the multi-reference frames is introduced. In the second part, the methodology of identification is presented. In the last part, experimental results are discussed.

MODELING THE 7-PHASE DRIVE IN A CONTROL FRAME

I.

A.

B.

Presentation of the 7-phase drive

The drive: it is composed of a 7-phase wye-coupled claw-pole synchronous machine with separate excitation. The machine is supplied with a Voltage Source Inverter (VSI). MOSFETS transistors (MOS) are used for the VSI switches. In normal operating modes (with high currents level), voltage drops of MOS are supposed to be a linear function of the current. The DC-bus voltage V DC is imposed by a 12V battery.

The experimental set-up: an experimental set-up has been developed in order to test control algorithms. It is described in Fig. 2. In this, the starter-alternator is directly mechanically connected to a brushless machine, used to simulate the ICE behavior. An electronic voltage source and an electronic load are used to replace the 12V battery. The whole system is managed using a dSPACE TM DS1006 control board.

Classical modeling of the 7-phase synchronous machine in the stator frame

In the stator frame, the 7-phase synchronous machine is commonly modeled with the set of equations ( 1) and (2). In Fig. 2, a classical scheme is used to introduce the common notations: v k , as the phase "k" voltage, and i k this phase current. The position of the rotor is θ and the excitation field current is i F . In (1) and ( 2), R S is the phase resistance value, , the magnetic flux linked by the phase k. The matrix [ is the stator inductance matrix, whose terms depend on the position θ. . In this paper, we focus on R s and :

[ ]

) (θ SS L
R S measurement in the stator frame: this measurement is quite difficult in the stator frame, because of its very low value (a few mΩ). Indeed, at first, to get a sufficient precision on the measurement, it is necessary to supply the windings with a sufficient voltage values. This implies very high value of current to flow through the windings (several hundreds of amperes). Moreover, the windings resistance value is close to the parasitic resistances (connection and wires resistances, MOS resistance …). So, it is necessary to use a methodology which takes into account all these parasitic resistances.

Terms of [ measurement in the stator frame:

This measurement is the most difficult, firstly, because of the very low value of the inductive terms, secondly, because of the evolution of the matrix terms. Indeed, is a full matrix, whose each term varies with the position (θ) and the magnetic state, which depends on the currents

(i ] ] ) (θ SS L [ ) (θ

SS L

F and i S ) values. So, it is difficult to determine the variation laws of the 49 inductive terms.

Finally, in order to control the machine, the needed parameters are time constants, i.e. ratios between inductive parameters and resistance. Then, a low precision on measurement of both inductances and resistances can lead to a completely false estimation of these time constants. In the second part, an original methodology is presented, which allows a direct estimation of these time constants.

[ ] [ ] [ ]

dt d i R v T T S T φ + = (1) [ ] [ ] [ ] [ ][] [ ] T SR T T T S T i f M i dt d i R v + = + = ) ( L SS θ φ (2) with [ ] [ ] 7 6 5 4 3 2 1 v v v v v v v v = ; [ ] [ ] 7 6 5 4 3 2 1 i i i i i i i i = ; [ ] [ ] 7 6 5 4 3 2 1 ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ = ; [ ] [ ] 7 6 5 4 3 2 1 f f f f f f f f = ; ) 7 / ) 1 ( 2 ( π θ - - = k f f k 1 2 3 i 1 i F v 1 θ 4 5 6 7 1 2 3 i 1 i F v 1 θ 4 5 6 7

C. Modeling in a control frame

The methodology of the electrical parameters identification presented in this paper lies on a control of currents in the generalized Concordia frame [START_REF] Figueroa | Generalized transformations for polyphase phase-Modulation motors[END_REF]- [START_REF] Levi | Multiphase induction motor drives-A technology status review[END_REF]. This is particularly interesting because the dynamic control of the alternative machines is generally established once the problem ( 1)-( 2) has been transposed in this frame. Then, studying the special case of a 7-phase machine, the mathematical transformation leads to a projection of the problem in 3 orthogonal dq-subspaces, named "S1-subspace", "S2-subspace" and "S3-subspace" (or S1, S2 and S3) [START_REF] Figueroa | Generalized transformations for polyphase phase-Modulation motors[END_REF]- [START_REF] Levi | Multiphase induction motor drives-A technology status review[END_REF]. This transformation yields a decoupled set of voltage equations, expressed with (3) in the S1, S2 and S3 subspaces. The first equation is related to the 0-sequence component, in a 0-sequence subspace, which is ignored because the electrical wye-coupling of the machine implies an always null i 0 current. Then, the 6 other equations yields a very classical and simple scheme (Fig. 3.) for controlling each dq-current expressed, in the 3 subspaces S1, S2 and S3. In these equations, the L S. terms are the cyclic inductances. As an example, L S1-d is the cyclic inductance related to the d-axis in the S1-subspace. These cyclic inductances are a linear combination of seven elements of [ ]

) ( L SS θ . As consequence, if these cyclic inductances are determined directly from elements of [ ] ) (θ SS L
, there is an accumulation of uncertainties.

⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ + + = + + = + + = + + = + + = + + = + + = - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - q S q S S q S q S q S d S d S S d S d S d S q S q S S q S q S q S d S d S S d S d S d S q S q S S q S q S q S d S d S S d S d S d S S e i R dt i d L v e i R dt i d L v e i R dt i d L v e i R dt i d L v e i R dt i d L v e i R dt i d L v e i R dt i d L v 3 3 3 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 / ) ( / ) ( / ) ( / ) ( / ) ( / ) ( / ) ( (3)
+ -

C(s) + + + - s K τ + 1

e compensation e i reference i

Controller 1 st order system with e as a perturbation + -+ -C(s) C(s)

+ + + + + + - + - s K τ + 1 s K τ + 1

e compensation e i reference i

Controller 1 st order system with e as a perturbation In the first section, two ways of modeling the 7-phase synchronous machine have been introduced: in the stator frame ( 1)-( 2) and in the generalized Concordia frame (3). The second one is generally used for a dynamic control of the machine (torque or speed control). The common methodologies for identifying the electrical parameters generally consist in the measurement of these parameters in the stator frame, i.e. in a frame where voltages and currents are directly measurable, and then to make the projection of these parameters in the Concordia frame.

Here, the methodology consists in directly measuring the parameters in the generalized Concordia frame.

These parameters are determined while controlling the six independent currents on both d-and q-axes, in S1, S2 and S3. Fig. 3 gives the model and the scheme for controlling only one of these currents: the parameters to be determined are K (the static gain, with K = 1/R S ) and τ (the time constant of the first order system, with τ = L S /R S ). "e" is the electromotive force, considered to be a perturbation in the model is compensated in the control by a feedforward action. Finally, C(s) represents the controller.

In the dq-spaces (3), the currents are then constant values in steady-state operation. The Matlab TM Simulink TM scheme used to control the whole system is given on Fig. 4 and is directly implemented in the dSPACE TM real time control board system. In this figure, it appears: (a), the step form references which can be imposed to each dqcurrent in the three different dq-subspaces, (b), the measurement of the "real" (in the stator frame) seven phase currents and their projection in the dq-subspaces, (c), the control structure of each current, (d), the estimation of the electromotive forces, in order to make the compensation and to minimize the perturbations on the different dq-axis , (e), the "inverse" transformation from the three dq-subspaces to the stator frame, and (f) the "real" system on which the control orders of the VSI are transmitted.

In order to identify time constant and gain relative to one of the six d-q axes, the approach is the following: a) Firstly, constant references of current are chosen in the five other axes and PI controllers are tuned in order to minimize interaction between the axes; b) secondly a proportional controller is chosen for the studied axis with a step waveform reference for the current; c) finally, measurements of closed loop static error and closed loop time constant allows us, according to Fig3., to determine open-loop K gain and τ time constant. Fig. 5 gives an example of the dq-currents waveforms in S1, S2 and S3-subspaces, obtained when the square wave reference is imposed to the q-axis in S1. So, this test aims at determining R S and L S1-q . For this example, the average working point is specified as followed: I F = 3A and <i S1-q > = -200A. The square wave reference, with a 20A (Fig. 5-(S1)) magnitude, is imposed around the average working point. The measured current along the qaxis in S1 is shown in Fig. 5. Its waveform corresponds to a first order system response, which allows the use of the model described in Fig. 3. All other currents are null. In this example, the proportional gain used for controlling i S1- q is K p-S1-q = 0.02. The measured static gain is 9.4/20 = 0.47. The deduced value of the stator resistance is R s = 22.6mΩ. The direct measurement of the closed loop time constant is τ S1-q-closed-loop = 1.14ms. This implies a value for the open loop time constant τ S1-q = 2.15ms and for the inductance L S1-q = τ S1-q R S = 48.5 µH. A.

GLOBAL EXPERIMENTAL RESULTS

In this part, robustness of the method is firstly examined on one axis. Secondly, the uncertainty on the stator resistance is evaluated by using the results obtained on the six d-q axes. Finally, evolution of the six time constants with the excitation current is given. The rotation speed is always null and the square wave reference magnitude is always 20A, around an average current on the studied axis of 200A.

Repititiveness of the measurement

As a first test, we want to examine the repetitiveness of the measures. Fig. 6 shows the evolution of the resistance (measured along the q-axis in S1, Fig. 6-(a)) and of the dand q-time constants in S1 (Fig. 6-(b)). For each measurement, two curves have been plotted: "meas1" and "meas2", in exactly the same conditions. Moreover, for this test, the proportional gain used for the controller on the studied axes (S1-q axis for Fig. 6-(a); S1-d and S1-q axes for Fig. 6-(b)) is the same gain as in the test described in Fig. 5, i.e. K p = 0.02. These measurements have been done for different excitation currents.

Finally, as the described identification is supposed to be independent of K p , tests have been done for different values of K p . So, for two excitation currents, i F = 3A and i F = 7.5A, the measurement has been made for: K p = 0.01, K p = 0.02 and K p = 0.04. For each case, two measurements have been made. This means that for these two excitation values, six measurements are compared Tau-S1d-meas1 Tau-S1d-meas2 Tau-S1d-average-meas Tau-S1q-meas1 Tau-S1q-meas2 Tau-S1q-average-meas B. Measurement of the resistance According to the model in the generalized Concordia frame (3), there is just one defined resistance R S . Indeed, the seven phases are supposed to be identical (from the electrical point of view). There are 6 ways to measure R S (two measurements on each d-or q-axis, in three subspaces S1, S2 or S3). The same value of resistance is supposed to be found for each measurement, in S1-d/q, S2-d/q and S3-d/q. It must be remarked that the measured resistance takes into account the equivalent resistance of MOSFET transistor of the VSI and all the parasitic resistances of wires and connections. Fig. 7 shows the evolution of the resistance as a function of the excitation current i F : the resistance appears to be quite constant while i F is changing. The global mean value is R S = 21.74mΩ. Table 1 gives the "partial" (related to each axis) mean values of the resistance, for each measurement. It also gives the relative error (in %) of these partial mean values relating to the average value. The maximal error, 10.7%, is given for the measurement on S3-q axis. A model using a single constant resistance R S appears then correct.

Axis: S1-d S1-q S2-d S2-q S3-d S3

C. Measurement of the time constants

Fig. 8 gives the evolution of the time constants measured in the three subspaces S1, S2, S3 (on d-and qaxes), as a function of i F . These time constants (and so, the inductances) decrease while the excitation current increases. This is due to the evolution of the magnetic materials properties, which proves their non-linear characteristic. This magnetic state depends on the currents (phase currents and excitation current). Here, the characteristic is just given as a function of the excitation current at a specific operating point. The evolution of the time constant with the stator currents evolution is not given. Nevertheless, the same approach has been lead at other operating points, and at different rotation speeds. The same measurements can be established as described in this paper without any major difficulty. The presented results, at a null rotation speed, correspond to the easiest measurement, because of a good decoupling between the six axes.

Finally, another point must be observed with the results shown at this operating point, particularly in S1: the inductance on the d-axis is much different with the inductance on the q-axis. This is due to variable reluctance effects. For this machine, at this operating point, the reluctant effects affect mainly the S1 subspace. For the other S2 and S3 subspaces, the d-and q-inductive parameters are almost the same. Tau-S1d

Tau-S1q

Tau-S2d

Tau-S2q

Tau-S3d

Tau-S3q CONCLUSION In order to establish the control of a 7-phase starter alternator drive and to build virtual models, a good knowledge of its electrical parameters and time constants is needed. Due to the low voltage and the high number of phases of the studied drive, indirect determination of the time constants, using classical inductive measurements in the stator frame, implies uncertainties.

Magnitude (mΩ)

In this paper, an experimental methodology for identifying directly in the Concordia frame, in which are determined the cyclic inductances, has been presented. Numerous measurements have been carried out to prove its robustness. Moreover, this methodology allows us to take into account all parasitic resistances (of the MOSFET transistors and at electrical connections), which are not negligible in this low voltage automotive application: it is really the identification of the drive and not of the electrical machine that has been achieved. The results have been given for a specific operating point, at null rotation speed. Nevertheless, the developed methodology has been also used at other operating points. Concerning more specifically the studied drive, the variations of the time constants with the excitation currents have been given, with a noticeable difference between the three characteristic subspaces S1, S2 and S3: only S1 subspace presents two different L d and L q cyclic inductances.

Figure 1 .

 1 Figure 1. Experimental bench picture In the stator frame, electrical parameters to be determined are then: the stator resistance R S , the electromotive forces and the terms of the matrix [ ] ) (θ SS L. In this paper, we focus on R s and :

Figure 2 .

 2 Figure 2. Schematic representation of the 7-phase synchronous machine with separate excitation

Figure 3 .

 3 Control of one single axis of a dq subspace II. IDENTIFICATION OF THE ELECTRICAL PARAMETERS IN CONTROL SUBSPACES

Figure 4 . 7 -

 47 Figure 4. 7-phase synchronous machine control structure; Matlab TM Simulink TM representation for dSPACE TM implementation III.

S1Figure 5 .

 5 Figure 5. Example of experimental results for a specific operating point: N = 0rpm, i F = 3A;Step reference is imposed to the q-axis in S1subspace and null references in the other axes

Figure 6 .

 6 Figure 6. Evolution of the resistance in S1-d (a) and the time constants in S1 (b) as functions of the excitation current; and comparison of different values of K p (circled zones)

Figure 7 .

 7 Figure 7. Evolution of the measured resistance in S1, S2 and S3, as a function of the excitation current

Figure 8 .

 8 Figure 8. Evolution of the time constants measured in S1, S2 and S3, as a function of the excitation current

This work was supported by the French car supplier Valeo and the regional council of France Region-Nord-Pas-De-Calais.


L Drouen

F Hauville

J F Charpentier

E Semail
email: eric.semail@lille.ensam.fr

S Clénet

A coupled electromagnetic / hydrodynamic model for the design of an integrated rimdriven naval propulsion system

This paper presents an analytical multi-physic modeling tool for the design optimization of a new kind of naval propulsion system. This innovative technology consists in an electrical permanent magnet motor that is integrated into a duct and surrounds a propeller. Compared with more conventional systems such as pods, the electrical machine and the propeller have the same diameter. Thus, their geometries, in addition to speed and torque, are closely related and a multidisciplinary design approach is relevant. Two disciplines are considered in this analytical model: electromagnetism and hydrodynamics. An example of systematic design for a typical application (a rim-driven thruster for a patrol boat) is then presented for a set of different design objectives (efficiency, mass, etc). The effects of each model are commented.

I. INTRODUCTION

The rim-driven permanent magnet (PM) propulsion system is a novel and emerging technology for the vessels propulsion. It consists in a synchronous PM machine that surrounds a propeller and is integrated into a duct. The permanent magnets are stuck on a soft magnetic material ring surrounding directly the blades. This assembly constitutes the rotor of the propulsion PM motor. The stator of this motor is inserted into the duct of the propeller. With this configuration, the gap between the rotor and the stator can be immerged in the sea water. In this case, active parts (windings, magnets, magnetic cores) are insulated from the sea water thanks to an epoxy resin. Compared with more traditional electrical propulsion system, it presents some interests such as a better hydrodynamic efficiency, the blades protection, a smaller electrical motor and the possibility to increase its rated power above the limits of the traditional pod thrusters [START_REF] Lea | Scale model testing of a commercial rim-driven propulsor pod[END_REF]. This kind of solution is now technologically mature and experimental studies, from industrial or academic laboratories, have already been performed in the last decade [START_REF] Abu-Sharkh | Structurally integrated brushless PM motor for miniature propeller thrusters[END_REF], [START_REF] Krøvel | Design of an integrated 100kW Permanent Magnet Synchronous Machine in a Prototype Thruster for Ship Propulsion[END_REF]. It can also be used in marine current energy harnessing [START_REF] Drouen | Study of an innovative electrical machine fitted to marine current turbines[END_REF] as shown in fig. 1. However, few multi-physic models for the design of those specific systems have been presented for the moment.

With this particular technology, the propeller and the electrical machine have the same diameter, torque and speed. For this reason, a coupled multi-physic design model is proposed in order to avoid a sequential approach that would be less relevant and imply time consuming calculus. To be inserted in a systematic design process, it seems necessary to develop a multi physic model which is accurate and simple enough to minimize the calculation time. In addition, the results given by the model must be insensitive to any mesh variation (due to a geometrical variation, as in numerical models). This is the reason why an analytical approach has been chosen for this work. This kind of model allows a fast and good convergence of such systematic design process. In this paper, a separate description of two models is given. The first one is an analytical first order electromagnetic (EM) model that is particularly relevant for this specific structure of electrical machine (section II). The second one is a model of propeller that is well known in the field of propeller design (section III). It is based on a set of typical ducted propeller data obtained from tests in ship model basins. The accuracy of each model is evaluated and both models are coupled (section III). Finally, a systematic design using the coupled model is achieved for a patrol boat propeller (section IV). The influence of each sub-system on the choice of the overall characteristics is discussed

II. ELECTROMAGNETIC AND THERMAL MODEL

The electrical model is used to deduce the dimensions and performances of the electrical machine for a given set of specification. It is a first order analytical model that permits a fast but fairly precise calculation of the characteristics of the machine and is adapted to an optimization work. It is directly inspired by equivalent analytical models that can be found in [START_REF] Polinder | Comparison of direct-drive and geared generator concepts for wind turbines[END_REF] and [START_REF] Grauers | Design of direct-driven permanent-magnet generators for wind turbines[END_REF] for instance. The input parameters of the model are the machine inner diameter D int (m) and rotational speed  (rad/s) as well as the propeller's torque Q (N.m). This model contains a certain number of variables that are fixed to relevant values by the designer. They may also vary in a given range if considered as key variables (such as current density, electric load, flux density, number of poles, etc).

The electrical machine is a synchronous PM radial flux machine. It is connected to an AC/DC Pulse Width Modulation voltage converter that can control the current wave into the stator windings. If the electrical motor is fed by a sinusoïdal current (with an appropriate control strategy), the medium EM torque is expressed as follow

  cos ). 4 / . . .( . . . 2 2 1 1 L D B A k T L b EM  (1)
where k b1 is the winding factor, A L (A/m) is the stator rms electric load, B 1 (T) is the peak value of the fundamental of the magnets flux density at the stator surface, D (m) is the gap diameter, L (m) is the iron axial length and  is the angle between the stator current and the electromotive force induced by the rotor.

A linear relationship between the airgap height h G and the gap diameter D is proposed

D k h G G .  (2)
where the coefficicent k G takes into account magnetic, mechanical, hydrodynamics and thermal considerations.

The relationship between B 1 , h G and magnets height h M is expressed as follow ) / . 1 /( .

1 M G r r h h B k B     (3)
) 2 / sin( ). / 4 (

 

  k where B r is the remanent flux density of the magnets,  is the magnet to pole width ratio and  r is the magnets relative permeability. This formula is of fair accuracy only in the case of a radial magnetic flux in the gap. An alternative and more exhaustive formula is given hereunder. This expression is derived from a 2D model proposed in [7] that solves the governing field equations by separating the polar variables. It predicts the open-circuit field distribution anywhere in the airgap of a slotless surface mounted PM machine and B 1 is expressed as a function of p

) 1 )( 1 ( ) )( 1 ( ) 1 /( 2 ). ) 1 ( 2 1 ( 2 2 2 2 2 2 1 1 1 p sm p rm r p rm p sm r p rm p rm p sm r R R R R p p R p R p R B k B                 ) 2 / /( 1 G M rm h D h R    ( 4
) ) / 2 1 /( 1 D h R G sm  
In addition, a coefficient k s , which takes into account the slotting effect, is applied to the airgap and magnet heights

) / /( . 1 r M G e o s h h R k      (5)
Two formulas are proposed for the reluctance R e in [START_REF] Matagne | Contribution à la modélisation des dispositifs électrotechniques en vue de leur modélisation[END_REF]. The first one shall be used in the case of a thin airgap

)))) 1 ( /( 1 1 .( /( ) . ( 1 t o r M G e k h h R          (6) ) ) 2 / ( 1 ln( ) / 2 ( ) 2 / ( tan ) / 2 ( 2 ' ' 1 M s s G M s h w w h h w       
where k t is the proportion of teeth, h' M is the magnetic height and w s is the slot width. A second formula shall be used in the case of a thick airgap

) 4 /( )) ln( ) 2 ln( ) 2 (( p mS D k k k k R pp o t t t t e      (7)
where S pp is the number of slots per pole and per phase and m the number of phases.

The slots and teeth height h S =h T depends on the rms electric load A L , as well as the rms current density J(A/m 2 ) in the slot conductors and the slott fill factor k

f )) 1 ( . /( t f L S k k J A h   (8)
The rotor and stator yoke minimum heights h Y (min) are chosen such that the flux density into the iron is lower than a maximum value B max (that generally corresponds to the saturation limit of the magnetic material). The following formula is determined by considering both superposed effects of magnets (height h YM ) and windings (height h YW ) on the iron flux density. The flux density in the airgap is assumed to be radial.

YW YM Y h h h   (min) (9) ) 4 /( . max 1 B pk B D h YM     ) 2 18 /( . . . 2 ' max 2 2 p h B D A h M o L YW   
It is important to note that the expression of h YW (effect of the windings) is given for the particular case of a three phase regular winding with one slot per pole and per phase.

In addition, the relationship between the gap diameter and the rotor inner diameter D int (m) is reminded

G M Y H h h h h D D 2 2 2 2 int      (10)
where h H (m) is an additional thickness that ensures a mechanical integrity to the rotor. The iron losses calculation is based on classical estimations of global losses p Fe per unit mass in each part of the stator magnetic circuit

c Feo Fe b o Fe Fe B B f f p p o ) / ( ) / .(  (11)
where f (Hz) and B Fe (T) are respectively the electrical frequency and flux density in the iron, p Feo (W/kg) is the iron losses per unit mass at a given frequency f o and flux density B Feo , b=1.5 and c=2.2, using typical medium quality Fe-Si laminated steel datasheets. For the calculation of the total losses P Fe , the flux density amplitude is supposed to be B max everywhere in the stator. The relationship between the EM torque T EM and the mechanical torque T Meca is

   / Fe Meca EM P T T (12)
In this relationship, we assume that the iron losses are mainly caused by the rotation of the rotor. The thruster's torque Q being an input data, if the mechanical losses are ignored for this study, then T Meca =Q. Now, let's summarize the sequential principle of this model that aims at determining the dimensions and performances of the electrical machine for a given set of specification. The propeller input data are the torque Q=T Meca , rotational speed  and diameter D P . For the special case of a rim propeller D P =D int . In addition, the following electrical parameters are fixed such that a unique solution to the equations can be found: A L , J, B 1 , p, k t ,  and k f . Assuming D≈D int , which tends to be a fair approximation for a rim-driven propeller (a thin machine with a large inner diameter), the airgap height is deduced from [START_REF] Abu-Sharkh | Structurally integrated brushless PM motor for miniature propeller thrusters[END_REF]. The slot height is deduced from ( 8) and the magnet height is roughly deduced from (3). The yokes height and the gap diameter are then deduced from ( 9) and [START_REF] Kuiper | The Wageningen propeller series[END_REF]. Finally, by ignoring the iron losses, the axial length is deduced from (1), the winding coefficient being set to 1 for this study (S pp =1) and the angle  being set to 0. Once all the dimensions are determined, it is then possible to calculate the iron losses and, thus, deduce the real EM torque from (12). Additionaly, the real flux density B 1R in the gap is deduced from (4). The real current loading A LR and density J R , are then deduced from (1) and [START_REF] Matagne | Contribution à la modélisation des dispositifs électrotechniques en vue de leur modélisation[END_REF]. The current in the conductors is then deduced from ) 2 /( . .

s LR mpn D A I   (13)
with n s the number of winding turns per phase and per pair of pole. In addition, the winding resistance (for each phase) is calculated.

cond S cond L w R /   ( 14
)
S cond is the conductor section which can be easily determined by the knowledge of k f , n s , and the slot dimensions. L cond is the total length of a winding conductor, it is composed of two elements: the axial resistance length and the end windings resistance length which can be estimated for each conductor in each end as a w pole diameter half-circle (for S pp =1).  (ohm.m) is the conductor resistivity. From equ. ( 13) and (14), it is then possible to calculate the copper losses If =0, the effect of the magnets and windings on the teeth saturation can be considered separately. Thus, the teeth saturation by the magnets is not reached as long as

) /( max 1 B k B k t   ( 17
)
In addition, the teeth saturation by the windings is not reached as long as

) 2 3 /( . . . . ' max M o L t h B D A k p    (18)
Again, this relationship is given for the particular case of a three phase regular winding with S pp =1. An additional constraint concerns the tooth shape that must follow the following criterion

max / R w h T S  ( 19
)
where w T is the teeth width and R max is a ratio that represents a limit in terms of teeth mechanical integrity.

Similarly, the magnets shape is chosen such that the ratio magnet height on magnet width remains realistic

2 max ) 2 / /( R p D h M   (20)
The constraint concerning the demagnetisation of the magnets is expressed as follow

cj G o G r L H h h B p D A   ' / ) / . ) 2 3 /( . . (   (21)
where H cj (A/m) is the coercive field of the magnets.

A thermal model has been built in order to limit the temperatures in the conductors to reasonable values. The constraint on the conductors maximum temperature T Cu (max) is simply expressed as follow max (max) T T Cu 

where T max is a limit temperature that depends on the conductor class. For a question of clarity, the thermal model is not detailled in this article but can be found in [START_REF] Drouen | Study of an innovative electrical machine fitted to marine current turbines[END_REF]. This model estimates roughly the temperatures in the different parts of the structure thanks to the dimensions and losses evaluated previously. It is based on a simple steady state thermal resistance network directly derived from the heat transfer equations under steady-state conditions. Additional mechanical constraints are considered. The first one is the thickness of the electrical machine h EM that must be lower than the duct thickness h duct that can be considered as directly proportional to the propeller diameter, i.e. h duct =k hduct .D P with k hduct <1

Fig. 2. Sketch of a cross section of 2 poles of the electrical machine

P hduct H Y S G M Y H D k h h h h h h h .        (23)
Secondly, the total length of the electrical machine L mach (end windings included) must be lower than the duct length L duct , considered as directly proportional to the propeller diameter, ie L duct =k Lduct .D P with k Lduct <1

P Lduct pp t D k m S p D k L . ) . . 2 /( ) 1 .( 2     (24)
Finally, an ultimate constraint concerns the voltage converter electrical frequency that must remain in a realistic range of values (max)

conv conv f f  (25)
The electrical frequency of the converter is directly dependant on the electrical frequency of the machine f mach , i.e. f conv > k fconv1 .f mach ) 2 /( .

 p f mach   ( 26
)
It also depends on the electrical machine time constant

 mach i.e. f conv > k fconv2 .1/ mach  ) /( ew a mach mach R R L    (27)
where L mach is the machine synchronous inductance (including the slot leakage inductance) that is calculated following classical equations [START_REF] Polinder | Comparison of direct-drive and geared generator concepts for wind turbines[END_REF], [START_REF] Grauers | Design of direct-driven permanent-magnet generators for wind turbines[END_REF].

The presented EM model has been validated in several typical sets of dimensions with numerical 2D and 3D simulations. This validation step is not presented in this paper for conciseness reasons.

III. PROPELLER MODEL

First of all, it seems important to remind the basic principle of a propeller. If we consider a section of blade (fig. 3) with a water inflow velocity V o (m/s) and a rotational speed =2n (rad.s -1

), then the relative velocity V R has an angle of attack  with the chord of the foil. It generates two forces: a lift force dL, normal to V R , and a drag force dF v , in the direction of V R . They both contribute to the thrust dT and torque dQ on the blade, that are the projections of the lift and drag forces on, respectively, axes x (the vessel trajectory) and z (the propeller plane). When working sufficiently far away from the free surface, the complete thrust T on a propeller of diameter D P can be expressed, on an exhaustive manner, as follow [START_REF] Carlton | Marine propellers & propulsion[END_REF] ) , , ( .

Q T o P K K J    (29)
For this study, we propose to use the Ka-N19A Wageningen series [START_REF] Kuiper | The Wageningen propeller series[END_REF] that give the performances of specific ducted propellers already tested at the Netherlands Ship Model Bassin. Torque and thrust coefficients K T and K Q , as well as efficiency  P can be expressed with polynomials in terms of advance coefficient J o and pitch ratio P/D for given blade area ratio and number of blades.

    2 ) 2 , 1 ( 2 1 ) 2 , 1 ( ) / .( . k k k k o k k T T D P J K  (30)
The interest of this model is that it is of good accuracy as directly derived from tests. Furthermore, thanks to this non dimensional approach, it is possible to deduce the performances of a propeller whatever its dimensions, which is particularly interesting in the case of a systematic design process. The only constraint is to keep a homothety on the whole propeller dimensions. Figure 4 represents the geometry of this type of propeller (with 5 blades).

The input parameters are the propeller thrust T, the water inflow velocity V o , both directly fixed by the vessel specification, the propeller diameter D P and the rotational speed =2n. T and V o are fixed parameters whereas D P and  can be considered as two free variables. Knowing T, V o , n and D P , it is then possible to calculate the thrust coeffcient K T and the advance coeffcient J o of the propeller. Figure 5 represents K T as a function of J o for different values of P/D varying between 0.5 and 1.4: obviously, one single pitch ratio P/D is fitted to a point (J o ,K T ). The model uses an iterative process to determine this value. Once the correct pitch ratio has been determined, it is then possible to deduce the torque Q and the efficiency  P of the propeller from the calculation of K T and K Q . The torque value is the output data that insures the main link between both propeller and electrical models.

Additionally, the mass of propeller is evaluated thanks to a reference masse M ref given for a Ka-N19 propeller of diameter D Pref . As all the dimensions follow a homothetic law, the masse M(D P ) for any diameter D P is 3 .

) / .( ) (

ref P P ref P D D M D M  (31)
IV. MODELS COUPLING Both models are coupled thanks to the following common input/output parameters: propeller diameter D P , torque Q and rotational speed =2n. In addition, the water speed V o can be considered as a common parameter as it is used in the thermal model (convective effects). The following scheme (fig. 6) describes, in a simplified form, the way both hydro. and EM/thermal models are coupled.

It is then interesting to study the influence of the propeller diameter and rotational speed on the overall perfomances of the system. An optimization of the system efficiency  elec × P , dimensions or mass are possible. Ideally, with this structure of model, it should be possible to add some more models that would represent other physical phenomenon that may have an influence on the performances of the machine. Phenomena such as mechanical distortion or viscous torque in the gap could be represented in future works.

V. EXAMPLE OF DESIGN OPTIMIZATION

For this study, a typical application is chosen and the set of specification of a patrol boat is considered. The water inflow velocity on the propeller is V o =19.87 knots. The thrust that must be delivered by the propeller is T=15.17 tons. It must be noted that this approach is simplified as, in the reality, the vessel speed V and hydrodynamic resistance R should be the real starting points of the study. Unfortunately, ratios V/V o and R/T are not constant and depend on the thruster geometry (thus vary with D P ). This could be taken into account thanks to a specific model of the vessel hull.

The diameter and rotational speed of the propeller vary on a realistic range of values, i.e. 1.0m ≤ D P ≤ 2.0m and 400 RPM ≤ 60n ≤ 600 RPM. Additional constraints such as maximum blade tip speed are considered in order to take into account the risks of cavitation max .

V D n p   (32)
for this specific study, V max is set to 50m/s which corresponds to a realistic value.

The proposed propeller has 4 blades and a fraction of surface A e /A o =0.75. It is important to note that the hydro. and EM models are not equivalent. Indeed, the hydro. model, based on a specific propeller geometry, has in reality few free variables: number of blades, fraction of surface but also chord and blade thickness are fixed. This is not the case of the EM model where the whole geometry of the electrical machine can be, theoretically, optimized.

As an exemple, the systematic design of the global efficiency of a rim-driven Ka-N19 propeller is presented. The aim is to maximize the global efficiency of the system. The permanent magnets are bonded NdFeB magnets: B r =0.6T, H cj =9.5.10 5 A/m and  r =1.20. Figure (7) gives the rim efficiency  RIM versus (D P ,60n) on a specific range of values. The optimum efficiency of the RIM  RIM =60.3%, is reached at (D P , 60n) = (1.62m, 517RPM). However, a non negligible range of solutions (D P , 60n) results in acceptable efficiencies ( RIM >60%). It should, potentially, permit a multi-objective optimization process (the minimization of the mass of the rotor for instance). It must be noted that unrealistic combinations are not displayed on the graph: either propellers of low diameter and speed that can't supply the requested power, or propellers of high diameter and speed that result in high blade tip speeds which leads to cavitation phenomenon.

For each point (D P , 60n), the following key variables are adjusted in order to optimize the electrical machine: A L , J, B 1 , p and k t . At the optimum point, the optimum values of In addition, the mass of the rotor (i.e. the propeller + the electrical rotor) of those machines is evaluated. The propeller mass is evaluated thanks to (31). Concerning the electrical rotor, the mass of the magnets and yoke is simply calculated thanks to the geometries evaluated by the EM model. The results are shown on fig. 8. It clearly reveals that the optimum in terms of efficiency is not ideal in terms of rotor mass (M rotor =2196kg). As it may be necessary to minimize the rotor mass to limit the constraints on the bearings of the thruster, a compromise between efficiency and mass seems necessary. The rotor mass is essentially dictated by the electrical rotor that tends to become lighter for higher speeds (lower torques for a given power). As an example, a possible alternative solution could be (D P , 60n) = (1.52m, 578RPM) where M rotor = 1939kg for a global efficiency =60.1%. The mass is reduced by 11.7% for a loss of efficiency of 0.3%, which corresponds to an electrical power increase of about 5.7kW (the mechanical power delivered to the propeller is T.V o =1521kW).

VI. CONCLUSION

In this paper, a multi-physic model for the design of an integrated rim driven propeller is presented. The rim-driven thruster is made of two main sub-systems (the propeller and the electrical machine) that are closely related. As a consequence, a coupled multi-physic model for the design of this machine seems essential. This paper gives a description of two hydrodynamics and electromagnetic / thermal models as as the way to couple them. The simplicity and the accuracy of the proposed analytical coupled model allow an easy insertion in a systematic design process. An example of systematic design is proposed where the efficiency and the rotor mass of the system are optimized. It is shown, on this particular example, how to select appropriately the characteristics of the propeller. This example highlights the interest of a coupled model on the design of a rim propulsion system. Some future improvements on the model should include the development of a more exhaustive propeller model, adapted to a more accurate optimisation process.

Fig. 1

 1 Fig. 1 Schematic view of a rim-driven system

  as saturation, demagnetisation and manufacturing constraints are also considered. If those constraints are correctly defined, it results in a very robust tool that eliminates any unrealistic solution.

  the advance coefficient, R n and  o are the Reynolds and cavitation non dimensional numbers. It is a common design practice, for given Reynolds and cavitation conditions, to express the non dimensional number K T coefficient, as a function of the advance coefficient J o . In the same way, the torque coefficient K Q =Q/ W n 2 D 5 is expressed as a function of J o for given Reynols and cavitation conditions. Those two non dimensional numbers characterize the performances of the propeller in terms of torque, thrust, as well as efficiency  P

Fig. 3 A

 3 Fig. 3 A blade profile

Fig. 5

 5 Fig. 5 KT (continuous lines) and  (dashed lines) function of Jo for different pitch ratios (P/D varies between 0.5 and 1.4)

Fig. 6

 6 Fig. 6 Coupled hydro and EM/thermal models


T Henneron

F Locment

E Semail

F Piriou

Predetermination of Currents and Field in Short-Circuit Voltage Operation for an Axial-Flux Permanent Magnet Machine

Keywords: 3D Finite Element Method, Short-circuit fault, seven-phase machine, Axial Flux Permanent Magnet Machine, Soft Magnetic Composite

Risk of irreversible magnet demagnetization during short-circuit fault is analyzed in case of an axial-flux dual-rotor machine, using a three-dimensional finite-element method (3D-FEM). In order to validate the numerical model, calculated waveforms of the currents are compared with experimental results for short-circuit at low speeds. Then currents and magnetic flux density inside the magnets are computed for short-circuit at higher speeds in order to predetermine the maximum admissible speed for the machine.

I. INTRODUCTION

Fault tolerance of electrical machine is an important issue for embedded systems as offshore wind generators. Permanent Magnet machines with more than three phases are attractive [START_REF] Levi | Multiphase induction motor drives: a technology status review[END_REF] [2]. Since they can still work with several open-circuited phases. However, it is also necessary to confirm their ability to withstand short-circuit fault currents without irreversible magnet demagnetization [START_REF] Gyu-Hong | Analysis of irreversible magnet demagnetization in line-start motors based on the finite-element method[END_REF]. Besides, new material and geometrical constraints yields interest forward Axial Flux Permanent Magnet synchronous machines (AFPM) [START_REF] Gieras | Axial Flux Permanent Magnet[END_REF] [START_REF] Chen | Axial-flux PM wind generator with a soft magnetic composite core[END_REF].

The paper tackles with the predetermination of shortcircuit currents in a seven-phase AFPM generator with two external rotors and iron-powder stator core. Even if thermal effect of short-circuit current is acceptable, it is necessary to verify that, during the time between the occurrence of a shortcircuit fault and the action of the overload protections, irreversible demagnetization of the magnet is not reached. Thus, determination of magnetic flux density inside the magnets is of great interest.

In section II, the small scale AFPM prototype, whose three-dimensional characteristics yields a 3D-FEM modeling, is described. In section III, the 3D numerical model that takes into account coupling with external electrical circuit is presented. In section IV, numerical model is validated using experimental results for different cases: without short-circuit; with a short-circuit fault at low speed. Finally, the numerical model is used for determination of short-circuit currents and magnetic flux density in magnets at higher speeds.

II. PRESENTATION OF THE MACHINE

The AFPM generator is composed of a stator and two externals six-pole rotors with rare-earth permanent magnets (Fig. 1). Since the magnetic flux paths are in all three directions of the magnetic circuit, a Soft Magnetic Material (SMC) with isotropic magnetic properties has been used for the stator core. The seven phases are obtained with toroidal coils distributed into 42 slots. The two rotors are identical but with an angular shift of 360/84 degrees between them in order to reduce the cogging torque. The spatial magnet repartition implies non-sinusoidal electromotive forces as it is usual for multiphase machines. All these specificities imply that analysis by FEM is necessary. Rated torque value of 65Nm torque is obtained for a 5.1A RMS current, rated power of 5.1kW rated power, for a 750RPM speed. More precise data are given in [START_REF] Locment | Soft Magnetic Composite Axial Flux Seven-Phase Machine[END_REF].

B. Study of flux magnetic density inside a permanent magnet

By using a numerical model, it is possible to determine the drop of the magnetic flux density inside the magnet due to the demagnetizing short-circuit currents. In Fig. 9, the magnetic flux density at the center of a permanent magnet, with and without short-circuit fault, is given at 30rad/s. In both cases, the effect of the slots (7 slots for 60°) on the variations of the flux density is obvious. In case B, it appears also the demagnetizing effect of the short-circuit current in two adjacent phases. The two coupled phenomena imply a complex waveform for the magnetic flux density. Nevertheless, it is easy with the obtained results to check that the minimum value that yields irreversible demagnetization is not reached. 10 presents the variation of the lowest value of the magnetic flux density inside a permanent magnet versus the rotation speed for both studies. In the case of the study B, the magnitude of the short circuit current on the sixth and seventh winding is more important compared to the others windings. The effect of short-circuit currents for an AFPM with two shifted rotors and open slots has been calculated with a 3D-FEM. Comparison of calculated currents with experimental ones in case of short-circuit at low speed has validated the model that has been then used to predetermine the complex magnetic flux density in the center of a magnet. It is thus possible to check that irreversible demagnetization is not reached before the action of the overload protection devices. The numerical model can be still improved by modeling the SMC with non-linear characteristics that becomes necessary at high speeds with high short-circuit currents.

Fig. 1 .Fig. 2 .Fig. 7 .Fig. 8 .

 1278 Fig. 1. Studied machineFor the numerical model, the small shift between the rotors imposes the mesh (791506 tetrahedron elements and 141937 nodes) of one sixth of the machine represented in Fig.2.

Fig. 9 .

 9 Fig. 9. Magnetic flux density inside a permanent magnet

Fig.

  Fig.10presents the variation of the lowest value of the magnetic flux density inside a permanent magnet versus the rotation speed for both studies. In the case of the study B, the magnitude of the short circuit current on the sixth and seventh winding is more important compared to the others windings.

Fig. 10 .

 10 Fig. 10. Lowest value of the magnetic flux density inside a permanent magnet V. CONCLUSION


P Dupont

A C Bayeul-Lainé

A Dazin

G Bois

O Roussette

Q Si

Leakage flow simulation in a specific pump model

This paper deals with the influence of leakage flow existing in SHF pump model on the analysis of internal flow behaviour inside the vane diffuser of the pump model performance using both experiments and calculations. PIV measurements have been performed at different hub to shroud planes inside one diffuser channel passage for a given speed of rotation and various flow rates. For each operating condition, the PIV measurements have been trigged with different angular impeller positions. The performances and the static pressure rise of the diffuser were also measured using a three-hole probe. The numerical simulations were carried out with Star CCM+ 8.06 code (RANS frozen and unsteady calculations). Comparisons between numerical and experimental results are presented and discussed for three flow rates. The performances of the diffuser obtained by numerical simulation results are compared to the performances obtained by three-hole probe indications. The comparisons show few influence of fluid leakage on global performances but a real improvement concerning the efficiency of the impeller, the pump and the velocity distributions. These results show that leakage is an important parameter that has to be taken into account in order to make improved comparisons between numerical approaches and experiments in such a specific model set up.

Introduction

This work is an extended study that had originated in LML laboratory concerning the well-known SHF (Société Hydrotechnique de France) pump. A lot of experimental data have already been produced (tests in air and in water) on that geometry and these results are used as databases for the validation of numerical calculation results and for new experimental results [START_REF] Wuibaut | Experimental analysis of velocities in the outlet part of a radial flow pump impeller and the vanelessdiffuser using particle image velocimetryProceedings[END_REF][START_REF] Wuibaut | Etude par vélocimétrie par images de particules des interactions rouediffuseurdansunepompe centrifuge[END_REF][START_REF] Cavazzini | Experimental and numerical investigation of the rotor-stator interaction in radial turbomachines[END_REF][START_REF] Cavazzini | Analysis of the rotor-stator interaction in a radial flow pump Houille blanche[END_REF][START_REF] Cavazzini | Analysis of unsteady flow velocity fields inside the impeller of a radial flow pump : PIV measurements and numerical calculation comparisons[END_REF][START_REF] Cherdieu | Data reduction problems using a 3 -hole directional pressure probe to investigate mean flow characteristics in the vaneless gap between impeller and diffuser radial pump[END_REF][START_REF] Bayeul-Lainé | Numerical and experimental investigations in a vaneddiffuser of SHF impeller : fluid leakage effect 21 èmeCongrèsFrançais de Mécanique[END_REF][START_REF] Bayeul-Lainé | Fluid leakage effect on analysis of a vaneddiffuser of SHF pump[END_REF] Flow behaviour in a radial machine is quite complex and strongly depends on rotor stator interactions and operating conditions [START_REF] Adamczyk | Wake induced unsteady flows :their impact on rotor performance and wake rectification[END_REF][START_REF] Arndt | -Experimental Investigation of Rotor -Stator Interaction in a Centrifugal Pump With Several VanedDiffusers[END_REF][START_REF] Eisele | Flow analysis in a Pump Diffuser Part 1 : LDA and PTV Measurements of the Unsteady FlowTransactions of ASME[END_REF].

In numerical simulations, two aspects have to be considered: -The first one concerns the governing equations which are solved in the model. Three kinds of numerical calculations are currently used in turbo machinery (frozen rotor calculations, mixing plane calculations, unsteady calculation). The frozen rotor calculation is the steady state method which uses the rotating reference frame to save the computational resources by converting inherently transient turbo-machinery flow into steady state. The difference between a frozen rotor and a mixing plane calculation is that the mixing plane mixes the flow and applies the average quantities on the interface for upstream and downstream components, while a frozen rotor will pass the true flow to downstream and vice versa. So if the wake effect on the downstream component performance is necessary, the frozen rotor method has to be used. A big disadvantage is that it gives the solution at a single relative position. The true transient method gives the wake effect on the downstream component for all relative positions (as it happens in reality) [START_REF] Culver | Plane Method for Flutter Computation in Multi-stage Turbomachines47[END_REF][START_REF] Petit | Numerical Investigations of Unsteady Flow in a Centrifugal Pump with a vanedDiffuserInternational[END_REF][START_REF] Tamm | Experimental and 3-D Numerical Analysis of the Flow Field in Turbomachines-Part[END_REF].

-The second aspect concerns the geometrical model. Some geometrical simplifications are currently used. For example, the leakage flows are often neglected. It is obvious that a complex model (fully unsteady, with leakage flows) will be more time consuming but closer to the real physics. In the presence of a gap, a leakage flow appears due to the difference of static pressure at the clearance between the shroud and the outer casing and between the hub and the outer casing.

In this respect, the effects of fluid leakage due to the gap between the rotating and the fixed part of this pump model are analysed and discussed.

In this paper, some limits of the numerical models are pointed out. To do so, several numerical calculations have been carried out: i-Frozen rotor without leakage, ii-Unsteady RANS calculations without leakage, iii-Frozen rotor with leakage and iv-Unsteady RANS calculations with leakage These numerical results are compared to experimental results for three flow rates to try to evaluate the effects of the numerical models on the prediction of the performance and on the local behaviour of the flow. 2) was the same as the one already developed for the PIV measurement already described in the previous papers, especially in reference [START_REF] Cavazzini | Analysis of the rotor-stator interaction in a radial flow pump Houille blanche[END_REF]. Test pump model and PIV measurements conditions have been already described in several papers [START_REF] Cavazzini | Experimental and numerical investigation of the rotor-stator interaction in radial turbomachines[END_REF][START_REF] Cavazzini | Analysis of the rotor-stator interaction in a radial flow pump Houille blanche[END_REF] and main pump characteristics are given in table 1.

Experiments

Test and apparatus

Impeller

Two sources of fluid leakage occur:

-The first one at the impeller inlet (leakage 1, figure 6) -The second one between the impeller outlet and the diffuser both for the hub and the shroud sides (leakages 2 and 3, figure 6).

PIV measurements

PIV measurements have been performed at different hub to shroud planes inside one diffuser channel passage for a given speed of rotation and various flow rates. For each operating condition, the PIV measurements have been trigged with different angular impeller positions (figure 3). For each angular position, several instantaneous velocities charts have been obtained, covering the space between inlet and outlet diffuser throats. This makes a rather good evaluation of phase averaged velocity charts possible. The PIV results are extracted from the thesis of G. Cavazzini, [START_REF] Wuibaut | Etude par vélocimétrie par images de particules des interactions rouediffuseurdansunepompe centrifuge[END_REF].

Three holes probes

A directional three holes probe has been used to make hub to shroud traverses [START_REF] Cherdieu | Contrôle du décollementdans un diffuseuraubé de turbomachine centrifuge[END_REF]. Using a specific calibration one can get total pressure, static pressure, absolute velocity and its two components in radial and tangential direction.

In order to well represent the flow field, twenty-three probe locations are defined as it can be seen in figure 5. For each location, ten axial positions are registered (b*=0.125, 0.2, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 0.925, 0.975 from hub to shroud). The present analysis focuses only on locations 19 to 23 along the blade to blade diffuser channel.

Calculations

Calculations were carried out with Star CCM+ 8.06 code. The sketches of the gaps are really taken account (figure 6). A complete meshing was set up with the real geometry of the gaps and external domain far upstream the labyrinth. In this respect, boundary conditions include the relative casing motion of the model.

The calculation domain is divided into four domains:

-Inlet domain, -Impeller domain, -Diffuser domain, -Exterior domain as can be seen in Figure 6.

As thickness gaps are very small (1 mm), a thin meshing model has been chosen. A thin meshing model allows thin regions in the geometry to have a prismatic type volume mesh. Using this kind of mesh improves the overall cell quality and reduces the cell amount when compared to an equivalent polyhedral type core mesh. The number of layers in the thin mesh has been set to 10. In the other parts of the domain, a polyhedral with prism layers is used for all calculations (5 prism layers for a total prism layer thickness of 1 mm). The target size is 3 mm and the minimum size 0,5 mm . The size of the grid is about 18. 10 6 cells for all the four domains. For frozen rotor calculations, line probes are plotted in the whole domain of the diffuser (figure 4): each probe is duplicated seven times in order to obtain all parameters (pressure, total pressure, radial, tangential and axial velocity and velocity magnitude) for different relative positions of the diffuser comparatively to the impeller (figure 3). This is equivalent to azimuthal positions equal to n* 360/(Zi*Zd) with n=0,1 ,2, 3, 4, 5, 6 and 7.

For unsteady calculations, point probes are plotted in the blade to blade channel of the diffuser as can be seen in figure 4.

The convergence criteria are less than 1.e-4. The values of y+ are below 15 in the whole computational domain.

Results and analysis

Let us first examine the global results: the non-dimensional static pump head, the non-dimensional isentropic pump head, the impeller efficiency, the diffuser efficiency and the pump efficiency.

These calculations show few influence of fluid leakage on global performances (figure 7). Global results of static theoretical head pump are in good agreement between one dimensional approximation using EULER equation with no inlet flow swirl and frozen rotor calculation. Unsteady calculation results, performed for the three non-dimensional flow rates, give the same level, as it can be seen in figure 7. The numerical non-dimensional isentropic pump head is higher when the fluid leakages are modelled in case of frozen rotor and in case of unsteady calculations. This leads to a decrease of impeller efficiency and pump efficiency due to fluid leakages (figure 8). It can also be noticed that the pump efficiency falls at high flow rate. The diffuser efficiency is higher at low flow rate and decrease with flow rate with and without fluid leakages for frozen rotor calculations. This is due to the fact that the diffuser is not well adapted at high flow rate. It is well adapted at Q*=0.766 as can be observed in figure 8 Fluid leakages are calculated for the three different locations: inlet (leakage 1), hub side (leakage 2) and shroud side (leakage 3) with frozen rotor and unsteady calculations. The impeller rotation effect was taken into account in the numerical solution.

Wuibaut [START_REF] Wuibaut | Etude par vélocimétrie par images de particules des interactions rouediffuseurdansunepompe centrifuge[END_REF] and Cherdieu [START_REF] Cherdieu | Contrôle du décollementdans un diffuseuraubé de turbomachine centrifuge[END_REF] calculated analytical solutions, based on friction losses and localized minor losses of charge, for each leakage. The impeller rotation effect was not taken into account in this analytical solution. By writing the Bernoulli's theorem, each leakage flow rate has been calculated knowing the total difference pressure between the inlet and the outlet of the gap.

Considering the leakage at the impeller inlet, only localized losses of charge in the labyrinth are taken account [START_REF] Wuibaut | Etude par vélocimétrie par images de particules des interactions rouediffuseurdansunepompe centrifuge[END_REF]. For the gap between the impeller and the diffuser on the shroud side, the losses are due to the contraction and the expansion through the gap for the localized part and also due to friction on the length of the ring [START_REF] Cherdieu | Contrôle du décollementdans un diffuseuraubé de turbomachine centrifuge[END_REF]. For the gap between the impeller and the diffuser on the hub side, Cherdieu considered that the losses are the sum of two elbows, the contraction and the expansion [START_REF] Cherdieu | Contrôle du décollementdans un diffuseuraubé de turbomachine centrifuge[END_REF]. As can be observed, at the inlet gap, all methods lead to quite same results. In case of Star CCM+, fluid leakages are smaller at high flow rate.

Concerning, leakages at hub and shroud sides, the amount of flow rate depends on the flow rate in the pump. At low values of Q*, a positive leakage between impeller and diffuser was observed. At high values the tendency inverses as shown in figures 10 and 11.These tendencies were also observed in previous studies [START_REF] Bayeul-Lainé | Numerical and experimental investigations in a vaneddiffuser of SHF impeller : fluid leakage effect 21 èmeCongrèsFrançais de Mécanique[END_REF].

Considering fluid leakages flow rates, the analytical results seem to be sufficient enough to predict these flow rates as also shown by Khelladi et al [START_REF] Khelladi | Study of a high rotational speed shrouded centrifugal, fan: aerodynamics and effects of a shroud-associatedProceedings of the Institution of Mechanical Engineers Part A[END_REF]. So new calculations, without external domain but with fluid leakages flow rates based on analytical model, can be used to compare local numerical results and experimental results.

Local absolute velocity angle relative to radial direction and non-dimensional radial velocity are plotted in figures 12 to 17 for the three flow rates.

The examination of the evolution of the absolute velocity angle relative to radial direction for the two numerical calculations, for the PIV measurements and for the three holes measurements show similar results (figures 12 to 14). There is no real influence of taken account the fluid leakages on the evolution of the angles.

The examination of radial velocity (figures 15 to 17) shows that fluid leakages have a real influence for Q*=0.776 and Q*=0.973. For Q*=1.173, numerical results with fluid leakages are better than those without leakage, but some differences still remain with experimental results in particular flow at low values of b*.

Conclusion

In this paper, the SHF pump was studied numerically and experimentally. The computational models with and without gaps (between the impeller, the fixed inlet domain and the fixed diffuser domain) were set up for frozen rotor and unsteady solutions in order to determine the leakages and to make comparisons with experimental results. A one dimensional loss model was also used to calculate the fluid leakages flow rates. The effects of the gaps on the calculated flow field were analysed. Local velocity components evolutions in the diffuser blade to blade channel have been calculated and analysed taking the leakage effects into account.

It has been found that the amount of fluid leakages corresponds to the one dimensional model, even if the effect of rotating wall wasn't considered in this model. This result is very interesting knowing that the model with gaps is more consuming in grid size and in computational time than the case without gaps.

It has been shown that leakages may affect the local flow characteristics depending on volume flow rates at the impeller inlet. The calculations show little influence of fluid leakage on the global performance but present a real improvement concerning velocity distributions. The local fields of velocities were analysed inside the blade to blade channel using the numerical results, the PIV measurements and the three holes probes. These results show that leakage is an important parameter that has to be taken into account in order to make improved comparisons between numerical approaches and experiments.

The comparison of pump performances between the calculations and the measurements showed a favourable tendency.

Nomenclature

Figure 1 .

 1 Figure 1. Apparatus: PIV measurement, laser sheet.Figure 2. Apparatus: three-hole probe.

Figure 2 .

 2 Figure 1. Apparatus: PIV measurement, laser sheet.Figure 2. Apparatus: three-hole probe.

Figure 3 .

 3 Figure 3. Impeller different angular positions relative to the diffuser vanes.

Figure 4 .

 4 Figure 4. Probes for frozen rotor calculations.

Figure 5 .

 5 Figure 5. Diffuser measurement locations for probe traverse and unsteady calculations.

Figure 6 .

 6 Figure 6. Regions modelled. Details of fluid leakages.

Figure 7 .

 7 Figure 7. Performances: static and isentropic .Figure 8. Impeller and pump efficiency 

Figure 8 .Figure 9 .Figure 10 .Figure 11 .

 891011 Figure 7. Performances: static and isentropic .Figure 8. Impeller and pump efficiency 

Figure 12 .Figure 15 .

 1215 Figure 12. for Q*=0.776 Figure 13. for Q*=0.973 Figure 14. for Q*=1.134

Table 1 .

 1 Pump characteristics.

	Diffuser