Julien Fontchastagner
email: julien.fontchastagner@n7.fr

Frédéric Messine
email: frederic.messine@n7.fr

Yvan Lefèvre
email: yvan.lefevre@n7.fr

A New Rational Way Combining Analytical and Numerical Models with a Deterministic Global Optimization Algorithm for the Design of Electrical Rotating Machines

Keywords: analytical model, numerical model, inverse problem of design, finite elements methods, deterministic global optimization, interval Branch and Bound algorithm

The purpose of this paper is to present a new method of design which is an extension of previous works combining analytical models with exact global optimization algorithms. The idea here is to associate analytical and numerical models of electrical machines with an extension of the exact global optimization algorithm previously used in order to solve more general design problems. Some examples validate this approach by comparing the optimal solutions found with the analytical models and those combining analytical and numerical models.

I. INTRODUCTION

N OWADAYS, the problem of the design of electrical machines is understood and formulated as an inverse problem. The direct problem of design can be defined as follows: From an electromagnetical actuator where the structure, the dimensions and the composition are known, compute some characteristic values; for example, the flux density, the torque, etc. The corresponding inverse problem of design is: From the characteristic values given by the schedule of conditions (for example the torque), get the structure, the dimensions and the composition of the required actuator. Such a problem is ill-posed in the Hadamart sense. Indeed, the existence and the uniqueness of the solution cannot be guaranteed. Furthermore, this problem may generate an infinity of solutions. Many counter-examples can be found, like the machine considered in [START_REF] Messine | Optimal Design of Electromecanical Actuators: A New Method Based on Global Optimisation[END_REF].

In [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF], the resolution of such inverse problems are considered by formulating the inverse problems of design as global optimization problems. Therefore in [START_REF] Fitan | A general analytical model of electrical permanent magnet machine dedicated to optimal design[END_REF], [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF], a rational way associating analytical combinatorial models of electrical rotating machines and an exact global optimization algorithm, named IBBA, was proposed and studied. IBBA is an exact global optimization method based on a Branch and Bound technique which uses interval analysis tools, see [START_REF] Hansen | Global Optimization using Interval Analysis[END_REF], [START_REF] Keafott | Rigorous Global Search: Continuous Problems[END_REF], [START_REF] Ratschek | New Computer Methods for Global Optimization[END_REF] for details on such algorithms and on its rigorous convergence to the global optimum. We propose in this work to find the exact solution of an optimal design problem which directly satisfies using numerical tools the imposed schedule of conditions. Indeed, when the so-obtained optimal solutions are validated by the means of numerical tools (such as finite element methods), we denote some differences about the values of the electromagnetic torque. Thus, this can involve some adjustments of the parameters of the so-obtained machines. These adjustments can be done by solving iteratively the direct problem of design until the schedule of conditions will be satisfied (using a finite element method).

In order to make easier, the resolution of the direct problem of design by using a finite element method, we have proposed a numerical tool, named NUMT, which can automatically mesh, draw and compute the torque of an electrical machine only defined by its parameters of design, [START_REF] Lefèvre | Building a CAD system for educational purpose based only on a mesh tool and a FE solver[END_REF]; the computation of the electromagnetical torque can be performed without the drawing of the motor.

The purpose of this paper is to extend the algorithm IBBA by inserting some steps of NUMT in order to directly solve a more interesting inverse problem of design where the solutions do not need the validation phase. In Section II, we recall the rational methodology published in [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF]. We then present the numerical tool NUMT in Section III, [START_REF] Lefèvre | Building a CAD system for educational purpose based only on a mesh tool and a FE solver[END_REF]. Section IV is dedicated to the new methodology combining IBBA and NUMT. In Section V, the new algorithm named IBBA+NUMT is validated on some examples of design and the solutions are discussed and compared with those produced by IBBA alone. [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF] The purpose of the paper [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF] was to propose a rational methodology for solving the inverse problem of design. Thus, new analytical models, named combinatorial analytical models, allowed us to take into account a lot of possible rotating electrical machine with permanent magnets. These combinatorial analytical models were done by introducing discrete variables into dimensional analytical model; for example the number of pole pairs, the kind of structure (internal or external rotor configuration) and the kind of materials used for the magnets. We then obtain a combinatorial analytical model which represents a large part of electrical rotating machines with magnetic effects, [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF]. Combining this general combinatorial model with an efficient exact global optimization code, named IBBA for Interval Branch and Bound Algorithm [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF], [START_REF] Messine | A Deterministic Global Optimization Algorithm for Design Problems, chapter in Essays and Surveys in Global Optimization[END_REF], some exact optimal solutions were found by minimizing the magnet volume, the active part volume, the total volume, the mass or a combination of them with an imposed torque, see [START_REF] Fitan | A general analytical model of electrical permanent magnet machine dedicated to optimal design[END_REF], [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF] for some examples. This methodology is perfect in the first steps of design of an electrical machine answering to an imposed schedule of conditions.

II. THE RATIONAL METHODOLOGY OF DESIGN OF

The inverse problems, that we consider, are formulated in [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF] as mixed constrained global optimization problems:

           min x∈I R n ,z∈I N m , σ∈ Q l i=1 K i ,b∈B r f (x, z, σ, b) subjected to g i (x, z, σ, b) ≤ 0 ∀i ∈ {1, ..., p} h j (x, z, σ, b) = 0 ∀j ∈ {1, ..., q} (1)
where f is a real function, K i represents an enumerated set of categorical variables, for example the type of magnet, and B = {0, 1} the boolean set which is used to modelize the fact that an actuator is with or without slot(s) for example. IR and IN are respectively the real and the positive integer sets. This formulation is called optimal design optimization problem and answers perfectly to the inverse problem of the design of electromechanical actuators, see [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF], [START_REF] Messine | Optimal Design of Electromecanical Actuators: A New Method Based on Global Optimisation[END_REF] for more details about this formulation.

To solve these problems (1), we must use an exact global optimization algorithm in order to characterize the solution of the problem which can established that one structure is most efficient than another (with respect of the dimensions). For the use of IBBA, all the functions must be explicitly defined.

A. IBBA Algotihtm

Interval analysis was introduced by Moore [START_REF] Moore | Interval Analysis[END_REF] in order to control the propagation of numerical errors due to floating point computations. Thus, Moore proposes to enclose all real values by an interval where the bounds are the two closest floating point numbers. Then expanding the classical operations -addition, subtraction, multiplication and division -into intervals, defines interval arithmetic. A straightforward generalization allows computation of reliable bounds (excluding the problem of numerical errors) of a function over a hypercube (or box) defined by an interval vector. Moreover, classical tools of analysis such as Taylor expansions can be used together with interval arithmetic to compute more precise bounds [START_REF] Moore | Interval Analysis[END_REF]. Other new bounding techniques include combining linear bounds at all vertices of the box [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF] or use affine arithmetic, [START_REF] Messine | Extension of Affine Arithmetic: Application to Unconstraineded Global Optimisation[END_REF], [START_REF] Messine | Exact and Rigorous Global Optimization Algorithm based on Affine Arithmetic and its Extensions to Quadratic Forms[END_REF]. Extensions of these methods are proposed in [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF], [START_REF] Messine | A Deterministic Global Optimization Algorithm for Design Problems, chapter in Essays and Surveys in Global Optimization[END_REF] in order to solve mixed (discrete and continuous) problems of type [START_REF] Bastos | Electromagnetic Modeling By Finite Element Methods[END_REF].

The principle of IBBA is to bisect the initial domain where the solution is sought for into smaller and smaller boxes, and then to eliminate the boxes where the global optimum cannot occur:

• by proving, using interval bounds, that no point in a box can produce a better solution than the current best one; • by proving (with interval arithmetic), that at least one constraint cannot be satisfied by any point in such a box.

To accelerate the convergence, constraint propagation techniques are used in some steps of IBBA, see [START_REF] Messine | Deterministic Global Optimization using Interval Constraint Propagation Techniques[END_REF] for details. The principle is to use, a priori, the implicit relations between the variables which are induced by the constraints in order to reduce the size of a box. Such interval Branch and Bound algorithms guarantee to produce an -global optimal solution, where (> 0) is the maximal error on the objective function value. For details and rigorous convergence analysis of these deterministic global optimization methods based on interval analysis, the reader is invited to consult to the three following books [START_REF] Hansen | Global Optimization using Interval Analysis[END_REF], [START_REF] Keafott | Rigorous Global Search: Continuous Problems[END_REF], [START_REF] Ratschek | New Computer Methods for Global Optimization[END_REF]. For details on IBBA dedicated to solve electromagnetical rotating machines and other actuators, see [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF], [START_REF] Messine | A Deterministic Global Optimization Algorithm for Design Problems, chapter in Essays and Surveys in Global Optimization[END_REF], [START_REF] Messine | An Interval Branch and Bound Method Dedicated to the Optimal Design of Piezoelectric Actuators[END_REF], [START_REF] Messine | Optimal Design of Electromecanical Actuators: A New Method Based on Global Optimisation[END_REF].

B. Combinatorial Models for Electrical Machines

Hereafter, the analytical equations of the magnetical model are recalled for the electromagnetical torque, see [START_REF] Fitan | A general analytical model of electrical permanent magnet machine dedicated to optimal design[END_REF], [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF] for details.

All the parameters of a rotating electrical machines are: -D(m) represents the bore diameter, -L(m) is the length,l a (m) the thickness of the permanent magnets, -E(m) the winding thickness, -C(m) the thickness of yoke, -β the polar arc factor, -g(m) the thickness of the mechanical airgap, -p the number of pole pairs, J(σ m ) the magnetic polarization which depends on the categorical variable σ m representing the type of permanent magnet, another categorical variable denoted by σ mt defines the type of magnetic conductor, etc, for details see [START_REF] Fitan | A general analytical model of electrical permanent magnet machine dedicated to optimal design[END_REF], [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF].

Γ em (D, L, ...) = k Γ D [D + (1 -b e )(2b r -1)E] LB e K S , K S (D, L, ...) = k r Ej b e a a + d + (1 -b e ) , k Γ (D, L, ...) = π 2 b f [1 -K f ] β +(1 -b f ) √ 2 2 sin(β π 2 ) , K f (D, L, ...) = 1.5pβ E + g D (1 -b e ).b f , B e (D, L, ...) = 2J(σ m )l a (2b r -1)D ln D+2E(2br-1)(1-be) D-2(2b r -1)[l a +g] 1 k c , k c (D, L, ...) = 1 1 -b e N e a 2 5πD.g+πD.a , . . .
where the generic expression of the electromagnetic torque is denoted by Γ em . K S represents the current electric loading. According to the considered kind of armature (non-slotted or slotted), K S is identified with two distinct functions. In the case of non-slotted machines, this function is written K S = k r Ej whereas for slotted machines it is written K S = k r Eja/(d + a); where j is the current density. A generic formulation of the current electric loading can then be elaborated by introducing a boolean variable (zero or one) b e .

When b e has value zero, non-slotted machines are considered, and when it takes value one, slotted machines are taken into account. k Γ is the torque coefficient, the expression of which depends mainly on the kind of waveform which has been chosen (sinusoidal or rectangular). This coefficient is written π 2 (1 -K f ) √ β for rectangular waveform machines and π √ 2 4 sin(β π 2 ) for sinusoidal waveform machines. The elaboration of a generic expression is proposed by introducing a new boolean variable b f . If b f is equal to zero, sinusoidal waveform machines are considered or if b f is equal to one rectangular waveform machines are rather taken into account. Concerning non-slotted machines, a semi-empiric magnetic leakage K f is proposed. B e represents the no-load magnetic radial flux density to the bore diameter neighborhood, which is supposed purely radial in the airgap. An analytical expression for non-slotted machines with internal rotor has already been elaborated [START_REF] Messine | Optimal Design of Electromecanical Actuators: A New Method Based on Global Optimisation[END_REF]. A generic formulation of B e allows on the one hand to take into account the kind of armature by using the b e boolean variable, and on the other hand to take into account the rotoric configuration (internal or external) by introducing an extra boolean variable b r . For details on such a model see [START_REF] Fitan | A general analytical model of electrical permanent magnet machine dedicated to optimal design[END_REF], [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF].

The other relations define the volumes of the active parts of the machines, the permanent magnet volume V m , the yoke volume V c , the teeth or wedge volume V d , the electrical conductors volume V co and the global volume V g , by considering that the slots, wedges and magnets have a radial geometrical form.

V m (D, L, b e , ...) = βπLl a [D -S(σ r )[2g + l a ]] , (2) V c (D, L, b e , ...) = 2πLC [D + S(b r ) [E -g -l a ]] , (3) V d (D, L, b e , ...) = πLE [D + S(b r )E] [[1 -β]b e +[ d d + a ](1 -b e ) , ( 4
)
V co (D, L, b e , ...) = k r πLE [D + S(b r )E] [β.(1 -b e ) +[ a d + a ]b e , (5)
V g (D, L, b e , ...) = πL 4 b r [D + 2 [E + C]] 2 +(1 -b r ) [D + 2 [g + l a + C]] 2 . ( 6
)
A last relation is dedicated to the mass of the active parts M a :

M a (D, L, b e , ...) = V m .ρ PM (σ m ) + V c .ρ CM (σ mt ) +V co .ρ co + V d [ρ CM (σ mt ).b e + ρ Al .(1 -b e )] (7)
where ρ Al , ρ co , ρ CM (σ mt ) and ρ P M (σ m ) are the densities of respectively the aluminum, the copper, the magnetic conductor σ mt and the permanent magnet σ m . Relations linked to volumes and weight will be the objective functions f of the corresponding global optimization problems.

The strong equality constraint is about the torque Γ em which is fixed to a value by the schedule of conditions, denoted in this paper Γ; i.e. Γ em (D, L, b e , ...) = Γ.

III. A NUMERICAL TOOL TO VALIDATE OPTIMAL SOLUTIONS Before the phase of prototype making the optimal solutions obtained by the methodology described in [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF] need to be validated by using numerical tools, such as finite elements methods EFCAD [START_REF] Bastos | Electromagnetic Modeling By Finite Element Methods[END_REF] or ANSYS for example. Some differences between the analytical and numerical values are denoted concerning the electromagnetical torque and then the optimal solution found by the rational methodology proposed in [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF] must be adjusted.

An in-depth analysis shows that the problem of comparing analytical and numerical results is a very complicated one. Indeed, the general analytical model is based on some restrictive assumptions which are taken into account in order to develop its equations. This model comes from the electromechanical conversion and the flux conservation by assuming that the magnetic induction in the airgap is purely radial. The respective permeabilities of magnets and iron are fixed as unity and infinity. Firstly, analytical models for non-slotted machines are developed and then they were extended to slotted machines thanks to the function K S (which gives current electric loading, with b e = 1) and to the Carter coefficient k c , see subsection II-B and [START_REF] Fitan | A general analytical model of electrical permanent magnet machine dedicated to optimal design[END_REF], [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF].

The magnetic flux density computation using finite element methods [START_REF] Sadowski | Finite element torque calculation in electrical machines while considering the movement[END_REF], is more accurate than the analytical one. Nevertheless, we must do some other assumptions. At a design stage, waveforms of the flux and the feeding currents of the electrical machines are assumed to be ideal: rectangular, trapezoidal or sinusoidal. So the performances or the characteristics of electrical machines can be deduced from flux computations. For instance, for a permanent magnet machine, the no-load flux in windings due to magnets (Φ 0 ) and the flux in windings for two types of load currents (longitudinal and transversal, which give the longitudinal and transversal inductances L d and L q ) are computed. From these three values, the torque, flux and voltage can be calculated for any type of sinusoidal currents. The electromagnetic torque can be expressed as follows:

NUMT(D, L, ...) = 3p(Φ 0 Icosψ - L d -L q 2 I 2 sin2ψ) ( 8
)
where I is the circuit current and ψ the phase angle difference between the current and the electromotive force.

To make easier the validation phase, we developed a numerical tool, named NUMT, [START_REF] Lefèvre | Building a CAD system for educational purpose based only on a mesh tool and a FE solver[END_REF]. This algorithm is able to translate the values of parameters issued from IBBA or given by the user and then NUMT draws and meshes automatically the corresponding machine. The meshing is performed using simple laws which divided the different regions of drawing in a well adapted way, before calling "Triangle" a free 2D mesh generator [START_REF] Shewchuk | Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator[END_REF]; two examples of so-obtained meshes are shown on Figure 1. The flux computations follow and can be performed with the drawing of the machine. This tool is very interesting to validate our analytical global optima, see [START_REF] Lefèvre | Building a CAD system for educational purpose based only on a mesh tool and a FE solver[END_REF] for details.

IV. THE NEW METHODOLOGY OF DESIGN

The purpose of this work is to answer to the question: is it possible to replace in Problem (1), the first strong equality constraint Γ em (...) = Γ (corresponding to h 1 (...) = 0) by NUMT(...) = Γ?

           min x∈I R n r ,z∈I N n e , σ∈ Q n c i=1 K i ,b∈B n b f (x, z, σ, b) g i (x, z, σ, b) ≤ 0 ∀i ∈ {1, . . . , n g } h j (x, z, σ, b) = 0 ∀j ∈ {2, . . . , n h } NUMT(x, z, σ, b) = Γ (9)
i.e. find the solution which satisfies the value of the torque by using a finite element method. Such a constraint NUMT(x, z, σ, b) = Γ is named a black box constraint because it depends on an algorithm for computing it. Such a new problem ( 9) is impossible to be solved using IBBA because for using interval analysis tools, all the expressions of the objective or constraint functions must be given explicitly. Indeed nowadays, in our knowledge, it is impossible to compute bounds (in a polynomial time) using interval analysis or other tools for the function NUMT over a box of the initial domain of research. Note that for the equality constraints the index j starts from 2. This is due to the fact that comparing to (1) the first equality constraint corresponding to the torque

(Γ em (x, z, σ, b) = Γ) is deleted and replaced by the numerical one NUMT(x, z, σ, b) = Γ.
Therefore, is it possible to consider the more interesting inverse problem of design than (1), which is defined as follows?

               min x∈I R n r ,z∈I N n e , σ∈ Q n c i=1 K i ,b∈B n b f (x, z, σ, b) g i (x, z, σ, b) ≤ 0 ∀i ∈ {1, . . . , n g } h j (x, z, σ, b) = 0 ∀j ∈ {2, . . . , n h } (1 -pc) × Γ ≤ Γ em (x, z, σ, b) ≤ (1 + pc ) × Γ NUMT(x, z, σ, b) = Γ (10)
The purpose of this work is to extend the code IBBA by introducing some steps of NUMT in order to solve problems of type [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF] applied to rotating machines with permanent magnets. Indeed, the combinatorial analytical model permits to lead such an algorithm to the determination of the deterministic global optima. The idea is to find a solution which satisfies numerically in place of analytically the equality constraint of the torque. The analytical computations of the torque is used to determine the domain where some numerical evaluations must be performed:

(1 -pc) × Γ ≤ Γ em (x, z, σ, b) ≤ (1 + pc ) × Γ,
where pc and pc are real values in [0, 1[ which permit to define the domain of research. Therefore, each optimal solution which is found by using the combination of IBBA and NUMT, named IBBA+NUMT below, satisfies numerically the equality constraint on the fixed torque: NUMT(x, z, σ, b) = Γ. The so-obtained solutions are the exact global ones of Problem [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF]; attention must be paid for the user of the algorithm IBBA+NUMT to the definition of pc and pc (during the following numerical experiments pc and pc are fixed to 0.1). In fact, the analytical model is just used to determine a small zone where the numerical solution is sought for, hence if this zone is too reduced then the true numerical optima of problem ( 9) cannot be reached and if the zone is too large the algorithm could not converge.

In the following, one recalls and extend IBBA algorithm published in [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF] in order to solve more general inverse problems of type [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF].

Algorithm IBBA+NUMT:

1) Set X := the initial domain in which the global minimum is sought for,

X ⊆ IR n r × IN n e × nc i=1 K i × B n b . 2) Set f := +∞. 3) Set L := (+∞, X). 4) Extract from L the lowest lower bound. 5) Bisect the considered box chosen by its midpoint, yielding V 1 , V 2 . 6) For j:=1 to 2 do a) Compute v j := lb(f, V j ) (a lower bound of f over V j ). b)
Compute all the lower and upper bounds of all the analytical constraints on V j ; deduction steps using the analytical constraints permit to reduce V j , [START_REF] Messine | Deterministic Global Optimization using Interval Constraint Propagation Techniques[END_REF]. c) if f ≥ v j and no analytical constraint is unsatisfied then

• insert (v j , V j ) in L. • set m the midpoint of V j • if m satisfies all the
analytical constraints and then if the numerical constraint NUMT(x, z, σ, b) = Γ is also satisfied then f := min( f, f (m)).

• if f is changed then remove from L all (z, Z) where z > f and set ỹ := m.

7) If f -min (z,Z)∈L z < (where z = lb(f, Z)) then STOP. Else GoTo Step 4.
we call the analytical constraints, all the constraints excepted the last one NUMT(x, z, σ, b) = Γ. Because the algorithm stops when the global minimum is sufficiently accurate less than , one -global numerical solution is reached: ỹ corresponding to f. However, it can be possible that a better solution exists in the sub-boxes remaining in the list L at the end of the algorithm. Nevertheless, even in the best case, its corresponding minimal value will not be less than f -. Therefore, by correctly fixing the value, the so-obtained solution (ỹ, f) is sufficiently interesting and can be considered as the global solution (more precisely the -global solution) of the considered problem [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF]. For details on IBBA, the way to bisect a box, to compute bounds, to propagate the constraints, to stop the algorithm, etc., see [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF]- [START_REF] Keafott | Rigorous Global Search: Continuous Problems[END_REF], [START_REF] Messine | Deterministic Global Optimization using Interval Constraint Propagation Techniques[END_REF], [START_REF] Messine | A Deterministic Global Optimization Algorithm for Design Problems, chapter in Essays and Surveys in Global Optimization[END_REF], [START_REF] Ratschek | New Computer Methods for Global Optimization[END_REF].

The way to define correctly the parameters pc and pc is not so easy. In this first study, we consider that variations about 10% around the analytical value of the torque are significant enough; i.e. pc = pc = 0.1.

V. EXAMPLES OF DESIGN

The global optimization algorithms used here, IBBA and IBBA+NUMT, are implemented in Fortran 90/95. All computations were performed on an isolated PC computer from our laboratory with 1.8Ghz, hard disk drive with 7200 rpm and 512Mb of RAM.

To illustrate the use of our new algorithm IBBA+NUMT, let us consider the following optimal design problems of a general electrical slotted rotating three-phase machine with permanent magnets. Some parameters such as the diameter D, the length L, and the thickness of the magnets l a are variable. They are all listed in Table I. And some other parameters, such as the current density, the kind of materials and waveform, or the winding fitting factor are fixed for this study. The mechanical airgap is also fixed to the value of 1 mm; if this value were free then it always decreases to its lower bound when we minimize some volumes or the mass. The magnets are a modern NdFeB and the magnetic circuits (yoke and teeth) are made in stamping.

Three problems are solved, the first two ones deal with the minimization of a single criterion, the global volume (V g ) and the mass (M a ), and the third one is a multicriteria problem (M ulti) given by : M ulti = M a 2.9 + V g 9.0 10 -4 [START_REF] Messine | An Interval Branch and Bound Method Dedicated to the Optimal Design of Piezoelectric Actuators[END_REF] where the weight factors are the inverses of the rounded optima which are obtained during the minimizations of the global volume and of the mass, see Table I. The equality constraint is true when the torque is equal to 10 (+ or -0.2) N.m. For IBBA+NUMT, the zone defined using the analytical model is between 9 and 11 N.m; i.e., pc = pc = 0.1.

The numerical results presented in Table I correspond to the solving of respectively problem (1) in columns named IBBA, and problem [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF] in columns named IBBA+NUMT.

Even if all the solutions have the same structure (8 pairs of pole and inverse rotoric configuration), the behavior of our new algorithm is never the same in the three cases. First, for the minimization of the global volume, we can note that the constraint upon the torque is not satisfied numerically for the solution obtained by IBBA. Considering the optimum given by IBBA+NUMT, only the two geometrical parameter D and L change to reach a correct value for the torque. An expert of the domain would be able to do these adjustments or equivalent ones. Nevertheless, the torque value has increased about 6.9% and the volume only about 2%.

The second minimization is different. The torque value of the solution obtained by IBBA satisfies numerically the equality constraint. The use of IBBA+NUMT algorithm confirms the analytical solution. No new optimum has been found in the defined zone using the numerical constraint. This validates, in such a case, the methodology of design proposed in [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF] combining analytical combinatorial models with IBBA algorithm.

For the multicreteria minimization, we note that the numerical torque of the analytical solution is close to 9 N.m. The changes of some parameters of design (for example, D and L) are not so obvious than those for the minimization of the global volume. In this case, IBBA+NUMT found a new optimum which differs from the analytical one for all the geometrical parameters. Hence, for an expert, this solution (or an equivalent one) should be difficult to obtain. This last case show perfectly the real efficiency of our new methodology of design. We can note that the increase of the numerical torque and the volume is respectively 8% and 6.4% whereas the mass has decreased of 2% compared to IBBA. These two solutions are drawn in Figure 1.

Furthermore, we can notice in Table I that the values of the criterion of the optima obtained using the IBBA methods are better than those found using IBBA+NUMT. However, the numerical value of the torque for a solution obtained using IBBA method does not satisfy the main constraint of the torque. Therefore, the analytical solutions of design are less interesting than the numerical optima. We can also note in Table I that the analytical values are very close to the numerical ones for the solutions obtained by IBBA+NUMT.

Moreover, the computation times are rather short, less than four minutes for the slowest problem which owns an optimisation phase with 49 runs of NUMT.

These first results on small-sized design problems are very encouraging for future works. We can expected to solve more complex problems which will take into account different kinds of material or of waveform. We also should be able to include some thermal considerations by the way of introducing new constraints in problems of type [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF].

VI. CONCLUSION

In this paper, one solves for the first time a more general and interesting inverse problem of design (formulated by [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF]) than those already solved in [START_REF] Fitan | A general analytical model of electrical permanent magnet machine dedicated to optimal design[END_REF], [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF]. The so-obtained solutions using our new algorithm IBBA+NUMT satisfy now the constraint about the torque (which is fixed by the schedule of conditions) by a numerical way. Thus, the so-generated solutions are directly validated numerically. Of course, problems of type [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF] are much more complicated than its corresponding problem considering the analytical equation in place of the numerical constraint [START_REF] Bastos | Electromagnetic Modeling By Finite Element Methods[END_REF]. Thus, in this first work, only smallsized problems has been solved. However, some differences are denoted comparing solutions found for (1) and [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF]. Regarding the CPU-times in Table I, we can expect to solve using IBBA+NUMT more general problems of design with more parameters such as those presented in [START_REF] Fitan | A general analytical model of electrical permanent magnet machine dedicated to optimal design[END_REF], [START_REF] Fitan | The Electromagnetic Actuator Design Problem : A General and Rational Approach[END_REF]. This emphasizes the interest of this first work. In the field of global optimization, this is, in our knowledge, the first time that problems with a black-box constraint are solved by an exact algorithm.

Fig. 1 .

 1 Fig. 1. Meshes of the two solutions corresponding to the minimization of the multicriteria

TABLE I COMPARISONS

 I OF THE PREVIOUS AND THE NEW ALGORITHMS

		Parameters		Min Vg	Min Ma	Min M ulti
	Name	Bounds	Unit	IBBA	IBBA+NUMT	IBBA	IBBA+NUMT	IBBA	IBBA+NUMT
	D	[0.01, 0.3]	m	0.1330	0.1310	0.1400	0.1400	0.1400	0.1330
	L	[0.01, 0.3]	m	0.0474	0.0496	0.0496	0.0496	0.0451	0.0519
	l a E	[0.01, 0.3] [0.005, 0.03]	m m	0.0047 0.0074	0.0047 0.0074	0.0039 0.0074	0.0039 0.0074	0.0039 0.0074	0.0047 0.0075
	C	[0.003, 0.02]	m	0.0049	0.0049	0.0039	0.0039	0.0050	0.0039
	β	[0.7, 0.9]		0.89	0.89	0.74	0.74	0.89	0.74
	k d p	[0.4, 0.6] [[3, 10]]		0.5043 8	0.5043 8	0.4978 8	0.4978 8	0.5043 8	0.5021 8
	m	{1, 2}		1	1	1	1	1	1
	σr	{-1, 1}		-1	-1	-1	-1	-1	-1
		V g M a Multi	m 3 kg	8.881 × 10 -4 3.208	9.063 × 10 -4 3.309	9.716 10 -4 2.939	9.716 10 -4 2.939	9.067 10 -4 3.099 2.0761	9.473 10 -4 3.048 2.1032
	Analytical Torque	N.m	9.810	9.992	9.823	9.823	9.859	9.846
	Numerical Torque	N.m	9.311	9.950	9.802	9.802	9.125	9.831
	Computation time	min	1'36	2'20	2'30	4'40	2'39	3'33
	Numerical Computations		-	13	-	38	-	49


Y Lefèvre
email: lefevre@leei.enseeiht.fr

J Fontchastagner

F Messine

Building a CAD system for educational purpose based only on a mesh tool and a FE solver

the article discusses the use of FEM software to help lecturer to teach design of electrical machines. Nowadays, lecturers have less and less time to teach the design of electrical machines. Special tools based on FEM can help them in this task. In past years, such tools needing sophisticated environments have been realized for designers in industry. Assuming minimum equipment, requirements, guidelines and a realization of a specific tool are presented.

I. INTRODUCTION

In one hand, we know that there is a growing problem to educate students engaged in an electrical engineering curriculum. Electrical engineering students have to achieve a very large spectrum of knowledge in science and technology ranging from fundamental sciences such as mathematics, physics (electromagnetism, material sciences, and mechanics) to applied science like control, data processing, electronic, and power electronics. So students in electrical engineering who want to be expert in electrical machines have less and less time to study what should be their main subject: the design of electrical machines.

In the other hand, we know that industry is lacking in designers of electrical machines because we are witnessing a growing development of the applications of electrical motors and actuators. Industry wants them more and more efficient and well adapted to their applications. So design of electrical machines becomes a more and more difficult task.

We think that finite element method (FEM) should be more and more used in academic area to help lecturers to teach electrical machine design. This paper presents a discussion about the use of CAD system based on finite element code for teaching electrical machines and a realization of such a system to help teaching design.

II. FEM AND TUTORIAL ELECTRODYNAMICS

Commercial finite element (FE) softwares are more and more used in industry. They are more and more introduced in academia. These softwares can be a very good help to teach electromagnetism and FEM which is now one of the most powerful tool in electrical engineering. They are very useful to help to explain many notions such as leakage flux, traveling field, air-gap induction, slots effects. They, especially those who can take into account the movement of the rotor, can also be a valuable help to analyze the electromechanical behavior of classical machines such as synchronous machines, induction machines or reluctance motors [START_REF] Sadowski | Finite element torque calculation in electrical machines while considering the movement[END_REF].

But their main drawback for students is that most of them are not easy to use and students must spend some hours, days or weeks before being familiar with them. It is due to the fact that these softwares are of general purpose. They can solve many types of electromagnetic problem, for several different kinds of actuators or electrical machines. So to perform a simulation it is necessary to do several tasks (drawing, meshing, solving and post-processing) which are not very easy for beginners.

For some tutorials proposed to students, a tool that can do these different tasks automatically has to be provided. For instance, students may be asked to analyze the validity domain of an analytic model for sizing permanent magnet motors. They have not to worry about how to draw, how to mesh, how to solve and how to compute flux. They only have to concentrate on their problem.

In the past, many searchers have developed such a tool for designers in industry. Most of these tools are developed in a special environment depending on some particular FE software. There is for instance a parametric environment which has been developed around OPERA-2D [START_REF] Parker | Parametric Environment for EM Computer Aided Design[END_REF], a design process has been built up around FLUX-2D [START_REF] Wurtz | Methodological Guidelines for the use of Analytical and Numerical Models in a Design Process of an Electromagnetic Device[END_REF] and a way of building an integrated cad system around EFCAD has been attempted [START_REF] Henaux | First steps towards full integrated cad system for electrical machines[END_REF].

III. A TOOL FOR TEACHING ELECTRICAL MACHINE DESIGN

A.Assumptions

We are trying to specify the requirements of such a tool from the point of view of a user with the minimum tools at his disposal. So, we don't try to adapt it to any kind of up-to-date methods which require special environment. The environment considered are common operating system like UNIX, Linux or windows (in fact dos).

We based our discussion by assuming only the existence of a meshing tool that can be called from a language like FORTRAN or C and a module that can solve a 2D magnetic problem formulated in potential vector and compute the flux in windings.

We assume that the meshing tool is capable to mesh quadrangular regions of any kind: contour lines can be arc of circle or segment of straight line.

We assume that the solver has a standard data base. As input, it needs: a table of node coordinates; a table of element location; a table of element material properties; a table informing which winding each element belongs to; a table of nodes with homogenous Dirichlet conditions (Az=0), a table of nodes with odd or even cyclic boundary conditions. And as output, it gives: the third component of vector potential (Az) on each node. The flux is calculated from Az on each node and from the table which relates elements to windings.

B.Requirements

The main requirement of the tool is that it can consider several types of machines: induction motor, synchronous motor, reluctance motor. The motors can have several kinds of geometry configuration (the rotor can be outside or inside the stator), several kinds of structure (the numbers of phases, of pair of poles and slots can vary); several kinds of winding are considered (concentrated or distributed). It can consider machines of any dimensions according to the fact that these dimensions are compatible between themselves and with the type of the machines, their configuration and their structure.

C.From geometric and material data to the mesh.

The first step is to divide the study domain (Fig. 1) in quadrangular regions (not shown) characterized by contour lines, a material property and the number of the winding to which it belong if necessary.

A number of nodes are affected to each line of each region. This affectation is done automatically according to a simple law that can be described like this: the number of nodes on the lines (arc) in the air-gap is equal to twice the number of slots harmonics considered multiplied by the number of slots in the study domain. Then the numbers of nodes on other lines are calculated in function of their distance from the air-gap.

When all the data structures of each region are filled, the mesh tool is called to mesh each region. During this process the tables affected to each element are filled from the physical characteristics of each region. The elements and nodes of all elements are collected together and then renumbered in order to reduce the bandwidth of the stiffness matrix. The homogenous Dirichlet conditions are assigned to the nodes on the frontier of the study domain located at the minimum radius and maximum radius, cyclic boundary conditions that can be odd or even are assigned to those located at the minimum angle and maximum angle.

D.Computation of characteristics and performances

At a design stage, waveforms of the flux and the feeding currents of the electrical machines are assumed to be ideal. So the performances or the characteristics of electrical machines can be deduced from flux computations. For instance, for a permanent magnet machine with sinusoidal waveforms, the no-load flux in windings due to magnets and the flux in windings for two types of load currents (longitudinal and transversal) are computed. From these three values, the torque, flux and voltage can be calculated for any type of sinusoidal currents.

E.Verifications

To allow the students verify their results, the tool can be linked to existing FE softwares. This can be done at different stages with the help of files: drawing files (DXF), pre-mesh files [START_REF] Henaux | First steps towards full integrated cad system for electrical machines[END_REF] describing the regions, mesh files with or without the potential on each node.

The figure 2 shows flux density displayed by our own general purpose FE software. To obtain these results we have only given the sizes of the different parts of the motor and the material data. Our tool computes automatically the vector potentials on nodes and the flux in windings for any load currents.

IV. CONCLUSIONS

In this paper, we have discussed about the use of FEM for help to teach applied electrodynamics. We have presented the requirements of a tool adapted to tutorial designs of electrical machines. Guidelines to build such a tool are presented, a prototype has been realized and the first results obtained have been presented.

Fig. 1 .

 1 Fig. 1. Study domain Fig.2.: Flux density


D Sabino

D Fabre

J S Leontini

D Lo Jacono

Vortex-induced vibration prediction via an impedance criterion

Keywords: vortex shedding, parametric instability, flow-structure interactions

The vortex-induced vibration of a spring-mounted, damped, rigid circular cylinder, immersed in a Newtonian viscous flow a nd c apable o f m oving i n t he direction orthogonal to the unperturbed flow i s i nvestigated f or R eynolds n umbers R e i n the vicinity of the onset of unsteadiness (15 Re 60) using the incompressible linearised Navier-Stokes equations. In a first s tep, w e s olve t he l inear p roblem c onsidering an imposed harmonic motion of the cylinder. Results are interpreted in terms of the mechanical impedance, i.e. the ratio between the vertical force coefficient a nd the cylinder velocity, which is represented as function of the Reynolds number and the driving frequency. Considering the energy transfer between the cylinder and the fluid, we show that impedance results provide a simple criterion allowing the prediction of the onset of instability of the coupled fluid-elastic s tructure c ase. A g lobal stability analysis of the fully coupled fluid/cylinder s ystem i s t hen p erformed. T he instability thresholds obtained by this second approach are found to be in perfect agreement with the predictions of the impedance-based criterion. A theoretical argument, based on asymptotic developments, is then provided to give a prediction of eigenvalues of the coupled problem, as well as to characterise the region of instability beyond the threshold as function of the reduced velocity U * , the dimensionless mass m * and the Reynolds number. The influence o f t he d amping p arameter γ o n t he i nstability region is also explored.

Introduction

Flows over bluff bodies are often encountered in a natural environment as well as in many engineering applications. For the canonical flow past a circular cylinder, it is now well established that the two-dimensional symmetric flow undergoes a global instability at a critical Reynolds number of Re c ≈ 47, via a supercritical Hopf bifurcation [START_REF] Provansal | Bénard-von Kármán instability: transient and forced regimes[END_REF][START_REF] Dušek | A numerical and theoretical study of the first Hopf bifurcation in a cylinder wake[END_REF]. This bifurcation is at the origin of the onset of the time-periodic vortex shedding phenomenon.

One of the consequences of the vortex shedding phenomenon is the dramatic increase of the drag force, due to the low base pressure on the rear of the body. The unsteadiness in the wake may excite resonant structural vibrations, thus causing fatigue failure, early deterioration of the materials and consequent safety problems. These vibrations are driven by the fluctuations in the lift and drag forces acting on the body surfaces. Consequently, a coupling between the motion of the structure and the flow, commonly referred as fluid-structure interaction or flow-induced vibration, must be taken into account when studying the vortex shedding phenomenon. In particular, when the coupling between the structure and the flow leads to the cylinder oscillating with the same frequency as the vortex shedding, high-amplitude vortex-induced vibration (VIV) occurs. Further details about the VIV phenomenon are given in the reviews of [START_REF] Griffin | The vortex-excited resonant vibrations of circular cylinders[END_REF], [START_REF] Parkinson | Phenomena and modelling of flow-induced vibrations of bluff bodies[END_REF] and [START_REF] Williamson | Vortex-induced vibrations[END_REF].

All these features are fundamentally linked to instability mechanisms that arise above a critical Reynolds number, Re c , denoted as a supercritical regime, where the coupled fluid-structure system becomes unstable, leading to an unsteady flow.

This instability was previously studied by [START_REF] Cossu | On the instability of a spring-mounted circular cylinder in a viscous flow at low Reynolds numbers[END_REF] via a linear stability analysis (LSA) of a system where the cylinder motion was constrained to oscillate only transverse to the free stream (1DOF). This study argued that the critical Reynolds number may decrease to less than a half of the Re c associated with the fixed-case cylinder, when the mass ratio parameter is small. Subsequently, [START_REF] Buffoni | Vortex shedding in subcritical conditions[END_REF] carried out an experimental study showing that the unsteady vortex shedding could be triggered at Reynolds numbers as low as Re = 25 by forcing the cylinder displacement in the crosswise direction (1DOF) at certain frequencies and amplitudes. Later, a transverse and in-line vibration (2DOF) was investigated by [START_REF] Mittal | Vortex-induced vibrations at subcritical Re[END_REF], showing that vortex shedding and VIV could be found for Reynolds numbers as low as Re = 20 for certain values of mass ratio and natural frequency of the coupled spring system. Furthermore, similar conclusions were made by [START_REF] Meliga | An asymptotic expansion for the vortex-induced vibrations of a circular cylinder[END_REF] using a 2DOF asymptotic expansion approach and by [START_REF] Zhang | Mechanism of frequency lock-in in vortex-induced vibrations at low Reynolds numbers[END_REF] using a 1DOF reduced-order model approach. Subsequently, a study of the lock-in regime (match of the wake vortex mode and body oscillation frequencies) was conducted by Navrose & Mittal (2016) using a 1DOF global LSA, and the subcritical regime of the 1DOF cylinder was studied by [START_REF] Kou | The lowest Reynolds number of vortex-induced vibrations[END_REF] using a dynamic mode decomposition.

The forced cylinder configuration, where the motion of the body is externally controlled, has also been extensively studied. The pioneering work on this area was first accomplished by [START_REF] Bishop | The lift and drag forces on a circular cylinder oscillating in a flowing fluid[END_REF] and [START_REF] Feng | The measurement of vortex induced effects in flow past stationary and oscillating circular and d-section cylinders[END_REF]. This work was followed by the investigations of [START_REF] Williamson | Vortex formation in the wake of an oscillating cylinder[END_REF] where the 'lock-in' regions are studied. [START_REF] Karniadakis | Frequency selection and asymptotic states in laminar wakes[END_REF] also studied a forced cylinder configuration at low supercritical Reynolds numbers. The connection between the coupled and forced configuration has also been investigated, for example, in [START_REF] Morse | Employing controlled vibrations to predict fluid forces on a cylinder undergoing vortex-induced vibration[END_REF], Morse & Williamson (2009a) and Morse & Williamson (2009b). A review on vortex-induced vibrations of isolated circular cylinders was presented in [START_REF] Bearman | Circular cylinder wakes and vortex-induced vibrations[END_REF].

This study is presented in the framework of hydrodynamic instabilities, with the focus on the characterisation of the flow behaviour by a global linear stability approach at a Reynolds number, Re, near the bifurcation where the flow transits from a steady equilibrium state to an unsteady periodic state. When cylinder displacement is allowed, this phenomenon no longer occurs at Re c ≈ 47. Instead, the bifurcation is expected to occur at lower Re and can possibly be triggered by a different mechanism from the classical one encountered in the fixed cylinder case.

The objective of the present study is to fully characterise this transition for a spring-mounted cylinder and predict it with an impedance criterion, as well as to use an asymptotic development based on the forced case to predict the eigenvalue behaviour of the coupled case for both subcritical and supercritical regimes, i.e. cases at Re below and above the threshold value. Recently, a similar reasoning based on the acoustic impedance was also successfully applied to the oscillating flow through a circular hole in a thick plate [START_REF] Fabre | Acoustic impedance and hydrodynamic instability of the flow through a circular aperture in a thick plate[END_REF].

The cylinder instability is therefore first investigated by analysing the case where the cylinder is forced to oscillate sinusoidally transverse to the free stream direction for Reynolds numbers 15 < Re < 60. The results are then treated in terms of the impedance, the ratio between the vertical force coefficient and the cylinder velocity. The impedance is closely related to the definition of energy transfer between the fluid and the cylinder, enabling the forced and coupled cases to be linked at the marginal energy transfer point and in its vicinity. Following this, the coupled case is analysed, by modelling the 1DOF motion with a damped spring-mass equation. Linear stability analysis is carried out around the steady equilibrium state, and the stability of this state is determined based on the eigenvalues of the linearised problem for several (Re, m * , U * , γ ) arrangements. The results of the coupled case are then compared with the impedance-based predictions from the forced case, both for subcritical and supercritical flow regimes.

Mathematical formulation

2.1. Fluid formulation A Cartesian coordinate system is used with its origin at the centre of the cylinder, which travels with the cylinder motion. The î-axis is considered parallel to the streamwise direction, whereas the  -axis points laterally, as depicted in figure 1(a). The spatial and temporal domains are represented by Ω ∈ R 2 and [0, T], respectively. Following the formulation of [START_REF] Mougin | The generalized Kirchhoff equations and their application to the interaction between a rigid body and an arbitrary time-dependent viscous flow[END_REF], which consists of expressing the absolute velocity u in the relative frame moving with the cylinder, the flow around the cylinder can be described by the non-dimensional laminar incompressible Navier-Stokes equations, defined as

∂ t u + ((u -ζ  ) • ∇)u = -∇p + 2Re -1 ∇ • D(u), (2.1a) ∇ • u = 0, (2.1b) on Ω × [0, T], where D(u) def = 1/2(∇u + ∇u T )
is the rate-of-strain tensor, p the static pressure, ρ f the fluid density and ζ the instantaneous cylinder displacement. The Reynolds number Re is defined as Re = DU ∞ /ν, where D is the cylinder diameter, U ∞ the upstream unperturbed velocity and ν the kinematic viscosity.

The domain limits are depicted in figure 1(b). To take advantage of the problem's symmetry, a half-domain representation is implemented with the proper boundary conditions on the axis in order to obtain the symmetric or anti-symmetric (periodic) solutions. The remaining boundary conditions are posed as

u = -ζ  on Γ cy , (2.2a) u = U ∞ î on Γ in ∪ Γ lat ,
(2.2b)

[-∇p + 2Re -1 D(u)] • n = 0 on Γ out , (2.2c)
with n the outward unit normal to the boundary.

ω is the complex eigenvalue. The base flow q b is the solution of the time-independent version of the Navier-Stokes equations, which are detailed in appendix B. Similarly, the cylinder velocity ζ (t) and the dimensionless force C y (t) will be expressed under the same 'eigenmode form'

ζ (t) = ζ e -iωt + c.c. 2 , C y (t) = Ĉy e -iωt + c.c. 2 . (2.6a,b)
The eigenmode ansatz (2.5) and (2.6) will be used for both the forced and the coupled problems, so that a difference in the ω can be interpreted as follows.

(i) For the forced problem, ω is the frequency of the imposed motion of the cylinder and will generally be considered as real (and noted ω ≡ ω F ). (ii) For the coupled problem, on the other hand, ω is the eigenvalue and is generally expected to be a complex number (ω = ω r + iω i ). Note that the imaginary part ω i corresponds to the amplification rate, and its sign defines whether the solution decays or grows over time, representing a stability criterion of the eigenmode.

Linearised forced problem

For the forced problem, q will be noted qF and the linearised Navier-Stokes equations can be set under the following form:

(A ff -iω F B ff )q F = ζF Y , (2.7)
where A ff and B ff are the fluid and mass operators and the operator Y represents the effect of the cylinder motion on the fluid (which appears both as a volumetric term in the Navier-Stokes equations and in the boundary conditions). The weak form expressions for these operators are given in appendix B.

Once this problem is solved, the vertical force coefficient can be extracted as Ĉy = L(q F ) where L is the 'vertical force coefficient' operator, also defined in appendix B. The mechanical impedance Z(ω F ) can now be defined as the ratio between vertical force coefficient amplitude and the cylinder velocity amplitude

Z(ω F ) = - Ĉy ζ = Ĉy iω F ζ .
(2.8)

As shown in the next sections, the concept of impedance provides a useful way to characterise the problem in both the forced case and the coupled case. Moreover, it provides a link with problems in electricity and acoustics where such a concept is widely used. In particular, the real and imaginary part of the impedance (Z = Z r + iZ i ) can be identified with the resistance and reactance, respectively.

Note that, for the forced case, the mechanical impedance can be formally defined as

Z(ω F ) = -L(q F )[(A ff -iω F B ff )q F ] -1 • Y ,
so the computation of this quantity after discretisation of the operators is straightforward.

Linearised coupled problem

In the coupled problem, the one-way coupling between fluid and the cylinder structure is still represented by (2.7) (with a generally complex ω instead of the real ω F ), but the cylinder displacement ζ is now related to the vertical force coefficient Ĉy through the dynamical equation (2.3). The whole set of equations can now be ˆζ written in eigenvalue form considering the state vector qCP = (q, ζ ˆ, ζ , ) T as follows:

(A -iωB)q CP = 0, (2.9)

with

A =       A ff 0 -Y 0 0 0 1 0 0 0 0 1 2 πm * L - 2π U * 2 - 4πγ U * -1       and B =    B ff 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0    . (2.10a,b)
3. Energy considerations 3.1. Linking energy transfer and impedance Following the work done by Morse & Williamson (2009b) and [START_REF]Lock-in in vortex-induced vibration[END_REF], the net energy transfer from the fluid to the cylinder, over a period T, can be defined as

E transfer = T 0 [C y (t) ζ (t)] dt. (3.1)
Considering a forced problem (with a real ω F ) with displacement ζ (t) and dimensionless force C y (t) given by (2.6), the integration over one period (T = 2π/ω F ) leads to (see appendix A)

E transfer = -π| ζ | 2 ω F Z r (ω F ).
(3.

2)

The property of E transfer being positive for some real value of ω F has been proposed by Morse & Williamson (2009b) as an instability criterion for the fully coupled problem. This criterion is usually expressed in terms of the phase shift φ between the vertical force coefficient and the cylinder displacement

E transfer > 0 ⇐⇒ φ = arg( Ĉy / ζ ) ∈ [0, π]. (3.3)
The concept of impedance introduced in the previous section, allows to express this instability criterion in an alternative form

E transfer > 0 ⇐⇒ Z r (ω F ) < 0 for some real ω F . (3.4)
This second expression makes sense when recalling that the real part of the impedance is analogous to an electrical resistance; indeed the possibility of negative resistance is a well known source of instability in electrical circuits (see [START_REF] Conciauro | Meaning of the negative impedance[END_REF].

Link between forced and coupled problems

Although the criterion discussed in the previous section is physically intuitive, it is not totally rigorous as it uses the impedance computed for some real frequency ω F to predict an instability, which corresponds to a solution of the coupled problem with some complex ω. The situation can be clarified when considering the impedance defined by (2.8) as an analytical function of the complex variable ω. Plugging this definition into the dynamical equation (2.3) leads to the following expression:

-ω 2 - 4πγ iω U * + 2π U * 2 - 2iωZ(ω) πm * = 0. (3.5)
Mathematically, the eigenvalues of the coupled problem (2.9) correspond to the complex zeros of this constitutive relation. The objective of this section is to derive a prediction of these eigenvalues as a function of the most easily computed and physically meaningful values of Z along the real ω-axis (i.e. based on the forced problem).

We first note that a condition for the existence of purely real solutions of the relation (3.5) is that the damping rate is zero (γ = 0) and the impedance is purely imaginary (Z r = 0). This leads to the assumption of the existence of nearly real solutions in the asymptotic limit corresponding to the following hypotheses:

|Z r | |Z i |, (3.6a) γ 1. (3.6b)
Under these hypotheses it is relevant to look for nearly real eigenvalues with the form

ω = ω F + ω (1) , with |ω (1) | |ω F | and ω F ∈ R. (3.7)
Introducing these hypotheses into equation (3.5) and performing a Taylor development in terms of the quantities assumed small leads to

-ω 2 F + 2π U * 2 + 2ω F Z i (ω F ) πm * + -2ω F ω (1) - 4πγ iω F U * + 2ω F πm * - ∂Z ∂ω ω=ω F iω (1) -iZ r (ω F ) + Z i (ω F )ω (1) ω F + h.o.t. = 0, (3.8)
where h.o.t. denotes the higher-order terms neglected in the development. Inspection of this expression leads to the following conclusions.

(i) Firstly, the leading-order term (term in the first line) defines an implicit, and possibly multivalued expression when considering ω F as a function of U * and m * . However, it is explicit and particularly simple when solving for U * in terms of ω F and m * ,

U * = 4π 2 ω F ω F - 2 πm * Z i (ω F )
.

(3.9) (ii) Secondly, the first-order term (term in the second line) leads to an explicit expression for the first-order correction to the frequency ω (1) as follows:

ω (1) = -i Z r (ω F ) + 2π 2 γ m * U * πm * - Z i (ω F ) ω F + i ∂Z ∂ω ω=ω F . (3.10)
The imaginary part of this expression directly provides an estimation of the growth rate as function of γ , m * and U * (using the implicit definition of ω F given by (3.9)).

The condition on the imaginary part of the frequency, Im(ω (1) ) > 0, thus provides an improved version of the instability criterion compared to the one discussed in the previous subsection, equation (3.4). One can note that in the absence of damping (γ = 0), the criterion reduces to Z r (ω F ) < 0 provided that the real part of the denominator is positive. As it will be seen in § 6.3, if this latter term changes sign, the asymptotic prediction may lead to erroneous predictions, but the latter are easily identified as being out of the range of application of the asymptotic theory.

Numerical formulation

The numerical approach is based on a finite element method. All the equations solved in the paper are first rewritten in a variational formulation and then spatially discretised via a Delaunay-Voronoi algorithm, generating an unstructured triangular mesh, adapted with an automatic procedure to a given field. This adaptation relies on the AdaptMesh procedure of the FreeFem++ software (see, for instance [START_REF] Hecht | New development in FreeFem++[END_REF]). An adaptation via the structural sensitivity [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF], calculated via the solutions of the direct and adjoint eigenproblems (see [START_REF] Fabre | A practical review on linear and nonlinear global approaches to flow instabilities[END_REF] was used, in order to adequately discretise the wake region in the rear of the cylinder. The mesh validation for the same configuration is detailed in [START_REF] Fabre | A practical review on linear and nonlinear global approaches to flow instabilities[END_REF].

The unknown velocity and pressure fields (u, p) were spatially discretised using the classical Taylor-Hood basis (P 2 , P 2 , P 1 ) in order to satisfy the Ladyzhenskaya-Babuška-Brezzi (known as LBB) condition (see, for instance, [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]). All the discrete matrices resulting from the projection of the variational formulations onto the basis of finite elements were built with the FreeFem++ software and handled using the MUMPS library, performing a direct LU inversion. The loop over the parameters and data treatment were performed using the Matlab drivers of the StabFem project (https://gitlab.com/stabfem/StabFem).

For the forced regime system described in (2.7), the solution is readily computed by inverting the matrix system and multiply it by the forcing term. On the other hand, the linear stability problem results in a generalised eigenvalue problem, which was solved via a shift-and-invert method, a particular case of the Arnoldi method. This strategy allows eigenvalues to be obtained in the vicinity of some given complex shift. Influence of the choice on the adaptation field (and therefore spatial resolution), confinement due to the domain size and sensitivity to the boundary conditions has already been reported in [START_REF] Fabre | A practical review on linear and nonlinear global approaches to flow instabilities[END_REF].

Linear forced results

Validation

The results for the linear forced case, described in the following, are treated in terms of the real and imaginary parts of the impedance Z(ω F ) as a function of the prescribed oscillation frequency represented by the Strouhal number St F = ω F D/(2πU ∞ ) (see figure 2). For validation the linear forced results are compared to fully nonlinear direct numerical simulations (DNS) performed using a spectral-element code. This method is described in detail in [START_REF] Karniadakis | Spectral/hp Element Methods for Computational Fluid Dynamics[END_REF], and the implementation used here follows that described in [START_REF] Thompson | Three-dimensional instabilities in the wake of a circular cylinder[END_REF]. Previous studies using this code for oscillating and coupled vibrating cylinders include Leontini et al. (2006a), Leontini, Thompson & Hourigan (2006b) and [START_REF] Leontini | The flow-induced vibration of an elliptical cross-section at varying angles of attack[END_REF].

Z i Z i Z r Z r St F St F St F Z i Z i Z i Z r Z r Z r (a) (b) (c)
The main parameter that dictates the validity of the comparison between the linear and nonlinear results is the amplitude of the imposed oscillation, which must be small enough to respect the linear hypothesis. The comparison between different imposed amplitudes for the same range of St F at Re = 25 is presented in figure 2(b), for Z r (ω F ) and Z i (ω F ). As the amplitude decreases, the DNS solution approaches the linear solution. The results for the two smaller amplitudes give almost the same results, confirming that the considered displacement amplitude is small enough. Another two DNS simulations were carried out at Re = 15 and Re = 35 for |ζ F | = 0.002. Further, one can see from the plotted data that as Re increases, the DNS results tend to more accurately match the linear forced results.

Instability thresholds established using linear forced oscillations

A presentation of the main results for the linear forced case can be seen in figures 2 and 3(a). The latter sketches the results in a Nyquist diagram (Z r versus Z i ). Two instability zones are highlighted and discussed in the following paragraphs.

The first interesting phenomenon occurs at Re = Re c1 = 19.915 (shown with the solid green curve, figure 3a) where the sign of Z r first becomes negative at St F,c1 = 0.104, indicating a positive energy transfer from the fluid to the cylinder. The limit condition can be written as

At Re = Re c1 , ∃ St F ∈ R + | Z r (St F ) = 0, (5.1)
and this is represented on the plot as the lighter grey region.

The second interesting phenomenon occurs at Re = Re c2 = 30.349 (shown with the dashed red curve, figure 3a) with the curve that passes exactly through the origin of the Nyquist diagram, indicating that a forcing frequency exists where the impedance is exactly zero. The limit condition for this point can be written as

At Re = Re c2 , ∃ St F ∈ R + | Z(St F ) = 0, (5.2)
and the region where simultaneously Z r is negative and Z i is positive is represented by the darker grey region. The classical Nyquist criterion then states that this value of Re c2 presents a lower bound for values of Re that can produce a positive value of Z i -for values of Re beyond this, the Nyquist curve will encircle the origin in an anticlockwise direction. An upper bound exists for this behaviour, Re c3 , which is explained below. Physically, the condition Z(St F ) = 0 implies that there exists a solution of the linearised problem in which the value of C y,F is zero and therefore, following equation (3.1), the energy transfer between the forced case cylinder and the flow will also be zero.

The third critical phenomenon occurs at Re = Re c3 = 46.766 (shown with the solid purple curve, figure 3a). At this Re the steady base flow becomes hydrodynamically unstable. This means that perturbations are expected to grow even in the absence of any forcing, i.e. at a certain St F one will have Z(St F ) → ∞ and therefore ζ → 0. One observes that the curve tends to infinity, with the limit condition written as

at Re = Re c3 , ∃ St F ∈ R + | Z(St F ) = ∞.
(5.3)

For Re > Re c3 , the flow is predicted to be unsteady, regardless of the frequency and energy transfer characteristics of a small-amplitude harmonic forcing. So, the results presented in the Nyquist diagram of figure 3(a) can be interpreted as follows. For Re Re c1 = 19.915, there is a range of forcing frequencies which results in a positive energy transfer to the cylinder and therefore the flow is predicted to be unstable. For Re Re c2 = 30.349, there is a range of forcing frequencies that produce a positive imaginary component of the impedance Z i within the unstable zone Z r < 0. For Re > Re c3 = 46.766, the flow is unstable regardless of the form of any linear forcing.

Linear results for the coupled problem

The results for the coupled problem, described in the following, are analysed in terms of the eigenvalue ω and its variation with the problem's parameters, (Re, m * , U * , γ ). The value of the damping ratio γ will be considered to be zero throughout the following analysis, with the exception of § 6.6. The impedance-based predictions from the forced case at the leading order and at the first order are compared with the results of the LSA both for subcritical and supercritical regimes.

Eigenmode classification and behaviour

As reported in [START_REF] Cossu | On the instability of a spring-mounted circular cylinder in a viscous flow at low Reynolds numbers[END_REF], [START_REF] Meliga | An asymptotic expansion for the vortex-induced vibrations of a circular cylinder[END_REF] and Navrose & Mittal (2016) (the latter will hereinafter be denoted as NM16), the coupled transverse cylinder motion introduces two complex-conjugate eigenvalues to the physical spectrum and therefore two eigenmodes, referred to as a fluid mode and a structural mode (EM), as proposed in NM16. This is in contrast to the fixed case, where only FM can be found.

The classification of the two modes found in the coupled moving cylinder case is divided according to their nature and relation. This identification is based on the eigenmode characteristics in the limit of very large m * since for this case, the flowinduced vertical force becomes very small compared with the structural ones and the fluid-structure interaction drops to negligible levels [START_REF] Meliga | An asymptotic expansion for the vortex-induced vibrations of a circular cylinder[END_REF].

On one hand, the real part of the eigenvalue associated with the EM tends to the natural angular frequency of the cylinder system, i.e. 2π/U * ≡ ω N and the mode behaves similarly to a cylinder oscillating in vacuo. Further, the imaginary part of the eigenvalue tends to zero, matching this latter behaviour.

On the other hand, the eigenvalue associated with the FM tends to the leading eigenvalue computed for the fixed case, i.e. the mode behaves similarly to the modes in the fixed cylinder case.

However, for relatively low m * the distinction between FM and EM is not possible, since eigenvalues may coalesce for certain combinations of (Re, U * , m * ), especially for lower m * or when Re → 46.7. The two leading eigenmodes do not exhibit a clear distinction in terms of their affiliation to being either a FM or a EM. They are therefore referred as coupled modes and, following NM16 nomenclature, they are denoted as fluid-elastic mode 1 (FEM1) and fluid-elastic mode 2 (FEM2). Further at Re Re c3 , for small U * , the FEM1 resembles to the FM, whereas the FEM2 is similar to the EM. However, for large U * , the characteristics of FEM1 and FEM2 resemble the EM and the FM, respectively. For Re Re c3 , the inverse reasoning can be applied and can easily be confirmed in the figures of the following sections.

Validation

The validation of the results is first done by comparing the results with the fixed-case cylinder, a classical case of hydrodynamic instability. For that purpose, results are obtained for U * → 0. This asymptotic approach constrains the spring stiffness to large values, k → ∞ while m * remains finite. It is first noted that the real part of the eigenvalue associated with the EM tends to infinity, for all Re, as expected for the natural frequency of the structure-only system, 2π/U * . These singularity configurations were intensively studied in [START_REF] Fabre | A quasi-static approach to the stability of the path of heavy bodies falling within a viscous fluid[END_REF]. The attention of the present study is rather focused on the behaviour of the FM. The imaginary and real parts of the eigenvalue associated with the FM show a good agreement when compared with the fixed case results, reported in [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF], as shown in figure 3(b).

The results are also compared with data from NM16 and from Zhang et al. (2015) (hereinafter denoted as ZLYJ15). Figure 4 shows the real and imaginary parts of the eigenvalue variation as a function of U * = [3,11] for Re = 60 and m * = {5, 20}. Real parts are sketched in terms of St LSA = ω r /2π. One can notice that the St LSA associated with the FM remains nearly constant, similar to the St LSA associated with the fixed case, for all values of U * . Further, one notes that the growth rate of the most unstable mode is maximum within the U * range where the frequencies of the two modes are close to each other. This is observed both in this case as in the following ones in the next sections. While the real part presents no difference when comparing the present study and the cited references, the imaginary part of the eigenvalue shows a discrepancy between all sources. Although there is a qualitative match of the results, their values do not agree. However, the present study shows an ω i at U * = 3 that is close to the fixed-case cylinder value for the same Re, depicted in figure 3(b), although this is not verified by the values from NM16 or ZLYJ15. Further, the present LSA results are compared with the impedance-based predictions in the next section, showing a good agreement where the asymptotic theory is valid.

One of the possible causes associated with these differences could be the computational resolution. On one hand, NM16 used a stabilised finite element method based on streamline-upwind/Petrov-Galerkin (known as SUPG) and pressurestabilizing/Petrov-Galerkin (known as PSPG) techniques (Navrose & Mittal 2016). On the other hand, ZLYJ15 used a reduced-order model of the unsteady Navier-Stokes equations in a finite-volume formulation coupled with the structural motion equations. Both methods differ from that implemented here and could be the cause of the discrepancy. Since the influence of the mesh refinement, boundary conditions and domain size was already explored in the present investigation (see [START_REF] Fabre | A practical review on linear and nonlinear global approaches to flow instabilities[END_REF]), the latter does not seem to be the source of error.

It is clear that, for Re = 60, an eigenvalue for which ω i is positive will always exist at any reduced velocity. For the purpose of identifying the stable and unstable zones, a different analysis was carried out at lower Re values which is explored in § § 6.4 and 6.5.

Comparison with impedance-based asymptotic predictions

The eigenvalues computed by resolution of the coupled eigenvalue problem are compared with the asymptotic predictions of § 3.2 in figure 5. Figure 5(a) considers the case Re = 60, m * = 20. In this case, already considered in figure 4(a,c), the FM and the EM are distinct, and the asymptotic prediction accurately follows the eigenvalue of the EM. For Re = 60, m * = 5 (figure 5b), a situation close to that in figure 4(b,d) and characterised by coupled modes FEM1 and FEM2, it is observed that the asymptotic prediction first follows the FEM2 branch and then abruptly switches to the FEM1 branch. The same pattern is observed for Re = 60, m * = 0.1 (figure 5c), but the asymptotic prediction displays an abrupt behaviour at the transition, which does not fit with the stability results. Still, the prediction correctly reproduces the behaviour in the vicinity of the thresholds.

Consider now the situation for Re = 42, a value located below the threshold Re c3 ≈ 46.7. In this case, numerical results only provide a single unstable eigenmode (noted EM or FEM2), and a stable mode (noted FM, when corresponding to the fluid, or FEM1). For m * = 20 (figure 5d), the impedance-based prediction asymptotically fits the EM branch around the linearisation value. Note that in this case, the ω F -U * relation defined by (3.9) is multivalued in the range U * ∈ [9.640 -9.729]. The zoom detail in the figure shows that the asymptotic prediction still correctly fits in the vicinity of the neutral curve but quickly departs from it.

The case Re = 42, m * = 5 (figure 5e) is even more complicated. In this case, the asymptotic prediction is multivalued in the range U * ∈ [10. 478 -16.147] and the corresponding amplification rate crosses the horizontal axis four times (identified with circles in the figure). The results obtained by resolution of the coupled eigenvalue problem are comparatively much simpler, with a single unstable FEM2 branch in the range U * ∈ [5. 304 -11.043]. Note that both these values are correctly predicted by the asymptotic analysis. On the other hand, the thresholds predicted by the asymptotic prediction at U * = 10.539 and U * = 16.110 do not fit with any results. Inspection shows that these values correspond to a change of sign of the real part of the denominator in (3.10), while the numerator Z r is large. Hence, in these ranges the asymptotic prediction is used far away from its starting hypothesis in (3.6), which explains why it does not fit with any physical results. Finally, for Re = 42, m * = 0.1 (figure 5f ) the numerically computed branch FEM2 is unstable in the range U * ∈ [1.270 -∞]. Note that in this case there is no upper threshold for restabilisation. The asymptotic prediction displays the same behaviour, but the predictions does not fit at all with stability results for large U * . Again this range is far from the range of application of the asymptotic theory. Note finally that the asymptotic prediction is multivalued in the range U * = [1. 304 -1.355]. The zoom detail in the plot shows that the threshold obtained at U * = 1.270 (which corresponds to a zero of Z r ) matches with the lower threshold of the stability results, while the two others (which correspond to a change of sign of the real part of the denominator in (3.10)) are incorrect predictions. 6.4. Instability map: Re-U * plane The goal of the following section is to identify the stable and unstable zones in the Re-U * plane, for a fixed mass ratio. By determining the limiting curve where the amplification rate is zero, the region of the parameter space where the cylinder is predicted to vibrate (unstable) can be delineated from the region where the flow is predicted to be steady (stable).

For this purpose, an LSA was carried out over a wide range of reduced velocity U * , for several values of Re and at a fixed m * . The values of U * for which ω i = 0, for different values of Re, for constant values of m * = 50 and m * = 4.73 are first plotted in figures 6(a) and 6(b), respectively. The results are compared with those from [START_REF] Kou | The lowest Reynolds number of vortex-induced vibrations[END_REF] and the impedance-based predictions at the leading order, observing a very good agreement between results. The unstable zone starts at Re c3 for lower U * values and tends to a fixed value when U * → ∞. This asymptotic value tends to decrease as m * decreases. Therefore, for high values of m * , the curve will tend to Re c3 (figure 6a), whereas for low values of m * , the curve tends to Re c2 (figure 6b).

This behaviour is confirmed in figure 7, where the instability threshold is plotted for a very low mass ratio of m * = 0.1. Again, a very close match is observed between the predictions from the linear stability analysis and the impedance-based criterion from the forced case. Also marked on the plot with a dashed line is the value of Re c2 , which shows that the threshold for instability approaches this value as U * increases.

This can be explained by considering that for high masses, the oscillation frequency will be similar to the natural structural frequency, and the imaginary component of the impedance Z i will have little effect on the value of U * (as shown in (3.9)) and on the growth. Conversely for lower m * , Z i will impact the dynamics, hence the threshold Re c2 where the change of sign of Z i plays a role. A similar reasoning can be used to explain the influence of the damping ratio γ presented below in § 6.6.

Lastly, the instability threshold curves for the limit cases m * → ∞ and m * → 0 were computed using the impedance arguments at the leading order and are depicted in figure 8. One can see that for high values of mass ratio, the instability threshold curve tends to be the same in the Re-U * plane (figure 8a). On the other hand, as m * tends to zero, a relation in the plane Re-U * / √ m * is found (figure 8b).

Instability threshold at Re ≈ 20

The following analysis is carried out at lower Re, fixing m * and varying U * . When decoupled, the eigenvalues associated with the FM always present a negative amplification rate. When coupled, the eigenvalues of one of the modes also obeys this observation (namely the FEM1). The objective is then to determine for which Re the EM (or the other coupled mode, FEM2) presents an eigenvalue with a positive amplification rate and the value of U * at which this occurs.

The first value of Re for which a positive amplification rate ω i occurs is Re = Re c1 = 19.946. Note that this value is very close to the value of Re c1 = 19.915 predicted via the impedance-based analysis of the linear forced case presented above in § 5.2. Figure 9(a) shows the amplification rate ω i as a function of U * for a series of values of m * for a fixed value of Re = Re c1 . It is clear that regardless of m * , there is a value of U * at the turning point of each curve where the amplification rate is positive and the flow is predicted to be unstable. It is also clear that the value of U * -here dubbed U * crit -increases with the increase of the mass ratio. The variation of U * crit with m * is compared in figure 9(b) with the impedance-based predictions as defined in (3.9) -the match is almost perfect. Further, the value of U * c1 tends to U * c1 = 9.605 as m * → ∞, which approaches the value predicted for the natural shedding frequency f N of the fixed cylinder. Indeed, from figure 3(b) one can see that f N = 0.71921/(2π) = 0.1145 at Re = 40, which is close to the present critical frequency f * c1 = 1/9.6048 = 0.1041.

6.6. Damping factor effect To this point, all the discussion has been made by supposing γ = 0. This assumption intuitively sets the amplification rate of the eigenmodes to their highest values, i.e. for the same set of (Re, m * , U * ) values, the most linearly unstable mode will be the one with γ = 0. Any value of γ different from zero will damp the system and render it 'more stable'.

Here this effect is explicitly investigated by examining the effect of the damping factor on the first instability threshold, Re c1 . A relationship exists between a massdamping parameter (the product γ m * ) and Re c1 for m * 1. This relation is sketched in figure 10 for several mass ratios. For γ = 0 (or zero damping), all the curves converge to the Re c1 discussed in the previous sections. For high values of γ (or γ m * ), the threshold tends to a different value of Re depending on the mass ratio, in a similar way as shown for increasing m * for zero damping as presented in § 6.4. For m * 1, the threshold tends to the value of the fixed-case cylinder, i.e. Re c1 → Re c3 . As the mass ratio decreases, the asymptotic value of Re c1 also decreases until it reaches the value of Re c2 for m * 1, i.e. Re c1 → Re c2 . In conclusion, it can be said that the damping removes energy from the system, and hence the instability is only triggered at a higher Re.

Conclusions

This paper describes the stability of the flow past an elastically mounted cylinder, focussing on the transition from a steady flow to an unsteady one. Two methods have been described and compared, one using a linear stability analysis of the fully coupled fluid-structure system, and another investigating the energy transfer and its link to impedance when the cylinder motion is harmonically forced.

It has been shown that these two cases are intrinsically linked. The fully coupled system and the harmonically forced system result in the same flow and cylinder motion only when there is zero energy transfer from the flow to the structure and the steady flow is marginally stable. Therefore, the direction of the energy transfer which is linked to the impedance in the forced case can be related to the stability of the flow in the coupled case.

The results show that the instability could be triggered as low as Re c1 ≈ 20 for a specific set of parameters. As Re increases, the unstable zone became wider in terms of U * range. For low mass ratios, this range had a branch tending to infinity when passing Re c2 ≈ 30. For high mass ratio the singularity appears only at Re c3 ≈ 47. These critical values are well-predicted by the impedance-based analysis. The same reasoning is able to explain the variation of the stability threshold with the damping ratio.

In the vicinity of the marginally stable configurations, the results from the coupled eigenvalue problem can be predicted by an asymptotic development at first order, based on the impedance results from the forced case, either at subcritical or supercritical regimes.

FIGURE 2 .

 2 FIGURE 2. Components of the impedance as a function of the imposed oscillation frequency for (a) Re = 15, (b) Re = 25 and (c) Re = 35. Linear forced results are shown with a solid black line (--), as well as DNS results for |ζ F | = 0.2 (dotted-dashed blue line -• -• -), |ζ F | = 0.02 (dashed blue line ---) and |ζ F | = 0.002 (dotted red line • • • • • •). The convergence of the DNS results to the linear results with decreasing amplitude is clear for the Re = 25 case shown in (b). The light and dark grey zones represent the first and second instability zones, respectively, according to the Nyquist arguments.

  FIGURE 3. (a) A Nyquist diagram plotting the real and imaginary components of the impedance for the forced results. Each curve represents a constant Re, St F progresses along the curve. Values of Re are: Re = 19.915 = Re c1 (solid green line --); Re = 30.349 = Re c2 (dashed red line -----); Re = 35 (dotted black line -----); Re = 46.766 = Re c3 (solid violet line --); Re = 55 (dotted-dashed blue line -• -• -). The light and dark grey zones represent the first and second instability zones, respectively. (b) A comparison of the eigenvalues calculated from the coupled case for the fluid mode (FM) for U * → 0 at m * = 1000, and those from the fixed-case cylinder (from Giannetti & Luchini (2007) noted GL07). The grey zone represents the unstable region.

FIGURE 4 .

 4 FIGURE 4. Linear stability results: real and imaginary parts of the eigenvalue ω with U * at Re = 60 for m * = 20 (a,c) and m * = 5 (b,d). The zoom in the frequency plots ((c) and (d)) shows that the curves cross each other for m * = 20 but not for m * = 5. The grey zone represents the unstable region.

FIGURE 5 .

 5 FIGURE 5. Real and imaginary parts of the eigenvalue ω with U * for Re = 60 (a-c) and Re = 42 (d-f ) for the linear stability analysis of the coupled problem and the impedancebased asymptotic predictions based on the forced problem. The grey zone represents the unstable region.

FIGURE 6 .

 6 FIGURE 6. The Re-U * plane distinguishing the region of instability for (a) m * = 50 and (b) m * = 4.73. Results from the current linear stability analysis, the impedance-based predictions, and the coupled cylinder results of Kou et al. (2017) (KZLL17).

FIGURE 7 .

 7 FIGURE 7. Stability threshold for m * = 0.1. Compared are the coupled case results using linear stability analysis and the impedance-based predictions. The grey zone represents the unstable region. The dashed line represents the value of Re c2 , showing that at low m * the instability threshold converges to this value as U * increases.

  FIGURE 8. Impedance-based results on the threshold instability curve for (a) the case of m * → ∞ in the Re-U * plane, and (b) the case of m * → 0 in the Re -U * / √ m * plane. The grey zones represent the unstable regions.

  FIGURE 9. (a) Variation of ω i with U * for different m * at Re c1 = 19.946. The value of m * is noted above curves in the graphic. The grey zone represents the unstable region. (b) Variation of U * c1 with m * at Re c1 19.946, comparing the values extracted from the turning points of the curves generated via the linear stability analysis, and the prediction from the impedance from the forced oscillation case.

  FIGURE 10. The LSA results on the variation of Re c1 with the product γ m * for several values of m * . The grey zone represents the unstable region for m * 1. The horizontal dashed lines indicate the Re c2 and Re c3 .

† Email address for correspondence: diogo.ferreirasabino@imft.fr

Appendix A. Derivation of the impedance energy relation

This appendix shows the deduction of (3.2). The non-dimensional energy transfer (or mechanical work) from the fluid to the cylinder, over a time interval [t 1 -t 2 ] can be defined as

where the (• h ) subscript identifies the sinusoidal component of the instantaneous signal.

Using the former expressions of C y (t) and ζ (t), the energy equation (A 1) can be expanded as

or

If the non-dimensional energy transfer over a time period of the cylinder oscillation T cycle = 2π/ω r is now considered, then t 2 = t 1 + T cycle . If one excludes the exponential growth/decay of the eigenmode, characterising only a harmonic motion, equation (A 4) is reduced to

The latter equation shows the important role of the phase angle φ in the control of the energy transfer between the fluid and the structure. When the coupled cylinder oscillations have achieved a fully developed state, the energy transferred from the fluid to the cylinder matches the energy dissipated by the structural damping. If one considers the structural damping to be zero, energy transfer in the limit cycle is also zero. Thus, the value of φ can either be 0 • or 180 • for the coupled case cylinder (Morse & Williamson 2009b;[START_REF]Lock-in in vortex-induced vibration[END_REF].

The impedance concept, introduced in § 2.3.1 and applied to a general case, reads Z = -Ĉy / ζ . Developing the former formula, one arrives to a real and an imaginary part of the impedance given by

For a harmonic motion, where ω i = 0, one will have ω r ≡ ω F and

Introducing the Z r (ω F ) expression on (A 5), one arrives at the expression for the energy transfer over one cycle of oscillation of the purely harmonic case

where

For this problem, the boundary conditions are defined as

2. Perturbation equations The equations governing the perturbation q are obtained by setting q = q b + (q exp(-iωt) + c.c.)/2 in the Navier-Stokes equations and linearising. The cylinder motion is also taken as ζ = ( ζ exp(-iωt) + c.c.)/2. This leads to

with boundary conditions

3. Weak formulation To set these equations to a form amenable to numerical resolution using a finiteelement method, the test functions u o , p o are introduced, associated with the velocity and pressure. Equations (B 3a) and (B 3b) are multiplied by u o and p o , respectively, and integrated over the domain. The Dirichlet boundary conditions are incorporated by penalisation, i.e. by multiplying u o • û by 1/ and integrating along the corresponding boundaries. The Neumann boundary conditions are naturally taken into account by the finite-element formulation. After integration by parts of the pressure and viscous stress terms, this leads to a weak form of the equations that can be defined, for the forced problem, as

where the bilinear operators A ff and B ff and the linear operator Y are defined as follows:

B ff (q F , q o ) = These operators are also used for defining the coupled problem, as presented in (2.9) and (2.10). This is followed by the spatial discretisation, described in § 4. Finally, the 'vertical force coefficient' operator is defined as 7)