Ion Hazyuk
email: ion.hazyuk@insa-toulouse.fr

Marc Budinger

Florian Sanchez

Jean-Charles Maré

Stéphane Colin

Scaling laws based metamodels for the selection of the cooling strategy of electromechanical actuators in the early design stages

Keywords: natural convection heat transfer coefficient, forced convection, vertical cylinder, design of experiments, motor cooling

The integration of electromechanical actuation systems for aerospace applications requires additional thermal stress considerations during the preliminary design, which were not necessary within the hydraulic technology. The favored models during this design phase are simple, continuous and explicit mathematical expressions to enable easy analytical manipulations or implementations in worksheets or optimization loops. The paper shows how such models can be obtained for the desired applications, illustrating the example of a Limited Angle Torque Motor (LATM) with different cooling strategies. The used method for gathering the models is a regression process on data obtained by finite element simulations. The originality of the method stands in the physical insight of the regression function which is based on scaling laws and dimensional analysis. An example of preliminary design of a LATM illustrates the use of the obtained models for the comparison of two possible architectures intended to minimize the mass of the actuator by implementing different cooling strategies.

(left) Spoiler electromechanical actuator on Boeing 787 from MOOG [START_REF] Gilson | Moog next generation control and actuation[END_REF], and (right) electromechanical thrust vector control of the Vega launcher [START_REF] Descamps | Hi-rel electromechanical thrust vector actuation systems for european unmanned launch vehicles a challenge for the next generation[END_REF] and pneumatic power segments (fig. 1). The more electric aircraft concept illustrates this trend with systems for flight control, braking, engine start, air conditioning, deicing, etc [START_REF] Rosero | Moving towards a more electric aircraft[END_REF]. Aircraft engines can be then optimized while taking advantage of electrical network reliability and ease of maintenance of electrical systems. However, this technology also keeps its own weaknesses, among which are the hard thermal constraints. While within hydraulic technology the hydraulic fluid helps to dissipate the heat generated by the losses, in the case of electric technology the heat remains rather localized. This requires additional, thermal stress, considerations during the design, which may be critical in some cases. For instance, electric motors and their associated converters are often sized on thermal criteria [START_REF] Hendershot | Design of Brushless Permanent-magnet Motors[END_REF]. Electrical network on the aircraft is also strongly thermally constrained, especially if the neighboring airframe is in composite material [START_REF] Long | Future directions in aeronautical composites[END_REF]. Thus, the technological change requires the integration of the new constraints in all the stages of system design.

Necessity of simple algebraic models in preliminary design

As shown in the V-cycle diagram (fig. 2), different types of models are used during the developing process of a mechatronic or more electrical system [START_REF] Auweraer | Virtual engineering at work: the challenges for designing mechatronic products[END_REF]. In the descending phase of the Vcycle, the system is gradually defined and the detail level of the models increases accordingly. During the specific design stage, the various components of the solution concept are fine-detailed and their three-dimensional (3D) simulation models can be built up. The aim of these 3D models is the precise evaluation of the developed components by performing local computations (finite element method, etc). Finally, during the integration phase of the process the various components are assembled to form the solution concept and validate its performance. To do this, the highly detailed 3D simulation models are transformed into, less detailed, macro 0D/1D models to perform physical simulations.

Preliminary design is situated on the descending branch of the V-cycle. During this stage, the designer should be able to narrow down to a single candidate solution out of a bunch of possibilities. This process may involve a number of steps such as the choice of one or several candidate concepts/architectures, rough sizing of the components, multi-criteria optimization of the components, comparison of the different solutions and their evaluation with respect to the initial performance requirements. The aim of this stage is to take strategic decisions early in the project, thus avoiding to waste time and money for the detailed design of solutions that in the end may turn to be unfeasible or far from optimal. Therefore, the following features are crucial for the models used in this stage:

• simple -to gain rapidity and easy handling, i.e. simple algebraic or ordinary differential equations that can be easily implemented in worksheets or optimization loops;

• explicit -to enable analytic manipulations that provide direct access to the main integration parameters (mainly weight and geometrical parameters) disposing of very little input information (usually functional requirements);

• continuous -at least in the range of conventional parameter variation for the given application, which is especially useful for the optimization.

Since the operation of even a single component may touch on several physical fields, its sizing is often performed within an optimization process. During this process, particular attention is granted to various hard constraints, also called design drivers, because they restrain the space of possible solutions. They alone may eliminate some candidate solution, by this also significantly reducing the optimization task. In the case of electrical power transformation, thermal constraints are frequently among the design drivers [START_REF] Hendershot | Design of Brushless Permanent-magnet Motors[END_REF].

Note that today a wide range of tools are available for the simulation of the thermal behavior, which are based on finite element/volume methods. However, all of them require a precisely defined geometry. Consequently, they can be used in the bottom stage of the V-cycle, or even during the ascending phases via model reduction techniques, but they are totally unsuitable for descending phase, i.e. the preliminary design.

For the early design stage, a preferred approach may be the use of analytic/semi-empirical laws elaborated for basic geometries such as vertical/horizontal plane surfaces [START_REF] Armaly | Correlations for laminar mixed convection on vertical, inclined and horizontal flat plates with uniform surface heat flux[END_REF][START_REF] Churchill | A correlation for laminar free convection from a vertical plate[END_REF], cylinders [START_REF] Churchill | Correlating equations for laminar and turbulent free convection from a horizontal cylinder[END_REF][START_REF] Minkowycz | Local non similar solutions for natural convection on a vertical cylinder[END_REF][START_REF] Rohsenow | Handbook of heat transfer[END_REF], cubes [START_REF] Yovanovich | On the Effect of Shape, Aspect Ratio and Orientation upon natural Convection from Isothermal Bodies of Complex Shape[END_REF] etc. Two aspects may favor this approach. Firstly, at this stage the system component geometries can be replaced by very basic ones such as parallelepipeds and cylinders (see fig. 1 right: cylinders for the motors and parallelepiped for the converters). Secondly, these laws can relate some geometrical parameters and physical properties of the component/environment to the global parameters such as thermal flows, temperature differences, heat transfer coefficients, etc. However, the dimensional parameter is usually a single "characteristic length" of the body, which gives models valid for limited ranges of the body aspect ratios. At best, some authors provide tables of coefficients to be used in the formula for different aspect ratios [START_REF] Rohsenow | Handbook of heat transfer[END_REF]. Moreover, the given models are valid only for particular topological and boundary configurations. This is not advantageous in preliminary design because the dimensions are actually the searched/unknown variables. Then, since the model to be used in preliminary design is not the same for different variation ranges of the searched variable, it is impossible to know beforehand which model has to be used.

Another possibility is to use the same technique that is used during the integration phase of the design process, i.e. model reduction [START_REF] Del Barrio | Using model size reduction techniques for thermal control applications in buildings[END_REF]. Here, as shown in fig. 2, the 3D models are transformed in lumped parameter models for the assessment and validation of system specifications. Thus, one can use this knowledge to generate relations between the integration parameters (geometrical dimensions and masses) of interest and the functional parameters (torque, inertia, heat transfer coefficient, etc.). One way to do this is the use of the metamodeling concept [START_REF] Raymond | Response surface methodology : process and product optimization using designed experiments[END_REF][START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF]. It supposes to generate a Design of Experiments (DoE) 1 on a particular case study for which the integration parameters are varied, and then fit a regression model on the obtained data. The problem of the regression models such as polynomial, Gaussian RBF or Kriging is their nonphysical basis which in some cases gives poor predictive capability outside the DoE used for their construction [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF]. Nevertheless, a recent paper has shown that if a physical base is used for the regression, the fitting gives simple relations with acceptable predictive capability outside the DoE [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF]. This technique is therefore used in this paper.

For the sake of clarity, the basic principles of the used technique for the research of the regression expression are restated in the second section. The third section briefly describes the application for which the models will obtained and tested. The fourth section presents the thermal model construction. Finally, an example of the utilization of the obtained models is given in the fifth section.

Scaling laws, dimensional analysis and metamodels

This section describes the fundamentals of a particular form of metamodels: the Scaling-LAW-based MetaModels (SLAWMM), which in our opinion are essential for the understanding of the paper body. For an in-depth description of the regression process, the readers are invited to consult [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF]. 1 Design of experiments is a branch of statistics that deals with the planning of a number of tests to be conducted in order to obtain a good quality models after the regression analysis. Its result is a minimal number of optimally spread configurations for which experiments should be conducted. An in-depth explanation of this tool can be found in [START_REF] Montgomery | Design and Analysis of Experiments[END_REF]

Scaling laws and dimensional analysis

Scaling laws based on dimensional analysis, also called similarity laws or allometric models, have been very successfully used throughout the past decades for solving scientific [START_REF] West | The origin of allometric scaling laws in biology from genomes to ecosystems: towards a quantitative unifying theory of biological structure and organization[END_REF] and engineering problems [START_REF] Budinger | Estimation Models for the Preliminary Design of Electro-Mechanical Actuators[END_REF][START_REF] Mendez | Scaling Laws From Statistical Data and Dimensional Analysis[END_REF][START_REF] Pahl | Engineering Design: A Systematic Approach[END_REF][START_REF] Petersheim | Scaling of hybrid-electric vehicle powertrain components for Hardware-in-the-loop simulation[END_REF] and for presenting results in a compact form. The mathematical form of scaling laws is a power law:

y = kL a ( 1
)
where y is the secondary characteristic to be estimated, L the main characteristic length of the component, and k and a are constants. Equation ( 1) is valid only under material similarity, geometric similarity and uniqueness of the design driver assumptions, as demonstrated in [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF] using the Buckingham theorem [START_REF] Miragliotta | The power of dimensional analysis in production systems design[END_REF][START_REF] Murphy | Similitude in engineering[END_REF][START_REF] Taylor | Dimensional Analysis for engineers[END_REF]. To illustrate the construction and the use of such laws, a simplified example addressing brushless motors is described hereafter. The objective is to evaluate the mass of the motor as a function of the rated torque.

For magnetic aspects, the torque T can be linked to the main characteristic length and current density through

f           T, L, d 1 , d 2 , . . . , d n n+1 , J, B r m           = 0 (2)
which, according to the Buckingham theorem, leads to

f               T JB r L 4 , d 1 L , d 2 L , . . . , d n L n               = 0 (3)
The construction of the first dimensionless number T JB r L 4 is explained in Appendix B. By Assuming geometrical similarity (the aspect ratio is constant during the size change) the following relation is obtained:

T = k 1 JL 4 ( 4
)
where J is the current density, L the length of the motor, d i other geometrical dimensions, B r the remanent induction of the permanent magnet and T the electromagnetic torque. It is assumed here that the main design criterion for the motor is the winding temperature [START_REF] Hendershot | Design of Brushless Permanent-magnet Motors[END_REF]. If convective heat transfer characterized by a constant global heat transfer coefficient h is the dominant thermal phenomenon, the current density J can be linked to the dimensions through:

f           J, L, d 1 , d 2 , . . . , d n n+1 , ρ el , θ, h m           = 0 (5)
which, according to the Buckingham theorem, leads to

f               ρ el J 2 L hθ , d 1 L , d 2 L , . . . , d n L n               = 0 (6)
The construction of the first dimensionless number ρ el J 2 L hθ is explained in Appendix B. Assuming geometrical similarity: Figure 3: Brushless motor masses according to the nominal torque [START_REF] Parvex | NX series Motor Datasheet[END_REF][START_REF] Kollmorgen | [END_REF]. The reference point corresponds to the motor whose data where used in the scaling law (eq. ( 9)) to evaluate the constant k 5 .

J = k 2 L -0.5 (7)
where ρ el the resistivity of the copper and θ the maximal temperature rise for the winding insulation. Combining these two aspects, modeled by eqs. ( 4) and ( 7), one can obtain an estimate of the torque depending on the motor size:

T = k 3 L 3.5 (8)
The motor mass can be estimated assuming material and geometrical similarity:

M = k 4 L 3 = k 5 T 3/3.5 (9)
Figure 3 compares the eq. ( 9) with data from industrial catalogs. It shows that scaling laws can provide good trend for the quantities of interest for such components.

Scaling-LAW-based MetaModels (SLAWMM)

To maintain the physical meaning and the benefits of scaling laws, the regression model will be based on the power form from eq. ( 1). However, to eliminate the assumptions of geometric similarity and uniqueness of design driver, estimation models of the following, more general form will be sought:

y = f (L, π 1 , π 2 , . . . , ) = k(π 1 , π 2 , . . .)L a(π 1 ,π 2 ,...) with π i = d i L ( 10
)
where y is the parameter to be estimated, L is the characteristic length of the component, d i are the secondary geometrical dimensions of the component and π i are the dimensionless numbers representing the aspect ratios.

With geometrical similarity, the aspect ratios are constant (π i = C st ) and the form from eq. ( 10) simplifies into the classic scaling law form, eq. ( 1). The desired form of the model [START_REF] Minkowycz | Local non similar solutions for natural convection on a vertical cylinder[END_REF] lends itself to direct regression less easily than polynomial response surfaces [START_REF] Myers | Response surface methodology: process and product optimization using designed experiments[END_REF], product of pure power laws, radial basis functions [START_REF] Jin | Comparative studies of metamodeling techniques under multiple modeling criteria[END_REF][START_REF] Simpson | Metamodels for computerbased engineering design: survey and recommendations[END_REF], or kriging [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF][START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF] directly on the parameter y. Budinger et al. [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF] proposed an approach to determine the shapes of functions k(π 1 , π 2 , . . .) and a(π 1 , π 2 , . . .) representative of multiplicative and power coefficients.

Budinger et al. [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF] tried to find an expression linking a variable of interest to the main dimensions of a component. Therefore, they defined and used only dimensionless numbers that represent aspect ratios. The present paper, however, attempts to express the variable of interest which depends not only on geometrical dimensions but also on other physical variables and/or properties of the used materials. Therefore, the dimensionless numbers that will be used later in the paper will also include physical properties.

Description of the case study

Electromechanical actuator and its considered cooling strategies

Aircraft flight control surfaces such as ailerons, elevators, flaps, spoilers, etc. are manipulated by actuation systems. They generally employ electric motors somewhere in the power transmission chain. For instance, brushless motors can be used at various levels of the system such as :

• In compressors to generate the hydraulic power for hydromechanical actuators.

• In electro-mechanical actuators to drive the screw-nut system.

• In Direct Drive Valves (DDVs) to control hydrostatic actuators, etc.

This paper focuses on the last type of applications where the used motor is a LATM [START_REF] Bencheikh | Limited angel torque actuator for fine angular positioning[END_REF][START_REF] Tsai | Design and control of a brushless DC limited-angle torque motor with its application to fuel control of smallscale gas turbine engines[END_REF]. An example of such apparatus is illustrated in fig. 4. For the selection of a suitable motor, two main aspects should be considered: (i) the motor must be able to deliver a certain nominal torque, T , and (ii) since the DDV is itself a closed loop system, it must ensure a given bandwidth. In order to achieve the latter, the motor mechanical inertia should be carefully chosen. Actually, the time constant of the LATM dynamical model depends mainly on its mechanical inertia (more exactly it is the ratio between the inertia and the viscous friction). Thus, if its corresponding cutoff frequency (the inverse of the time constant) is smaller than the required closed loop bandwidth, the controller should have larger gains. This will give higher amplitudes for the command which, at best, will lead to larger energy consumption. At worst, the maximal torque of the motor will be smaller than the maximal command (the motor torque will enter in saturation) and the system will not ensure the required closed loop bandwidth.

The nominal torque of the LATM depends on its size and the maximal electrical current density. The maximal current density, in turn, is limited by the maximum winding temperature that prevents the insulation degradation. In other words, the nominal torque of the motor is that one which avoids the overheating of its windings. For ground-based applications, the nominal torque of a motor, provided by constructor catalogs, is generally calculated for an open space at a fixed environmental temperature, usually 25 Figure 4: Illustration of a DDV driven by a LATM [START_REF] Fhl | Excellence in aerospace actuation and support chain management[END_REF] the thermal and mechanical boundary conditions may differ significantly, which make catalog data unsuitable. Furthermore, different cooling strategies may be employed in order to improve the thermal behavior of a motor. The usually accepted cooling strategies in this area on which engineers may count during preliminary design are natural and forced convection. Heat removal by conduction toward the structure via the mounting brackets is generally prohibited. Thus, heat may be removed from the windings either by natural convection via its housing (or via a heat sink), and when this is not efficient enough (or the needed heat sink is too heavy) heat may be removed by forced convection via a fluid circulating through its air gap. In both cases, the estimation of the average heat transfer coefficient is crucial for the preliminary design stage. Moreover, since the mechanical inertia of the motor depends on its geometry, implicitly it is impacted by the cooling capabilities of the motor.

The particularity of this application is that it is used for combustion control of a helicopter turboreactor. Thus, since kerosene is available nearby, if necessary, it may be used for equipment cooling.

Modeling approach

As stated in the previous paragraph, the nominal torque delivered by a LATM is limited by the hot-spot temperature of its windings. Thus, to assist the engineer in the selection of a LATM with an appropriate cooling strategy, one would need one or several expressions relating the nominal torque of the motor to the hot-spot temperature of its windings for a given size and cooling strategy. Hence the need of reduced models to predict the electro-magnetic and thermal phenomena taking place in the motor during its operation. This task is achieved in three steps, as shown in fig. 5. First, electric Joule losses, Q, are evaluated from motor geometry, materials and current density. Actually, in electrical motors there are other sources of heat generation such as mechanical friction and iron losses (hysteresis, Foucault current). However, since LATMs function at small average velocities (common for servo-positioning actuators), these ones are usually low in comparison with electric Joule losses and thus can be neglected. Secondly, a model is established linking the winding temperature rise above the environmental temperature, ∆θ, and the input parameters with the corresponding losses Q. When only natural convection is involved, the environment represents the ambient air while in the case of forced convection, the environment is considered to be the temperature of the cooling fluid. In the latter case, the coolant flow rate is also an input of the thermal model. The output of this model helps to determine if the current density J avoids winding overheating. Lastly, a third model is established to relate the torque to the motor geometry, materials and current density.

Model construction and validation

Figure 6 illustrates the sectional views of the studied LATM. Finite element simulations were performed in COMSOL Multiphysics. For each thermal simulation, the hot-spot temperature of the motor windings was retrieved in order to use it for the regression process. Similarly, the torque was retrieved from electromagnetic simulations.

Electromagnetic modeling

As stated in the previous section, the torque of a motor depends on its geometry, material properties and current density:

f (T L , D, J, B r ) = 0 ( 11
)
where T L is the linear torque, D the motor diameter, J the current density and B r the magnet remanent induction. In the case of LATM geometry, the torque is very little influenced by the magnetic circuit saturation under the effect of the current density J. According to the Buckingham π theorem, the problem can be expressed as:

f (π 1 ) = 0 ( 12
)
where π 1 = T L JB r D 3 , which implies that the motor's electromagnetic torque can be expressed as: In order to obtain the parameter k by regression, finite element simulation were carried out in Comsol. Since we measured the linear torque, a planar electromagnetic model was used (fig. 6 bottom) where the air gap e varied with the same ratio as the motor diameter D (this is an observation that can be made on the Parvex motors range [START_REF] Parvex | NX series Motor Datasheet[END_REF]) and the remanent induction of the magnet was constant (B r = 1 T ). The mesh contained only triangular elements and its independence on the solution has been tested. During the motor size variation, the computational domain always consisted of about 30 000 elements. The value of the coefficient k in eq. ( 13) is obtained by regression on simulation data (k = 0.00418). These simulation data were obtained for a five level full factorial DoE on the current density J and motor diameter D. The variation range of the current density J was 10 6 -10 7 A/m 2 and that of the diameter is given in table 1. Figure 7 compares the electromagnetic torque gathered from simulations with that obtained by eq. ( 13) for two data sets (first set -"DoE" and second set -"Outside DoE" ranges, indicated in table 1); the first set is used for regression and the second set is generated only for validation purpose. As it may be observed, eq. ( 13) predicts the torque based on motor's geometry and current density with satisfactory precision.

T L = kJB r D 3 (13)
Similarly, the losses transformed in heat by Joule effect, are estimated by:

Q = ρ el τJ 2 ( 14
)
where Q is the heat generation rate, ρ el the electrical resistivity of the copper and τ the copper volume fraction in the copper-resin mixture (τ = 0.4 for all the sizes).

Natural convection modeling

A previous study [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF] have shown that if the windings of the LATM are flooded in potting resin, the thermal gradient appears essentially in the convective boundary layer. For usual configurations, the ratio between the convective (R cv ) and conductive R cd thermal resistances is very high: R cv /R cd ≈ 15. Thus, the temperature of the motor body is practically homogeneous. Consequently, it is acceptable to consider that the heat exchange can be characterized only by its convection heat transfer coefficient.

In day-to-day applications, the actuators are attached to a frame or to another component. In the studied application case, it is considered that the LATM is vertically placed on a horizontal flat surface. For the study of the worst case scenario, it is considered that there is no heat exchange between the LATM and the support surface, i.e. adiabatic boundary conditions are assumed on the bottom face of the motor. This is often imposed by the system engineer. If, for a different application, this assumption is invalid, an additional study must be conducted similarly to the one for the forced convection which is described in the next section.

For this study, the finite element model built in Comsol is axisymmetric (fig. 6 top). The simulated motor is placed in a confined space, whose size is ten times larger than the cylinder itself in order to eliminate the impact of confinement degree on the hot-spot temperature. The boundary conditions are a constant heat flux in the LATM windings and a bulk air temperature of 300 K. The following assumptions are considered for modeling and simulation:

• The air flow around the motor is assumed laminar;

• Given the considered range of the environment temperature, the fluid properties are also assumed to be constant.

Actually, when the density variation is small, the Boussinesq approximation can be made. It assumes that the density is constant everywhere in the equations except in the buoyancy term itself. For the latter it is considered that ρρ a = -ρ a β(θθ a ), where ρ and θ are the density and the temperature in the buoyancy term, ρ a and θ a are the ambient density and temperature respectively and β is the air thermal expansion coefficient [START_REF] Cebeci | Physical and Computational Aspects of Convective Heat Transfer[END_REF];

• The copper-resin mixture composing the LATM windings is homogeneous with a copper volume fraction of τ = 40%. Mixture properties are calculated as suggested by [START_REF] Laïd | Extraction de modèles thermiques simplifiés des machines électriques à partir dun calcul du champ de températures[END_REF];

• Ideal solid-solid contacts are considered in the LATM. Thermal contact resistances are thus neglected;

• At this stage, only convective heat transfer is considered.

The radiative heat transfer is considered to be decoupled and thus modeled by the Stefan-Boltzmann law. This assumption was verified and validated by additional simulations.

The mesh in the solid parts contained triangular elements and in the fluid part contained rectangular elements. Mesh independence on the solution has been tested. The studied problem depends on eight dimensional variables:

f (ρ, µ, λ, C p , gβ∆θ, L, D, h) = 0 ( 15
)
where ρ is the air density, µ is the air dynamic viscosity, λ is the air thermal conductivity, C p is the air heat capacity, g is the gravitational acceleration, ∆θ is the winding temperature rise above the bulk air temperature, β is the air thermal expansion coefficient, L is the LATM length, D is the LATM diameter and h is the global convective heat transfer coefficient. The combination gβ∆θ corresponds to the Archimedes force in Boussinesq approximation. Since the problem depends on eight dimensional variables and four independent fundamental quantities (mass, distance, time and temperature), according to the Buckingham theorem the problem can be reformulated in terms of four dimensionless variables:

Nu = f (π 2 , Pr, Gr) (16)
where

π 2 = D L is the aspect ratio, Gr = gβϕρ 2 L 4 µ 2 λ
Grashof number, Nu = hL λ Nusselt number and Pr = λ µC p Prandtl number. Since in this application the convection is generated by the heat produced in the cylinder, the Grashof number is defined in terms of the heat flow rate density ϕ instead of temperature difference ∆θ (via the Fourier's law). For the given parametric study, the Prandtl number is constant. Thus, according to [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF], the researched mathematical form of the convective heat transfer coefficient is:

h = k(π 2 , Gr)L a(π 2 ,Gr) (17)
By replacing the expression of the Nusselt number in eq. ( 16) gives: Thus, the power expression is a(π 2 , Gr) = -1 while k(π 2 , Gr) = f (π 2 , Gr) should be identified.

Nu = hl λ = f (π 2 , Gr) (18)
or h = f (π 2 , Gr)λL -1 (19)
The regression process is performed in two steps. First, the mathematical form of the multiplicative expression k(π 2 , Gr) is identified and then the regression itself is performed on the identified bases. To proceed for the regression, a three level full factorial DoE was generated in logarithmic scale, whose parameters range of variation are given in table 1 in the second column. Using this DoE, the data necessary for the regression process were obtained by finite element simulations.

Following the methodology described in [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF], the variation of the multiplicative coefficient k versus π 2 while Gr is constant is plotted on a logarithmic base in fig. 8 left. The same is done for the second dimensionless number Gr in fig. 8 right. Thus fig. 8 shows that the expression of the multiplicative coefficient can be expressed as a power law: k(π 2 , Gr) = Aπ B Gr C . By performing a nonlinear regression on the simulation data (DoE -table 1 second column), the obtained coefficients A, B and C are:

h SLAWMM = 0.61π -0.179 2 Gr 0.19 λL -1 (20)
The results obtained by applying the SLAWMM method are compared with those obtained by two other correlations provided by the literature for a representative geometrical shape. The considered geometry is a suspended vertical cylinder since this shape resembles most to the LATM geometry. Moreover, the same procedure is used by motor manufacturers to study the thermal behavior of their motors [START_REF] Kollmorgen | Limited Angle Torquers (LATs)[END_REF]. The first compared correlation, suggested by Minkowycz and Sparrow [START_REF] Minkowycz | Local non similar solutions for natural convection on a vertical cylinder[END_REF], evaluates the heat transfer coefficient by an analytical study of the boundary layer for different sizes and thermal boundary conditions. The second one, presented in the Handbook of Heat Transfer [START_REF] Rohsenow | Handbook of heat transfer[END_REF], gives an experimental correlation of the convective heat transfer coefficient for a suspended vertical cylinder. Figure 9 compares the convective heat transfer coefficients evaluated by SLAWMM, Minkowycz and Handbook methods. Two data sets were used for the SLAWMM method; one set was used for the regression procedure, and another data set is used for the evaluation of prediction capabilities of the correlation. For the second data set, the variation ranges of the parameters exceed those of the DoE used for regression (refer to the corresponding column 1). By analyzing the results obtained in fig. 9 it can be concluded that Handbook method does not offer sufficiently accurate predictions for parameter variations in the range given in table 1. Minkowycz method offers satisfactory prediction capabilities, except for some particular cases where the precision may be unacceptable (see fig. 9 left, circled areas). These values are obtained for π 2 > 2, that is to say the motor has a large diameter as compared to its height. In practice, such geometry may be particularly interesting when lightweight actuators with large torque are required (see eq. ( 13)). Actually, Minkowycz method is based on the analytic resolution of the convective boundary layer of the lateral surface of the cylinder. Thus, when motor diameter exceeds its height, convection on the top surface of the motor significantly impacts the overall mean heat transfer coefficient. As shown in fig. 10, the larger the aspect ratio π 2 , the higher the relative error. By contrast, in the SLAWMM method the correlation is obtained by regression on simulation data, so the impact of the top surface convection is still captured. Therefore, eq. ( 24) offers better prediction capabilities as compared to Minkowycz method (see fig. 10). Furthermore, in the case of SLAWMM method, the prediction error can be decreased even more by extending the DoE over the regions of interest, while this is impossible to achieve within the Minkowycz method.

The importance of a simple mathematical form of the obtained model in eq. ( 20) lays in its further use for preliminary design. During this phase, the modeling process described in fig. 5 is reversed; engineers tend to impose the maximal temperature rise in the motor winding and calculate the maximal admissible current density, which is then used for torque calculation. The current density is obtained from the Newton law, where the heat transfer coefficient should be used. Thus, since the heat transfer coefficient in eq. ( 20) depends on the heat flow rate (via Grashof number), eq. ( 20) cannot be used directly in this form. Nevertheless, after some simple mathematical manipulations, described in Appendix A, one is enable to express this coefficient in terms of motor geometry and winding temperature rise above 21). It should be noted that if the mathematical form of eq. ( 20) would not be so simple, probably it was more difficult, or even impossible, to obtain eq. [START_REF] Pahl | Engineering Design: A Systematic Approach[END_REF].

h SLAWMM = 1.05D -0.22 L -0.0753 ∆θ 0.234 (21)

Internal forced convection modeling

For this case study the same LATM is considered except for the fact that a fluid passes through the air gap of the motor to cool it form inside. The used fluid is kerosene -the same one used by the hydraulic jack (for certain applications in the proximity of the reactors). Several simulation tests revealed that even for the smallest flow rates usually encountered in this area of applications, the heat evacuated by natural convection toward the ambiance is insignificant as compared to the heat evacuated by the coolant fluid. Moreover, in the range of the studied geometrical dimensions and for the usually encountered flow rates, the hydraulic and thermal entry lengths are small (max. 10%) as compared to the size of the duct (L). Therefore, it can be considered that the flow is laminar and fully developed thermally and hydrodynamically. Therefore, for this configuration, adiabatic limit condition is assumed on the external surface of the LATM (see fig. 11). Given all these assumptions, the problem depends on nine dimensional variables:

f (ρ, µ, λ, C p , U, L, D, e, h) = 0 ( 22
)
where ρ -kerosene density, µ -kerosene dynamic viscosity, λ -kerosene thermal conductivity, C p -kerosene heat capacity, U -kerosene velocity along the lateral surface of the LATM, L -cylinder's length, D -cylinder's diameter, e -the air gap and h -convective heat transfer coefficient. Since the problem depends on nine dimensional variables and four independent fundamental quantities, it can be reformulated in terms of five dimensionless variables:

Nu = f (π 2 , π 3 , Pr, Re) (23)
where π 2 = D L and π 3 = e D aspect ratio, Nu and Pr Nusselt and Prandtl numbers defined as previously and Re = ρUe µ the Reynolds number. Since the air gap of the motor varies with the same rate as the diameter, the form factor π 3 is constant and thus it can be omitted from the analysis. Assuming small variations of the fluid temperature, the same applies to the Prandtl number. The dimensional analysis which gave these dimensionless numbers, is conducted by considering only the heat exchange between the lateral surface of the motor with the fluid. The heat exchanged between the fluid and the top/bottom surfaces of the motor are neglected at this stage. According to SLAWMM methodology the researched mathematical shape of the heat transfer coefficient is:

h = k(π 2 , Re)L a(π 2 ,Re) (24)
For the same reasoning as in the previous, natural convection, case the power expression is a(π 2 , Re) = -1, which gives:

h = k(π 2 , Re)λL -1 (25)
Applying [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF], the variation of the multiplicative coefficient k versus one dimensionless number while the other one is constant, is plotted in fig. 12. In this case, its mathematical form is a polynomial:

k(π 2 , Re) = A 0 + A 1 π 2 + A 2 π 2 2 + B 1 Re + B 2 Re 2 + Cπ 2 Re (26)
By performing the regression itself, the obtained coefficients are:

h SLAWMM = (77 -42.795π 2 + 2.487π 2 2 + 5.876Re -0.070Re 2 - -0.295π 2 Re)λL -1 (27)
Figure 13 compares the convective heat transfer coefficient obtained by eq. ( 27) with the one obtained by the Newton law, for two data sets. The first data set was used for the regression itself (DoE) and the second one is used for validation purpose (outside the DoE). It can be noted that there is a group of estimates (circled points in fig. 13) that are quite far from the simulated ones. This group is obtained for geometries with π 2 = 2.6. When the aspect ratio becomes larger, the estimation error increases even more. Since the dimensional analysis was conducted for the heat exchange between the lateral surface of the motor and the fluid, this explains the validity limits of the obtained correlation. That is to say, once the motor diameter is larger than its high, the heat transferred by the top and bottom surfaces becomes significant. In order to extend the validity range of the model for larger aspect ratios, the authors presume that the dimensional analysis should also consider the heat exchange via the top/bottom surfaces of the motor. This will give additional dimensionless numbers that potentially may help to minimize the estimation error for the usage beyond the DoE used during the regression.

Metamodel use for the selection of the cooling strategy

This section presents an example of the utilization of previously obtained metamodels in preliminary design of a LATM with its cooling strategy intended for the command of a DDV in aerospace applications. The usual requirements for such applications are the peak and/or Root Mean Square (RMS) torques, the maximal stroke and the actuator bandwidth. Given the application area, the mass of the actuator is also regarded. The peak torque of a permanent magnet motor is usually limited by the demagnetization, its RMS torque is limited by the overheating of the winding insulation [START_REF] Roos | Optimal selection of motor and gearhead in mechatronic applications[END_REF] and its bandwidth -mainly by the rotor mechanical inertia. Here, the RMS torque is considered to be the design driver because of the severe thermal environment (next to the reactor), while the mass and the inertia are surveyed

h Newton [W/m 2 K] h SLAWMM [W/m 2 K]

Outside DoE DoE

Figure 13: Comparison of the convective heat coefficient for the internal forced convection to remain below the acceptable limits. Thus, for the application the following input data and objectives are considered:

• the required RMS torque is 0.1 Nm

• the maximal winding temperature is 140 • C

• the environment temperature is 90 • C

• the LATM mass M should be smaller than 200 -250 g

• the LATM inertia I should be smaller than 10 -6 kg m 2

For the estimation of the mass and inertia of the LATM, the following scaling laws are used [START_REF] Budinger | Estimation models for the preliminary design of electromechanical actuators[END_REF]:

M = kD 2 L (28)
I = kD 4 L (29)
These laws can be used for two distinct cases:

• To check if a given motor will resist overheating in a different thermal environment from than which is specified in the catalog

• To find the dimensions of a motor for a given list of specification (torque, thermal constraints, cooling, etc.) within an optimization loop.

The sizing procedure of the LATM can be resumed in four steps, fig. 14: 1. Heat transfer coefficient calculation: the radiative and convective (natural or forced) heat transfer coefficients are calculated for a given motor dimension (D and L, chosen by the user or the optimization algorithm) and a maximal allowed temperature rise in the motor's winding. 2. Thermal losses calculation: the maximal heat flux density is calculated using the Newton law which gives the corresponding current density for the LATM. Figure 14: Preliminary LATM sizing by optimization and using the constructed metamodels. In the optimization problem, T > T * represents the inequality constraint and the eq. ( 28) represents the objective function to be minimized 3. Torque calculation: the RMS torque of the LATM is calculated using the previously obtained scaling law. If the obtained torque is lower than the required one, the optimization repeats steps 1-3. 4. Mass and inertia calculation: the motor's weight and/or rotor's inertia are calculated, and according to the optimization criteria the steps 1-4 may be repeated or not.

During the first try, with the LATM cooled by natural convection, when the LATM mass was optimized, the sizing gave a motor with the following characteristics: D = 5.56 cm, L = 2 cm, m = 252 g and I = 9.58 • 10 -6 kg m 2 . The obtained mass is acceptable but the inertia is too large. A second try, where the inertia of the LATM rotor was optimized gave D = 2.84 cm, L = 15 cm, m = 492 g and I = 4.9 • 10 -6 kg m 2 . The rotor's inertia is smaller but still too important while the motor is too heavy and its aspect ratio D/L is unrealistic (π = 0.16). Thus, a solution to minimize the motor's mass and inertia is to change the cooling strategy by passing a fluid in the air gap (internal forced convection). The considered fluid is kerosene at 90 • C and a flow rate of 3 • 10 -3 l/s. The minimization of the mass gave D = 2.9 cm, L = 2 cm, m = 68 g and I = 0.7 • 10 -6 kg m 2 . In this case all the design requirements are fulfilled.

Even if a motor with exact sizes obtained in the last example does not exist in catalogs, the aim is to give to the designer a rough idea about the achievable performances when different cooling architectures are used. For this example, the designer can see that by changing the cooling strategy, the motor can be more than three times lighter and its inertia can be about ten times smaller.

Conclusions

This paper attempts to apply the technique of scaling laws based metamodel construction to obtain thermal models adapted for preliminary design of electromechanical actuation systems. An adapted model for preliminary design should relate the integration variables to the other variables of interest via simple, explicit and continuous expressions. Given that during preliminary design the geometry of the components may be viewed as basic geometrical shapes (rectangles and cylinders), some existing semi-empirical models were tested. Although these models gave acceptable estimations for boundary conditions where they were defined, the results may diverge for slightly different topological configurations. For instance, the expression for the estimation of the convective heat transfer coefficient, applicable to a suspended cylinder, don't give accurate estimations for a cylinder, which is three times larger than its high, sitting on an adiabatic surface. Moreover, the mathematical expressions of the semi-empirical models are often laborious which may burden their integration in worksheets or optimization loops (although from scientific point of view this is not an issue, it may count in industrial environment).

The method presented in the paper is based on scaling laws metamodeling, which by their physical fundamentals provide simple mathematical functions while presenting interesting extrapolation capacities. Moreover, the advantage of this method is that it enables to adapt the researched model to a topology requested by the application, for which no anterior studies are available. In comparison with the previous study of the authors, Scaling-LAW-based MetaModels (SLAWMM) [START_REF] Budinger | Scaling-law-based metamodels for the sizing of mechatronic systems[END_REF], where only aspect ratios of the studied component are varied, this paper extends the modeling to cases where other physical variables may also be varied.

The models derived in the paper give the expression of the convective heat transfer coefficient for a Limited Angle Torque Motor (LATM) when it is cooled by the surrounding air or by a fluid passing through its air gap. The obtained expressions are simple enough to be easily handled in preliminary design stage, yet predicting with sufficient accuracy outside the Design of Experiments (DoE) range where they where constructed.

Finally, an example is given to illustrate the SLAWMM use for the selection and preliminary design of a LATM with its cooling strategy. The most important outcome of the obtained models is that they permit to the designer to have rapid rough estimations of the design when comparing several system architectures.

Figure 1 :

 1 Figure 1: Examples of electrical power transformation in aerospace applications:(left) Spoiler electromechanical actuator on Boeing 787 from MOOG[START_REF] Gilson | Moog next generation control and actuation[END_REF], and (right) electromechanical thrust vector control of the Vega launcher[START_REF] Descamps | Hi-rel electromechanical thrust vector actuation systems for european unmanned launch vehicles a challenge for the next generation[END_REF]

Figure 2 :

 2 Figure 2: V-cycle diagram for the design of a mechatronic system and its associated modeling

Figure 5 :

 5 Figure 5: N-squared diagram of the LATM modeling process

Figure 6 :

 6 Figure 6: Sectional views of the studied LATM

Figure 7 :

 7 Figure 7: Comparison of the electromagnetic torques obtained by Comsol simulation, T simulation and eq. (13), T estimate

Figure 8 :

 8 Figure 8: Plots in logarithmic bases of the multiplicative function k versus Grashof number (left) and form factor (π 2 ) (right)

Figure 10 :

 10 Figure 10: Impact of the aspect ratio on estimation precision of Minkowycz and SLAWMM methods

Figure 11 :

 11 Figure 11: Simulation model for the LATM in the case of internal cooling

Figure 9 :

 9 Figure 9: Comparison of the convective heat transfer coefficient for the LATM using Minkowycz (left) Handbook (center) and SLAWMM (right) methods

8 Figure 12 :

 812 Figure 12: Plots in logarithmic bases of the multiplicative function k versus Reynolds number (left) and form factor (π 2 ) (right)

Table 1 :

 1 Range of parameter variations for suspended cylinder and LATM study

		DoE range	Outside DoE range
	Variable	(for regression)	(for model validation)
	L(cm)	2 -6	2 -10
	D(cm)	2.5 -10	2.5 -12.5
	π 2	0.8 -2	0.8 -8
	Gr	2 • 10 6 -4 • 10 6 2.7 • 10 4 -6 • 10 8
	Re	12.9 -25.8	5 -50

Table B .

 B 2: Construction of the dimensionless number from eq.[START_REF] Long | Future directions in aeronautical composites[END_REF]

	Quantity	Symbol		Dimensions
			kg m	s	A
	Torque	T	1	2	-2	0
	Length	L	0	1	0	0
	Current density	J	0 -2	0	1
	Induction of PM	B r	1	0	-2 -1
		T B r T B r J T B r JL 4	0 0 0	2 4 0	0 0 0	1 0 0
	Table B.3: Construction of the dimensionless number from eq. (6)
	Quantity	Symbol		Dimensions
				kg m	s	K	A
	Current density		J	0 -2	0	0	1
	Length		L	0	1	0	0	0
	Electric resistivity		ρ el	1	3	-3	0	-2
	Temperature		θ	0	0	0	1	0
	heat transfer coefficient	h	1	0	-3 -1	0
			ρ el h ρ el J 2 h ρ el J 2 hθ ρ el J 2 L hθ	0 0 -1 3 0 -1 0 0	0 0 0 0	1 1 0 0	-2 0 0 0

Acknowledgements

This work has been partly funded by the European ACTU-ATION 2015 project (Modular Electro Mechanical Actuators for ACARE 2020 aircraft and helicopters) and by the French SYRENA project (SYstme de REgulation Nouvelles Architectures).

Appendix A.

This section illustrates how the convective heat transfer coefficient of a cylindrical shape, expressed and calculated in terms of heat flux rate density, may be expressed in terms of temperature difference. The general expression of the convective heat transfer coefficient, obtained in section 4.2, is

where π 2 = D L and Gr = ρ 2 gβϕL 4 µ 2 λ . As it can be noted, it depends on the heat flow rate density ϕ. According to the Newton law:

where Q is the heat flux to be dissipated and S the motor external surface. Equation (A.2) can also be expressed as:

By replacing ϕ from eq. (A.3) in the Grashof number and in eq. (A.1) gives:

(A.4) Thus, the convective heat transfer coefficient can be expressed from eq. (A.4) as:

This section illustrates how are obtained the dimensionless numbers in eqs. ( 3) and ( 6). These numbers can be obtained by several methods; we used the method described in [START_REF] Taylor | Dimensional Analysis for engineers[END_REF].

The procedure which helps to obtain the dimensionless numbers is based on the use of a table where are listed the variables (quantities) that may construct the dimensionless number, along with their dimensions (usually expressed in SI base units: length, mass, time, electric current, temperature, amount of substance and luminous intensity). Thus, for the construction of the dimensionless number from eq. ( 3), table B.2 can be constructed. In the first part of this table are listed the physical variables involved in the problem, with their associated dimensions. For instance, the dimension of the torque is Nm = kg m s 2 m = kg m 2 s -2 . The exponent of each dimension is noted in the corresponding column. This step is repeated for every variable. Since dimensionless numbers don't have dimensions, the idea is to find how can be combined the physical variables (divided or multiplied) so that the sum of the exponents for each dimension gives zero. These combinations are presented in the second part of the table B.2. We can see that in order to cancel the kg from the torque, only the induction of the permanent magnet can be used. Thus the first group that can be formed is T B r , which has the dimension m 2 A. In order to cancel the A from the first group only current density can be used. The second group is then T B r J , whose dimension is m 4 . Finally, in order to cancel the m from the second group we can use the length at forth power, which gives T B r JL 4 . All the dimensions have been canceled, so the last obtained group is a valid dimensionless number. Exactly the same procedure is followed in table B.3 to obtain the dimensionless number from the eq. ( 6).


Marc Budinger

Valérie Pommier-Budinger

Gael Napias

Arthur Costa Da Silva

Arthur Costa Da Silva

Ultrasonic Ice Protection Systems: Analytical and Numerical Models for Architecture Tradeoff

Keywords: de-icing systems, electromechanical actuator, ultrasonic actuator, piezoelectric transducer, pre-stressed Langevin actuator

come   L'archive ouverte pluridisciplinaire

INTRODUCTION AND STATE OF THE ART

Ice protection systems

Ice accretion on aircraft has been a well identified problem since the beginning of the 20 th century. It can lead to decreased lift, increased drag, reduced thrust reduction, and risk of stalling or even engine failure owing to ice ingestion.

Icing occurs both during flight and on the ground. It has led to many aviation accidents, such as the Air Florida Boeing B737 (1982), the American Eagle ATR 72 (1997), and the Air France Airbus A330 (2009). To ensure aircraft safety, regulatory agencies, such as the CAA (Civil Aviation Authority [START_REF]Aircraft Icing Handbook[END_REF]) and the FAA (Federal Aviation Administration [START_REF]Advisory Circular AC 91-74A -Pilot Guide: Flight in Icing Conditions[END_REF]), have established regulations for aircraft anti-icing and de-icing.

Current strategies for anti-icing and de-icing [START_REF]Aircraft Icing Handbook[END_REF] can be chemical, thermal or mechanical -each having different degrees of efficiency or environmental impact. The chemicals used for de-icing (ethylene, propylene glycol or diethylene glycol) can lower the freezing point but require large volumes of fluids and induce environmental issues and premature wear of the treated parts (especially corrosion). Thermal techniques are used for anti-icing and de-icing in flight and on the ground by liquefaction and vaporization of the ice, but require either a large amount of hot air under pressure to be bled from the engine or a large amount of energy to be provided by the electric grid for electro-thermal solutions.

Thermal solutions require around 4 kW/m² and the total amount of power required to de-ice a Boeing 787 with an electrical de-icing system has been estimated at 76 kW [START_REF] Venna | Mechatronic development of self-actuating in-flight deicing systems[END_REF]. Mechanical de-icing systems are low-energy solutions that aim to break the accumulated ice by applying a mechanical pulse or vibrations to the structure to be protected. Pneumatic systems are commonly used for their low cost but have a significant impact on the aerodynamics of the aircraft and require maintenance. Recent efforts to develop electromechanical systems are justified by the potential of such systems in terms of weight, durability and energy savings.

Electromechanical de-icing systems

Electromechanical deicing systems use electromagnetic actuators, piezoelectric actuators or shape memory alloys. These systems are presented in [START_REF]Aircraft Icing Handbook[END_REF] and [START_REF] Myose | Application of Shape Memory Alloys for Leading Edge Deicing[END_REF] and their advantages and drawbacks are analyzed in Table 1.

This article focuses on actuation with piezoelectric technology, and especially on resonant piezoelectric technologies, which have a better power/mass ratio than static solutions. The deeper insight into piezoelectric de-icing systems presented in Table 2 allows a more detailed analysis to be carried out according to the frequency range of the modes activated by the piezoelectric actuators and to the type of piezoelectric actuators.

Ramanathan et al. [START_REF] Ramanathan | De-Icing of Helicopter Blades Using Piezoelectric Actuators[END_REF] proposed the use of ultrasonic shear waves at very high frequency (1 MHz). They performed experiments with piezoelectric patches bonded to an isotropic plate with a layer of ice. The results indicate that the actuators were able to de-ice the aluminum plate by melting the ice at the interface. Kalkowski et al. [START_REF] Kalkowski | Removing surface accretions with piezo-excited high-frequency structural waves[END_REF] analyzed the frequency range for which wave-based technologies efficiently promote the delamination of ice with minimum power requirements.

Venna et al. [START_REF] Venna | In-Flight De-Icing Self-Actuating Wing Structures with Piezoelectric Actuators[END_REF][10] [START_REF] Venna | Mechatronic development of self-actuating in-flight deicing systems[END_REF] [START_REF] Venna | Piezoelectric Transducer Actuated Leading Edge De-Icing with Simultaneous Shear and Impulse Forces[END_REF] used piezoelectric ceramics bonded onto plates and on the inner flat surface of a leading edge structure to excite low frequency modes to delaminate ice (below 1000 Hz). They used analytical and numerical models to identify the first modes for which the shear stress produced in the ice was greater than the shear stress that would theoretically lead to delamination. The average de-icing time varied between 46 s and 280 s and increased as the icing temperature decreased. Palacios [START_REF] Palacios | Instantaneous De-Icing of Freezer Ice via Ultrasonic Actuation[END_REF] analyzed this result and found that the de-icing time seemed to show that the de-icing was more probably caused by thermal effects than by shear stress. S. Struggl et al. [START_REF] Struggl | A basic approach for wing leading deicing by smart structures[END_REF] conducted the same kind of analysis and experiments with piezoelectric ceramics bonded to a plate and on a leading edge structure to excite low frequency modes (below 500 Hz). They also performed tests in an icing research tunnel and the de-icing was successful at a frequency of 307 Hz.

Seppings [START_REF] Seppings | Investigation of Ice Removal From Cooled Metal Surfaces[END_REF] used a stack of thin piezo-electric discs held in compression by a bolt running through the center of the stack and showed that the pre-stressed actuator driven at 20 kHz was more efficient than piezoelectric patches.

Palacios initiated many studies on de-icing systems [START_REF] Palacios | Design, Fabrication, and Testing of an Ultrasonic De-Icing System for Helicopter Rotor blades[END_REF][17] [START_REF] Palacios | Ultrasonic De-Icing of Wind Tunnel Impact Icing[END_REF][19] [START_REF] Overmeyer | Rotating Testing of a Low-Power, Non-Thermal Ultrasonic De-icing System for Helicopter Rotor Blades[END_REF][21] and tested several technologies. For piezoelectric de-icing systems, he used piezoelectric patches to generate ultrasonic shear stress at high frequency (around a few tens of kHz). He performed tests on plates [START_REF] Palacios | Design, Fabrication, and Testing of an Ultrasonic De-Icing System for Helicopter Rotor blades[END_REF][17] and on leading edges [START_REF] Palacios | Ultrasonic De-Icing of Wind Tunnel Impact Icing[END_REF]. At such frequencies, the delamination of the ice was instantaneous. He also tested an original design of a shear tube actuator driven at 300 V and 436 Hz [START_REF] Palacios | Dynamic Analysis and Experimental Testing of Thin-Walled Structures Driven by Shear Tube Actuators[END_REF]. In [START_REF] Overmeyer | Rotating Testing of a Low-Power, Non-Thermal Ultrasonic De-icing System for Helicopter Rotor Blades[END_REF], experiments showed two main failures related to the bonding of the actuators: fracture of the piezoelectric ceramic at the solder joint location and delamination between the ceramic and the host. To tackle this problem, Palacios proposed optimizing the bonding of the supply wire on the ceramic and the bonding between the piezoelectric actuators and the substrate to avoid concentrations of stresses [START_REF] Overmeyer | Actuator Bonding Optimization and System Control of a Rotor Blade Ultrasonic Deicing System[END_REF]. In [START_REF] Tarquini | Investigation of ice shedding properties of superhydrophobic coatings on helicopter blades[END_REF], he also investigated the effect of hydrophobic coating combined with an ultrasonic deicing system and showed that the ice adhesion depended on the substrate roughness.

More recently, Strobl [START_REF] Strobl | Development of a Hybrid Ice Protection System based on Nanostructured Hydrophobic Surfaces[END_REF] used multilayer piezoelectric patches at frequencies around 4 kHz and icephobic coating to delaminate ice instantaneously on polished surfaces with a low supply voltage.

-4 - According to the previous analyses, the main features of piezoelectric de-icing devices are as follows:

-As regards the frequency range, previous studies on piezoelectric de-icing devices have shown that both low frequencies and ultrasonic high frequencies can be excited to generate shear stress and to break ice. However, no studies have been carried out to determine the best frequency range for which structural resonances maximize the shear stress in the ice while minimizing the stress in the substrate and in the actuators. Consequently, a specific study to determine the optimized frequency range is required and is proposed in this article.

-One other issue is the nature of the modes to be excited: both flexural and extensional modes can produce shear stress in a structure. Thus, a specific study to determine the kind of modes that produce maximal shear stress will be useful for the design of the de-icing system architecture.

-Regarding the actuation technology, this article will compare non-pre-stressed piezoelectric transducers (bonded ceramics) and pre-stressed piezoelectric transducers (Langevin transducer) in terms of force-to-density ratio, robustness, ease of integration in a curved leading edge and power consumption. Non-pre-stressed piezoelectric transducers have been widely tested and examples off architectures can be found in [START_REF] Venna | In-Flight De-Icing Self-Actuating Wing Structures with Piezoelectric Actuators[END_REF][15] [START_REF] Palacios | Design, Fabrication, and Testing of an Ultrasonic De-Icing System for Helicopter Rotor blades[END_REF][18] [START_REF] Strobl | Development of a Hybrid Ice Protection System based on Nanostructured Hydrophobic Surfaces[END_REF]. Figure 1 shows two different architectures with a pre-stressed piezoelectric actuator, one exciting flexural modes, the other exciting extensional modes. Both architectures lead to shear stress at the ice/substrate interface. The objectives of this paper are thus:

-to develop analytical models validated by numerical simulations and complemented by measurements in order to perform trend studies highlighting the nature of the resonance modes and the frequency ranges leading to lowconsumption, compact piezoelectric de-icing devices. -6 --to compare piezoelectric transducer-based architectures enabling an interesting compromise to be reached to generate shear stress at the ice/substrate interface and to promote ice delamination and cracking without damaging the actuator or the substrate.

FEASIBILITY OF PIEZOELECTRIC DE-ICING SYSTEMS BASED ON STRESS ANALYSIS

Study case

In this article, the study case is a simple plane rectangular plate of dimensions 290 × 200 × 1.5 mm 3 that is covered by a 2mm thick layer of ice. The panels of an aircraft airframe are usually supported by ribs and stiffeners and the boundary conditions of such panels are close to clamped on all sides. However, the plates will be studied in the pinned boundary condition in this section because this condition enables simple analytical expressions to be formulated to estimate the stresses produced at the ice/substrate interface and because, at the high frequencies used in the study, the deformations in the middle of the plate are similar for pinned and clamped boundary conditions (as illustrated in Figure 2). The mechanical properties of the ice considered for calculation were chosen so as to correspond to glaze ice and were selected among the values presented in [START_REF] Petrovic | Mechanical Properties of Ice and Snow[END_REF]. The ice and aluminum mechanical properties used for the analyses performed in this article are shown in Table 3. Table 3 -Mechanical properties of the study case materials

Stress generation

In order to choose the vibration modes that are the most suitable for de-icing by generation of stresses in a plate, four types of "failure" modes of the ice/substrate interface were considered: failure owing to excessive tensile stress; failure owing to excessive in-plane shear stress; failure owing to excessive out-of-plane shear stress and failure owing to excessive out-of-plane tensile stress. Figure 3 shows a simplified scheme of the mechanisms and the interface failure expected in each case. The values of ice strength are discussed in [START_REF] Scavuzzo | Structural Properties of Impact Ices Accreted on Aircraft Structures[END_REF], [START_REF] Jellinek | Adhesive Properties of Ice -Part II[END_REF], [START_REF] Laforte | Deicing Strains and Stresses of Ice Substrates[END_REF]. The parametric studies carried out by Scavuzzo et al. [START_REF] Scavuzzo | Structural Properties of Impact Ices Accreted on Aircraft Structures[END_REF] and by Jellinek [START_REF] Jellinek | Adhesive Properties of Ice -Part II[END_REF] and the experiments performed by Laforte et al. [START_REF] Laforte | Deicing Strains and Stresses of Ice Substrates[END_REF] indicate that the adhesion strength of ice depends mainly on the roughness of the accretion surface. Scavuzzos's study also suggests that the adhesion resistance to an out-of-plane stress (shear or tensile) depends, but to a lesser degree, on droplet momentum and surface temperature. Laforte's work shows that the adhesion resistance to an in-plane deformation (distortion or elongation) varies with the ice thickness and Loughborough observed dependence of the strength on the nature of the substrate [START_REF] Loughborough | Reduction of Adhesion of Ice to De-icer Surfaces[END_REF]. Finally, Table 4 exposes the values of the properties presented in the references discussed. In this study, the targeted application is the de-icing of aircraft flight control surfaces. For this kind of application, ice delamination by shear is more efficient than breaking by tensile stress. This is why we focus on out-of-plane shear stress generation. For the assessment of the different de-icing system architectures, we will compute the out-of-plane shear stress and ice will be assumed to be de-bonded when the out-of-plane shear stress exceeds the critical adhesive strength of ice (considered to be 1 MPa for numerical applications of this paper). Moreover, to assess de-icing systems, the stress induced in the substrate and in the piezoelectric ceramics will also be computed in order to quantify the risk of damage to the ceramics.

Comparison of resonance modes on 1D models

The plate resonance modes used in the literature (section 1.2) are essentially in-plane extensional modes and out-ofplane flexural modes. They generate stress corresponding to failure modes a) and c) in Figure 3. However, piezoelectric actuators have limited displacement capacity. Consequently the comparison of these modes is made here in terms of stresses generated for a given displacement. To simplify the analysis, several assumptions are used: support and ice are considered as a thin multilayer beam (1D model), mode shapes are assumed to be identical to those of a uniform beam and boundary conditions are simply supported as introduced in section 2.1. Figure 4 shows the beam under study where:

 x is the transverse position along the beam of length a ;

 n is the number of anti-nodes for the mode considered;

  is the pulsation of the mode considered; For flexural modes, the position of the neutral line h n can be obtained by assuming that the tensile force in a section is zero. This results in:

ℎ 𝑛 = 1 2 𝑐 𝑎𝑙𝑢 ℎ 𝑎𝑙𝑢 2 -𝑐 𝑖𝑐𝑒 ℎ 𝑖𝑐𝑒 2 𝑐 𝑎𝑙𝑢 ℎ 𝑎𝑙𝑢 + 𝑐 𝑖𝑐𝑒 ℎ 𝑖𝑐𝑒 (1)
Table 5 synthetizes the analytical equations derived from [START_REF] Timoshenko | Strength of Materials[END_REF] to compute strains and stresses in the ice, in the aluminum substrate and at the ice-substrate interface. The equations of the peak out-of-plane shear stress were obtained by isolating an element of ice of thickness dx subjected to elastic forces, inertial forces and shear forces at the icesubstrate interface. The resonance frequencies were obtained by the Rayleigh method [START_REF] Géradin | Mechanical vibrations: theory and application to structural dynamics[END_REF]

+ 𝐶 𝑖𝑐𝑒 3 ((ℎ 𝑛 + ℎ 𝑖𝑐𝑒 ) 3 -ℎ 𝑛 3 )

Table 5 -1D modes equations

These analytical equations show that:

for extensional modes, the maximum tensile stress in the substrate or the ice is located on the displacement nodes and the maximum out-of-plane shear stress at the ice-substrate interface is located on the displacement antinodes, for flexural modes, the maximum out-of-plane shear stress is located on the displacement nodes.

Figure 5 also compares the extensional and flexural modes, particularly their ability to generate shear stress at the icesubstrate interface. It also gives the tensile stress in the ice and in the aluminum substrate. It shows that, for a given frequency and for a given displacement, the shear stress level at the ice-substrate interface is smaller for extensional than for flexural modes, while the tensile stress in the aluminum substrate is almost the same. Thus, to generate ice delamination while minimizing the displacement of the piezoelectric actuators, it is more interesting to excite flexural modes. This is done for example by architecture a of Figure 1. Another conclusion that can be drawn from Figure 5 is that it is more favorable to work at high frequencies: the higher the frequency, the higher the shear stress. However, the frequency of use may be limited by other criteria such as power supply issues.

Note: the results of Figure 5 were obtained for a 2 mm thick ice layer. The expression of out-of-plane shear stress  xz shows that this stress would be maximum for a null position of the neutral line (h n =0), which implies an optimal ice thickness h ice :

ℎ 𝑖𝑐𝑒 = ℎ 𝑎𝑙𝑢 √ 𝑐 𝑎𝑙𝑢 𝑐 𝑖𝑐𝑒 (2)
Relation ( 2) gives an optimal ice thickness of 4.5 mm for a 1.5 mm thick aluminum substrate. However, for future calculations we chose a 2 mm thick ice layer, less attractive in terms of stress generation, but more realistic for the targeted application (de-icing of aircraft flight control surfaces).

Figure 5 -Comparison of extensional and flexural modes

Shear stress estimation on 2D models for flexural modes

The previous section gave 1D models of beam type structures. In this section, stresses are expressed for 2D plate type geometries. Following the conclusions of the study for 1D models, stresses will be computed only for flexural modes. For this type of movement, each point of the plate is considered to have a vertical displacement ww in its x and y directions such that 𝑤 =w = w(x, y, t). We assume that the displacement w(x, y, t) can be approximated by the analytical solution for flexural modes of a homogeneous plate [START_REF] Géradin | Mechanical vibrations: theory and application to structural dynamics[END_REF]. Table 6 synthetizes the analytical equations to compute the strains, the shear stress at the ice-substrate interface and the tensile stresses in the aluminum substrate and in the ice. n is the number of anti-nodes on the length for the considered mode and m the number of anti-nodes on the width.

Frequency range for piezoelectric de-icing systems

The equations of Table 6 will be used to study the feasibility of de-icing systems with piezoelectric actuators and to highlight the frequency ranges leading to efficient piezoelectric deicing devices. The study will be performed for the study case of section 2.1, with a number of anti-nodes in the length varying between 1 and 25 and a number of anti-nodes in the width varying between 1 and 15.

Figure 6 and Figure 7 show that the shear stress per µm of displacement at the ice/substrate interface and the tensile stress per µm in aluminum or ice increase with the number of anti-nodes and with frequency. One requirement for deicing systems is to remove ice without stressing the structure on which it is deposited. This means that de-icing systems must maximize stress at the ice/substrate interface while minimizing stress in the structure. The ratio of tensile stress in aluminum to shear stress at the ice/substrate interface represented in Figure 8 allows the frequency range to be found for -15 -which the de-icing systems are the most efficient in meeting this requirement and shows that ultrasonic frequencies, higher than 20 kHz, are more favorable.

EVALUATION OF DIFFERENT ARCHITECTURES OF PIEZOLECTRIC DE-ICING SYSTEMS

The previous section highlighted the type of modes (flexural modes) and the frequency range (ultrasonic) of the resonant modes to be excited to lead to ice delamination. This section aims to assess the type of piezoelectric actuators to be used. Two kinds of actuators are compared: patch type actuators, directly bonded on to the structure to be activated, and Langevin pre-stressed actuators in the configuration Figure 1a.

Methodology for evaluating piezoelectric deicing systems

The proposed methodology for evaluating the 2 different architectures of piezoelectric de-icing systems is divided into 2 main phases:

1. Computation of a reduced model (analytically or numerically) of the chosen architecture of the de-icing system connected to the surface with the layer of ice. Computation of the maximal tensile stress per µm within the PZT ceramics and of the shear stress per µm at the ice/substrate interface. All these results are computed for the resonance mode for which the coupling between the piezoelectric actuator and the structure is the best, i.e. for which the required voltage will be the lowest.

2. Computation of the displacement required to generate the minimal stress value and of the voltage that leads to the required displacement. As the piezoelectric de-icing systems are resonant systems, this result depends strongly on the damping of the structure with the ice.

Reduced model of systems with piezoelectric actuators

This section details the computations of the two design methodology phases. These details are extracted from [START_REF] Pommier-Budinger | Sizing optimization of piezoelectric smart structures with meta-modeling techniques for dynamic applications[END_REF].

The reduced model of a structure with piezoelectric actuators can be made for one mode with a mechanical equation and an electrical equation:

        V C Nq q NV Kq q D q M o C s    ( 3
)
where q is the modal displacement, M the modal mass, K the modal stiffness, q c the electrical charge, V the voltage, D s the modal damping,  the modal electromechanical coupling factor, and C o the modal turned-off capacity.

This model can be computed analytically for simple geometries or with multiphysics Finite Element software (such as COMSOL©, ANSYS© or ABAQUS©) that allows calculations with piezoelectric elements. If the computations are carried out with short-circuited piezoelectric patches, V=0 and the equations become:

       θq q NV Kq q D q M C s    (4)
which means that

𝑁 = 𝑞 𝑐 𝑞 (5)
and for the resonance:

𝑁 = 𝑞 𝑐 𝑞 (6)
For structures with low damping, 𝑄 𝑚 ≈ 𝐾 𝜔𝐷 𝑠 [START_REF] Ramanathan | De-Icing of Helicopter Blades Using Piezoelectric Actuators[END_REF] Consequently, if the displacement required to de-bond ice is known for a resonance mode, the voltage that generates this displacement is given by:

𝑄 𝑚 ≈ 𝐾 𝜔𝐷 𝑠 (8)

Design drivers of piezoelectric actuators

The electromechanical coupling factor N describes the ability of the piezoelectric actuators to generate a force on the degree of freedom q for a given resonance mode. Moreover, stresses are proportional to displacements of the degrees of freedom. We chose to compute the electromechanical coupling factor for the point characterized by the largest displacement (in-plane or out-of-plane depending on excited modes). The electromechanical coupling factor may be influenced by: the mode of the piezoelectric ceramic of the transducer (mode 31 for the bonded patch and mode 33 for the Langevin transducer), the dimensions and positioning of the ceramic relative to the nodes and antinodes, and the presence of adhesive.

A piezoelectric transducer may have several operational limits [START_REF] Berlincourt | Behaviour of Piezoelectric Ceramics under Various Environmental and Operation Conditions of Radiating Sonar Transducers[END_REF]: maximum voltage to avoid depolarization by application of an excessive electric field (400 V/mm of ceramic thickness), maximum stress in the ceramics to avoid failure (about 24 MPa in extension for dynamic applications and for non-pre-stressed ceramics) and maximum temperature to avoid reaching the Curie temperature, which leads to depolarization. The thermal limit will not be considered here because of the very low ambient temperatures.

As shown in Equation 30, the displacement at resonance depends on the damping of the structure, characterized by a mechanical quality coefficient. This coefficient is highly dependent on the boundary conditions and on the complexity of the assembly. To assess the different architectures, we assume a value of Q m equal to 100 here, which corresponds to a damping coefficient of 0.5%.

Piezoelectric de-icing system with bonded ceramics

The first configurations studied here are inspired directly from the work of Palacios ([17][18][20] [START_REF] Overmeyer | Actuator Bonding Optimization and System Control of a Rotor Blade Ultrasonic Deicing System[END_REF]) because he experimentally proved their ability to provide instantaneous de-icing even for untreated surfaces (without polishing, without hydrophobic or icephobic coating). The actuator configurations consist of one or more PZT disks glued onto the surface to be de-iced. Palacios explained that ice delamination occurs for frequencies near the first radial extensional mode of the piezoelectric disk.

Two simple configurations will be processed in this section:

 Free axisymmetric disks described in [START_REF] Soltis | Design and testing of an erosion resistant ultrasonic de-icing system for rotorcraft blades[END_REF] (Figure 9a). This study, which allows simple 2D Finite Element analyses because of the symmetry, will help the phenomena involved to be understood on a 2D model and allows the parameters that are negligible for future simulations in 3D (e.g. glue thickness) to be estimated, thus avoiding fine meshes of very thin layers and heavy 3D models.

 Clamped plates with the dimensions given in section 2.1 (Figure 9b). This geometry will allow a comparison with the architecture based on a Langevin transducer, which will be studied in the next section. A volume of ceramic similar to that of the Langevin transducer is chosen in order not to bias the comparison.

Figure 9 -Studied configurations with bonded ceramics (free disk and clamped plate)

We start by the study of the free axisymmetric disks. Figure 10 describes the geometry studied by Soltis under the direction of Palacios [START_REF] Soltis | Design and testing of an erosion resistant ultrasonic de-icing system for rotorcraft blades[END_REF]. We performed studies step by step (PZT ceramic only, aluminum plate only, PZT/aluminum plate, PZT/aluminum plate /ice) to understand the principle of ice delamination.

The study of the piezoelectric ceramic alone shows that the extensional radial frequency of the PZT disk is around 29.9 kHz. As Palacios recommends studying the modes around the extensional radial mode of the PZT disk, extensional and flexural modes around this frequency are preferentially studied.

We first look at the modal shape of these 2 modes. Figure 11(a)(b)(c) shows the modal shapes of the aluminum plate alone (case (a)), the aluminum plate with the piezoelectric disk (case (b)) and the aluminum plate with the piezoelectric disk and ice (case (c)). These modal shapes are represented for the extensional radial frequency close to that of the PZT disk. With a PZT disk bonded on the aluminum plate (with or without ice), the configuration is not symmetric about one plane normal to the axisymmetric axis and the radial mode thus generates both radial and axial displacements. Then Table 8 gives (i) the parameters of the reduced model for the flexural and extensional modes, including the electromechanical force factor, and (ii) the computation results of the shear stresses required to de-bond the ice. The extensional mode has a high electromechanical force factor; the transducer is well-coupled with the structure. For a quality factor of 100, the required displacement to de-bond the ice is 9.22µm, which induces a voltage of 124 V -which is reasonable -and a tensile stress in the aluminum plate of nearly 42 MPawhich is quite high. For the flexural mode, the electromechanical force factor is much smaller but, as the stiffness is also smaller, the required displacement to de-bond the ice is 7.72µm and the required voltage is 88 V. The tensile stress in the aluminum plate is also around 42 MPa, similar to the extensional mode. The next study is the study case (see section 2.1) of this article: a rectangular aluminum plate of dimensions 290x200x1.5 mm 3 . We consider a square piezoelectric ceramic of dimensions 60x60x2.5 mm 3 bonded in the center of the plate, the volume of which is similar to that of the bolt-clamped transducer that will be used in the next section. The extensional radial frequency of the PZT ceramic is around 31.8 kHz. The resonance extensional frequency of the aluminum plate/ice/PZT assembly close to that frequency is 33.2 kHz. For this extensional frequency, Figure 12 shows the in-plane and out-of-plane displacements. In-plane magnitudes are 5.6 times higher than out-of-plane magnitudes. The computed electromechanical force factor is 1.2 N/V. Obtaining sufficient stresses to achieve delamination requires a voltage of about 86 V and also displacement of 21 µm, which generates a very high stress level in the piezoelectric ceramic.

Palacios emphasizes this issue, which leads to the breakdown of ceramics, in numerous works ( [START_REF] Soltis | Design and testing of an erosion resistant ultrasonic de-icing system for rotorcraft blades[END_REF], [START_REF] Overmeyer | Actuator bondline optimization and experimental deicing of a rotor blade ultrasonic deicing system[END_REF]). In order to excite flexural modes of the surfaces to be protected, the architecture of Figure 1(a) is studied here. A boltclamped Langevin transducer (Figure 13) is thus connected to the plate described in section 2.1 through a spacer ring of length 10 mm. The extensional resonance frequency of the transducer alone is 40.44 kHz. For the Langevin transducer/aluminum plate/ice assembly, the resonance frequency for which the delamination of the ice occurs with the lowest supply voltage is 46.857 kHz. For that frequency, the electromechanical force factor is 1.25 N/V. Obtaining sufficient shear stresses to achieve delamination requires a voltage of about 123 V and the tensile stress in the piezoelectric ceramics is then 9.8 MPa, which is much less than for architectures with bonded piezoelectric ceramics.

Comparison of the 2 architectures

It is possible to compare the actuating architectures according to several criteria:  The limitations of use owing to stresses in PZT ceramics. The advantage of the architecture with Langevin transducers is that stresses in the ceramics are only 9.8 MPa, which is 7 times less than for the architecture with bonded ceramics. Moreover, the risk of mechanical failure is much lower for Langevin transducers than for bonded ceramics as their pre-stressed structure enables them to withstand higher stresses in operation.

 Power consumption is mainly due to mechanical losses during resonance and thus to the mechanical energy of elastic deformation. This energy is 7.8 times greater for the architecture with bonded ceramics because it requires much greater deformations (21 microns against 1.83 microns).

 The size of the power electronics is dependent on the power consumption but also on the capacitive energy stored in PZT ceramics. Power devices with good efficiency, such as resonant inverters, have inductors which are all the larger when the capacitive energy is high. This energy is 2.2 times greater for the architecture with bonded ceramics because, even if the voltage is lower, the capacitance is greater.

Experimental validation

Tests were performed on the solution with the Langevin transducer, which seems the most promising if we consider a criterion that minimizes the energies involved and the stresses in the piezoelectric ceramics. A clamped aluminum plate without coating was covered with 2 mm of ice (glaze-type ice obtained in a freezing chamber). The Langevin transducer was then supplied with a voltage of around 46kHz frequency and 150V and 180V amplitude.

Figure 16 shows the results achieved during the experiment. The delamination of the ice started to occur at 150 V and was even more visible at 180 V. The difference between the estimated voltage and the voltage measured during the experiment may come from the uncertainty on the quality factor, the non-linearity and the uncertainties on the shear stress leading to ice delamination (assumed to 1 MPa for computations). interface. These models also highlighted that ultrasonic frequencies over 20 kHz allowed the shear stress to be maximized for a given displacement, thus avoiding fatigue in the substrate and limiting breakdown of the actuators. The modeling of the substrate and the actuators by finite element models highlighted the interest of Langevin type actuators with regard to the resistance to mechanical stress and the mechanical and electrical energies involved in this kind of de-icing system.

Tests validated the efficiency of ice delamination with an architecture based on Langevin type actuators and favoring structural flexural modes. These experiments also confirm the interest of studying such solutions more intensively.

Figure 1 -

 1 Figure 1 -Architecture for a de-icing system with a pre-stressed piezoelectric transducer: (a) Excitation of flexural modes, (b) Excitation of extensional modes

Figure 2 -

 2 Figure 2 -Comparison of two boundary conditions (clamped versus pinned) for the study of a plate at high frequencies: displacements in the middle of the plate are sinusoidal in both cases

Figure 3 -

 3 Figure 3 -Modes of failure of the ice-aluminum interface. a) Failure due to excessive in-plane tensile stress. b) Failure due to excessive in-plane shear stress. c) Failure due to out-of-plane shear stress. d) Failure due to excessive out-of-plane tensile stress

Figure 4 -

 4 Figure 4 -Element of the beam under study

  ..) and tensile (___) Stress (MPa/µm) Stress vs frequency Shear at ice/substrate interface (flexural mode) Tensile in aluminum substrate (flexural mode) Tensile in ice (flexural mode) Shear at ice/substrate interface (extensional mode) Tensile in aluminum substrate (extensional mode) Tensile in ice (extensional mode)

Figure 6 -Figure 7 -FrequencyFigure 8 -

 678 Figure 6 -Flexural resonance frequencies and stresses for the study case versus the number of anti-nodes

Figure 10 -

 10 Figure 10 -Axisymmetric configuration with bonded ceramics (disk)

Figure 11 (

 11 d) shows the modal shape of the flexural frequency close the extensional radial frequency of the PZT disk for the case of aluminum plate with piezoelectric disk and ice. The axial displacements generated by the flexural mode are phasedshifted by a quarter of a wavelength compared to the axial displacements generated by the extensional mode. Such axial displacements generated by the flexural and extensional modes and phased-shifted explain why this configuration excited by a sweep in frequency is favorable to de-ice the entire surface: nodes of the axial displacements for which shear stress are important are spread over the entire surface.

( a )Figure 11 -

 a11 Figure 11 -Axial and radial displacements for axisymmetric configuration

Figure 12 -

 12 Figure 12 -Displacements for the extensional frequency (33.20 kHz) close to radial frequency of the ceramic (rectangular configuration with a bonded piezoelectric ceramic)

Figure 13 -Figure 14 -Figure 15 -

 131415 Figure 13 -Bolt-clamped Langevin transducer and its first extensional resonance mode

Figure 16 -

 16 Figure 16 -Experimental results for the configuration Langevin transducer+spacer+aluminum plate+ice and excitation of the flexural mode of 46.86 kHz frequency

Table 1 -Comparison of electromechanical deicing systems

 1

	Technological solutions	Main advantages		Main drawbacks
		Contactless possible, high	Low force density, electromagnetic pollution,
	Electromagnetic
		displacements possible		size of feeding electronics
		Force generation (good force/mass
	Piezoelectric (PZT)	ratio) even at high frequency,	Brittleness, small displacements
		energy consumption
	Shape Memory Alloy	Force generation in static	Low dynamics, response time, energy
	(SMA)				consumption
	Comparison				Main
	Technological solutions	Previous works	Main drawbacks
	criteria				advantages
	Static, Low Frequency			Power
			[9][10][11][12] [14]	Force of activation
		(Hz)			supply
	Frequency		[16][17] [18][19][20][21]
	Vibration (kHz)			Fatigue
	range		[22][23]
					Energy
		Waves (MHz)	[7][8]
					consumption
					Brittle
	Bonded Piezoelectric	[9][10][11][12][15]
	ceramics (PZT or	[16][17][18]	Easy to	Not easily
		multilayers)	[20][21][21][22][14][23]	implement	adaptable to
	Technology				curved surfaces
					Less brittle,	Frequency linked
	Pre-stressed piezoelectric
			[15]		stress-	to the size of the
		transducer
					resistant	transducer

Table 2 -Comparison of piezoelectric de-icing systems 1.3. Review of previous works, proposed architecture and objectives of present work

 2

Table 4 -Adhesion properties of impact ice to aluminum surface with matte finish and its strength properties ([25],[26],[27])

 4

  using the expression of the strains to estimate the kinetic and potential energies.

	Angular frequency	𝜔 = 𝜋	𝑛 𝑎	√ 𝜌 𝑖𝑐𝑒 ℎ 𝑖𝑐𝑒 + 𝜌 𝑎𝑙𝑢 ℎ 𝑎𝑙𝑢 𝑐 𝑖𝑐𝑒 ℎ 𝑖𝑐𝑒 + 𝑐 𝑎𝑙𝑢 ℎ 𝑎𝑙𝑢	𝜔 = 𝜋 2 (	𝑛 𝑎	)	2	√	𝐸𝐼 𝜌 𝑖𝑐𝑒 ℎ 𝑖𝑐𝑒 +𝜌 𝑎𝑙𝑢 ℎ 𝑎𝑙𝑢	with
																𝐸𝐼 =	𝑐 𝑎𝑙𝑢 3	(ℎ 𝑛 3 + (ℎ 𝑎𝑙𝑢 -ℎ 𝑛 ) 3 )
			Extensional modes			Flexural modes
	Displacement	𝑢(𝑥, 𝑡) = 𝑈(𝑥) 𝑠𝑖𝑛 𝜔𝑡 = 𝑈 0 𝑠𝑖𝑛	𝑛𝜋𝑥 𝑎	𝑠𝑖𝑛 𝜔𝑡	𝑤(𝑥, 𝑡) = 𝑊(𝑥)𝑠𝑖𝑛 𝜔𝑡 = 𝑊 0 𝑠𝑖𝑛	𝑛𝜋𝑥 𝑎	𝑠𝑖𝑛 𝜔𝑡
																𝑢(𝑥, 𝑡) = -𝑧	𝜕𝑤 𝜕𝑥	= 𝑈(𝑥)𝑠𝑖𝑛 𝜔𝑡
	Peak tensile strain	𝜀 𝑥 =	𝜕𝑈(𝑥) 𝜕𝑥	= 𝑈 0	𝑛𝜋 𝑎	𝑐𝑜𝑠	𝑛𝜋𝑥 𝑎	𝜀 𝑥 =	𝜕𝑈(𝑥) 𝜕𝑥	= 𝑊 0 𝑧 (	𝑛𝜋 𝑎	)	2	𝑠𝑖𝑛	𝑛𝜋𝑥 𝑎
	Peak ice tensile	𝜎 𝑥 = 𝑐 𝑖𝑐𝑒 𝜀 𝑥 = 𝑐 𝑖𝑐𝑒 𝑈 0		𝑛𝜋 𝑎	𝑐𝑜𝑠	𝑛𝜋𝑥 𝑎	𝜎 𝑥 = 𝑐 𝑖𝑐𝑒 (ℎ 𝑖𝑐𝑒 + ℎ 𝑛 ) (	𝑛𝜋 𝑎	)	2	𝑊 0 𝑠𝑖𝑛	𝑛𝜋𝑥 𝑎
	stress in ice
	Peak tensile stress	𝜎 𝑥 = 𝑐 𝑎𝑙𝑢 𝜀 𝑥 = 𝑐 𝑖𝑐𝑒 𝑈 0	𝑛𝜋 𝑎	𝑐𝑜𝑠	𝑛𝜋𝑥 𝑎	𝜎 𝑥 = 𝑐 𝑎𝑙𝑢 (ℎ 𝑎𝑙𝑢 -ℎ 𝑛 ) (	𝑛𝜋 𝑎	)	2	𝑊 0 𝑠𝑖𝑛	𝑛𝜋𝑥 𝑎
	in aluminum
	Peak out-of-plane shear stress at the	𝜏 𝑥𝑧 = (𝜌 𝑖𝑐𝑒 𝜔 2 -𝑐 𝑖𝑐𝑒 (	𝑛𝜋 𝑎	)	2	) ℎ 𝑖𝑐𝑒 𝑈 0 𝑠𝑖𝑛	𝑛𝜋𝑥 𝑎	𝜏 𝑥𝑧 = 𝑐 𝑎𝑙𝑢 (	𝑛𝜋 𝑎	)	3 (ℎ 𝑎𝑙𝑢 -ℎ 𝑛 ) 2 -ℎ 𝑛 2 2	𝑊 0 𝑐𝑜𝑠	𝑛𝜋𝑥 𝑎
	ice/aluminum
	interface

Table 6 -2D flexural modes equations

 6

	Flexural modes

These equations are validated by a Finite Element analysis performed for the geometry of the study case (h alu = 1.5 mm, c alu = 70 GPa, h ice = 2 mm, c ice = 9.7 GPa, h n = 0.48 mm, a=290 mm, b = 200 mm). Table

7

shows the mode for n=13 Ultrasonic ice protection systems: analytical models for architecture trade-off -14 -and m=5 of frequencies around 15 kHz, which is an intermediate frequency among the frequencies usually used for tests with piezoelectric ceramics bonded on substrate. The comparison of analytical and numerical results shows a difference of less than 8% and validates the analytical models developed.

Table 7 -Comparison of analytical equations and Finite Element Analysis

 7

Table 8 -Computation results for the axisymmetric configuration in the case of bonded PZT+Aluminum Plate+Ice

 8

	1,2
	1,0				Axial
					Radial
	Normalized	0	0,05	0,1	0,15	0,2
			Radial position (m)

Table 9 -Computation results for the rectangular configuration in the case of bonded piezoelectric ceramic+aluminum plate+ice 3.3. Piezoelectric de-icing system with Langevin transducers

 9

	Computed electromechanical parameters of the reduced model and stresses per µm
	(results of Finite Element Analysis)
				Extensional mode
				at 33.2 kHz
	Modal mass	M	Kg	0.01
	Modal stiffness	K	N/m	4.90E+08
	Modal electromechanical coupling factor	N	N/V	1.2
	Shear stress xz /µm at ice/substratum interface	Txz/µm	MPa/µm	0.05
	Tensile stress xx /µm in PZT	Txx/µm	MPa/µm	3.33
	Ratio Txx/Txz		-	70
	Displacements and voltages required to de-bond ice
	Shear stress required to de-bond ice	Txz	MPa	1
	Displacement required to de-bond ice	U 0	µm	21
	Tensile stress in PZT	T xx PZT	MPa	70

Table 10 -Computation results for the rectangular configuration in the case of Langevin transducer+spacer+aluminum plate+ice

 10

	Computed electromechanical parameters of the reduced model and stresses per µm
	(results of Finite Element Analysis)
				Results at 46.857 kHz
	Modal mass	M	Kg	0.0096
	Modal stiffness	K	N/m	8.4e+09
	Modal electromechanical coupling factor	N	N/V	1.25
	Shear stress xz /µm at ice/substrate interface	T xz /µm	MPa/µm	0.55
	Tensile stress xx /µm in PZT	T xx /µm	MPa/µm	5.3
	Ratio Txx/Txz		-	9.63
	Displacements and voltages required to de-bond ice
	Shear stress required to de-bond ice	T xz	MPa	1
	Displacement required to de-bond ice	W 0	µm	1.83
	Tensile stress in PZT	T xx PZT	MPa	9.8
	Quality factor	Q m	-	100
	Voltage supply	U	V	123