B Montagne
email: benoit.montagne@isae.fr

F Lachaud

E Paroissien
email: eric.paroissien@isae.fr

D Martini

FAILURE ANALYSIS OF COMPOSITE BOLTED JOINTS BY AN EXPERIMENTAL AND A NUMERICAL APPROACH

Keywords: Composite material, bolted joint, finite element analysis

Composite materials are widely used by aircraft manufacturers because of their high mechanical properties. These properties severely decrease when composite material are drilled and bolted. To avoid structural failure, the behavior of composite bolted joints has to be understood and mastered. An experimental database was constituted by Dassault Aviation with supported single-shear tests with countersunk head screws on different composite materials. This database is analysed to find the main parameters that lead to the failure of the joint. The end distance is the most important parameter. It has a strong influence on both the bearing stress at failure and the failure mode. Width and lay-up have a lower effect on joint failure. A two-dimensional finite element (FEM) model with elastic membrane behavior is then developed to simulate the different configurations in the database. Fiber failure criteria is computed and the dispersion of the results for each failure mode is compared to the dispersion of the experimental data to evaluate the quality of the model. The modification of the failure criteria by taking into account different ply orientations and shear stress permits to reduce the dispersion but it is still higher than the dispersion of the experimental data. This work highlights the limits of such linear 2D modelling and it confirms the needs for more complex modelling such as the local bending effects of the joints and/or the non linear behavior of the composite material.

Introduction

Aeronautical structures are made with a lot of pieces assembled together. These joints are mainly done with bolts or rivets. Bolted or riveted joints are easier to tear down than welding or bonding. It permits low cost maintenance operations. Dassault Aviation's aircrafts count hundreds of thousands of fasteners (250 000 on the Falcon 7X and 300 000 on the Rafale). This number could raise to a few millions on the biggest planes such as civil airliners. Nowadays, composite materials are more and more used by aircraft manufacturers. To satisfy the current performance and safety aircraft requirements, it is therefore necessary to master the composite bolted joints behavior. However these joints are difficult to model because of contact and geometrical or material non linearity issues. In a joint, the bolts transfer the load to the different parts. It results in mechanical degradations in the composite material, what could lead to different failure modes of the structure depending on the design parameters such as geometry or lay-up. Therefore, aircraft manufacturers have to understand these modes in order to optimize their design.

In this context Dassault Aviation in collaboration with Institut Clément Ader has started a research project to study composite bolted joint behavior and its modelling. The goal of this project is to define the necessary modelling complexity for characterizing the effects of these design parameters on the failure modes of composite bolted joints. This subject has been studied for many years. Hart-Smith [START_REF] Hart-Smith | Design and analysis of bolted and riveted joints in fibrous composite structures[END_REF][START_REF] Hart-Smith | Mechanically fastened joints for advanved composites -Phenomenological considerations and simple analysis[END_REF][START_REF] Hart-Smith | Bolted joints in graphite-epoxy composites[END_REF] developed a phenomenological approach that resulted in design rules for composites bolted joints. Analytical method based on electrical analogy proposed by Ross [START_REF] Ross | An electrical computer for the solution of shear-lag and bolted joints problem[END_REF] lead to several stiffness formulas [START_REF] Tate | Preliminary investigation of the loads carried by individual bolts in bolted joints[END_REF][START_REF] Swift | Fracture analysis of stiffened structure, damage tolerance of metallic structures : analysis methods and application[END_REF][START_REF] Huth | Influence of fastener flexibility on the prediction of load transfer and fatigue life for multiple row joints[END_REF][START_REF] Paroissien | Analytical Two-Dimensional Model of a Hybrid (Bolted/Bonded) Single-Lap Joint[END_REF] that often need to be adjusted to fit experimental data. Finite element method is another way to study bolted joints. In this paper, an analysis of the experimental database on composite bolted joints provided by Dassault Aviation is presented first for identifying the most important design parameters regarding to bolted joint strength. A two-dimensional FEM is then developed on Samcef v18.1 for improving the understanding of the design parameter influence on the failure of composite bolted joints. This modelling allows an accurate description of the geometry (diameter, width and end distance), the lay-up and the bolt/hole contact. At this point, the composite material behavior is assumed to be elastic and orthotropic.

Experimental results and database analysis

The present database was generated by Dassault Aviation during several experimental studies on composite bolted joints. This database consists in supported single-shear specimens with countersunk head screw with different materials, lays-up, widths and end distances. The specimens are made of composite material jointed with aluminium plates. The countersunk screw head is located in the composite part.In the following analysis the failure stress corresponds to the bearing stress when the specimen load is maximum. The failure modes observed for the different configurations were sorted following the usual classification [START_REF] Hart-Smith | Design and analysis of bolted and riveted joints in fibrous composite structures[END_REF] which is presented in Fig. 2. Theses failure modes are influenced by design parameters such as the end distance (P), the width (2D) or the lay-up.

Influence of the end distance

The bearing stresses at failure are ploted with respect to the end distance P divided by the hole diameter φ in Fig. 3. Three differents T700GC/M21 laminates are studied here. Their proportions of plies oriented in the 4 directions are presented in the following way: 0°/45°/-45°/90°. The shape of the points indicates the failure mode and its color the laminate. The graphs presented in this paper are normalized.

Bearing

Net section Cleavage Cleavage/net section Shearout P D From Figure 3, we observe that specimens fail in bearing mode at a relatively constant stress when the end distance P is large enough (> 2φ). When P decreases (φ ≤P≤2φ) the cleavage failure mode is observed with a failure stress reduction. This reduction is proportional to P. These results were also observed by Collings [START_REF] Collings | On the bearing strength of CFRP laminates[END_REF] and are coherent with industrial criteria. Hart-Smith [START_REF] Hart-Smith | Bolted joints in graphite-epoxy composites[END_REF] also observed constant failure stresses for the bearing failure of T300/N5208 double shear specimens with large P.

It is indicated that if the end distance is enough (P=3φ) for avoiding cleavage failure mode, the net section failure mode could then be observed when D is small as illustrated in figure 3. Hart-Smith [START_REF] Hart-Smith | Design and analysis of bolted and riveted joints in fibrous composite structures[END_REF][START_REF] Hart-Smith | Mechanically fastened joints for advanved composites -Phenomenological considerations and simple analysis[END_REF][START_REF] Hart-Smith | Bolted joints in graphite-epoxy composites[END_REF] advised to use laminates with a repartition of plies in the four directions (0°,45°,-45° and 90°) between 1/8 (12.5%) and 3/8 (37.5%) within the laminates in order to favor bearing failure. The three laminates studied here have the following repartitions of 90°-oriented plies:

-L1 = 45%/20%/20%/15% -L2 = 45%/24%/24%/7% -L3 = 35%/20%/20%/25%

The failure modes of these laminates agree with this design rule: the net section failure mode of L2 laminate leads to a significant reduction in the failure bearing stress compared to the bearing failure mode when it does not for L3 laminate.

Lastly, when the end distance is reduced (φ≤P≤2φ), the variation of the failure bearing stresses of the three laminates is reduced compared to this same variation when 2.5φ≤P≤3φ even if different failure modes are observed when P≤2φ. The pictures on Figure 4

Influence of the width

The Figure 5 shows the evolution of the failure bearing stresses as a function of D (W=2D) for two values of P. This D parameter does not seem to have a significant effect on stresses at failure except for the net section failure mode of the L2 laminate as mentioned before. Modifying D does not change the failure mode when the end distance is short (P=1.5φ) even if the L3 laminates failed by shearout.

When the end distance is larger (P=3φ) the reduction of D leads to net section failures but the failure bearing stresses are almost unchanged for well proportioned laminates.

Influence of the load orientation

The geometrical parameters are set to P=3φ and D=2.5φ for analysing the influence of the load orientation (i.e. the angle between the 0° plies and the direction of loading). In this part, the studied specimens are made in T800S/M21 with two different 50/20/20/10 stacking sequences. According to the Fig. 6 there is no effect of the load orientation on the failure bearing stresses or on the failure mode (bearing) excepted for the 90° angle. In this case, the failure mode is net section for a slightly smaller stress at failure (~10% drop compared to the 0° angle).

Figure 6. Bearing stress at failure according to the angle between the 0° plies and the loading direction

Conclusions

From experimental data analysis we conclude that the most influent design parameter on the joint strength is the end distance P. It has a strong effect on the failure modes and stresses. The D parameter and the loading direction only have a second order influence on bolted joints strength. The same results are obtained for IMS/977-2 and T800S/M21 materials. This conclusion is also in agreement with Hart-Smith's observations [START_REF] Hart-Smith | Design and analysis of bolted and riveted joints in fibrous composite structures[END_REF][START_REF] Hart-Smith | Mechanically fastened joints for advanved composites -Phenomenological considerations and simple analysis[END_REF][START_REF] Hart-Smith | Bolted joints in graphite-epoxy composites[END_REF] on other composite materials and should be therefore a intrinsic characteristics of composite materials bolted joints.

Two-dimensional finite element analysis

A two-dimensional FEM is then developed for improving the comprehension of the effect of the design parameter on the joint strength. This modelling allows a precise description of the geometry of the previsously presented joint configurations and therefore of the associated strain and stress fields in the hole neighbourhood. A first-ply failure analysis is then performed using the Hashin's failure criteria for predicting the differents failure modes presented in this work. The failure of the laminate is expected when at least one of its plies failed. The objective is to reduce the dispersion of the experimental results compared to the previous failure bearing stress analysis. The dispersion is characterized by the coefficient of variation (CoV) which is the standard deviation divided by the mean of the considered values

Model presentation

The FEM is developed on Samcef Software (SIEMENS PLM Software). The configurations are the previously mentioned ones (i.e. T700GC/M21, IMS 977-2 and T800S/M21 materials with different sets of P, D and lays-up). The materials of the finite element model are linear and elastic. The T700GC/M21 mechanical properties are listed in the Table 1. The Young modulus and stresses at failure for IMS977-2 and T800S/M21 are 10% higher than these ones. The bolt is modelled by a steel disk (isotropic linear elastic material). A drawing of the model is presented Fig. 7. The prescribed force F is equal to the experimental maximum force applied during the test thus it varies for each configuration. The center of the bolt is fixed and node to node contact elements are used for modeling the contact between the disk and the composite plate. A nominal clearance is also introduced between the bolt and the hole and soft springs are used to avoid rigid body modes. Lastly, composite membrane elements with two degrees of freedom (displacement) are used to model the plate and to get the strain and stress fields in each ply of the laminate.

Results

All the configurations in the experimental database were computed with the 2D FEM. The Figure 8 shows the fiber strains for three of the four plies orientation (0°, 45° and 90°). These fiber strains are in compression (blue zones on Fig. 8) where the fiber directions are normal to the hole edge. They are in tension (red zones on Fig. 8) where the fiber directions are tangent to the hole edge. Assuming a first-play failure, the Hashin's criteria could be set in the differents plies according to the expected failure modes. Hashin's criteria [START_REF] Hashin | Failure criteria for unidirectional fiber composites[END_REF] is presented in Table 2. In this analysis, we only consider the fiber failure criteria.

Table 2. Hashin criteria

Hashin criteria Tension Compression

Fibers

Matrix

The values of the stress at failure are reported in Table 1. For bearing failure analysis, the Hashin criterion for fiber failure in compression is used for the 0° and ±45° plies. For cleavage failure mode, the tension fiber criterion in 90° plies is used. Finally, for net section failure mode, the tension fiber criterion is analysed in the 0° plies. A maxima stress criterion in fiber direction is also computed to study the influence of shear in the Hashin's criteria. Its expression is presented in Eq. 1 (T:Tension and C : Compression).

𝜎 11 𝑅𝑇/𝐶 𝑋 𝑇/𝐶 ≤ 1 (1)

Cleavage

The cleavage appears when the end distance is strictly smaller than 2,5. The Figure 9 shows the results obtained for the maximal stress and Hashin's criteria for tension fiber failure in 90° plies. Both criteria maximum values are smaller than 1 thus they do not predict the cleavage failure of the specimens. When the edge distance P is large, it is obvious since the experimentally observed failure is bearing. Nevertheless, this result is the same when P reduces. Moreover, the maximum criterion locations are different between maximum and Hashin's criteria: the Hashin's criterion is maximum where shear stress is maximum, at 40° from the load orientation, when the maximum stress criterion is maximum at 0° as illustrated in Fig. 9. Consequently, Hashin's criterion may over-estimate matrix failure influence on fiber failure. Eleven configurations of T700GC/M21 failed in cleavage. The maximum stress criterion gives an average of half the Hashin's criterion average. Their CoV are of 21% for the maximum stress criterion and of 18% for the Hashin's criterion while the experimental data CoV is 8.8% assuming a linear evolution of the failure bearing stress as a function of P. Both of the criteria give worse dispersions than the experimental one.

(

Net section

Only three specimens of T700GC/M21 have failed in net section in the database. Both maximum stress and Hashin's fiber failure criteria were computed in the 0° plies. They respectively give 36% and 28% CoVs while the CoV of the experimental data equals 14% assuming a constant net stress at failure. Finally, these two criteria gives worse dispersions than the experimental data.

Bearing

Bearing is a complex failure mode including several phenomena such as fiber failure in compression, delamination or matrix damage. The assumption here is that the bearing failure of the specimens occurs when 0° or ±45° plies fail due to fiber compression. The Hashin's criterion for fiber failure in compression is the same as the maximal stress criterion. But, it was shown in literature [START_REF] Budiansky | Compressive failure of fiber composites[END_REF] that shear can significantly reduce compression strength in fiber direction. Thus, a shear dependency was introduced in the Hashin's fiber criterion in compression as it is for tension. The maximal stress criterion gives a CoV of 10% for T700GC/M21 specimens. Let's note that failure mainly occurs in ±45° plies where maximum fiber compression is higher than in 0° plies. This result agrees with Gohorianu's observations [START_REF] Gohorianu | Interactions entre les défauts d'usinage et la tenue en matage d'assemblages boulonnés en carbone/epoxy[END_REF]. The Hashin's criterion gives a higher average value than the maximum stress criterion and a 7.4% CoV for the same specimens. The first-ply failure is balanced between 0° and ±45° plies depending on the geometries and lays-up. The CoV of the experimental data assuming a constant failure bearing stress is 6.3% so those criteria have worse dispersions than the experimental data. For IMS 977-2, the maximal stress criterion gives a CoV of 5.9% and the Hashin's criterion 7.7% with still the highest average value. For this material the CoV of both criteria are better than the CoV of the experimental data (12%).

Conclusions

The experimental database provided by Dassault Aviation was analysed to identify the influence of the design parameters on the failure of composite bolted joints. We observed that the end distance has a strong effect on joint failure mode and strength. The width and the loading direction are second order parameters for joint failure modes and strength.A 2D finite element modelling was then performed for improving the comprehension of the joint failure modes and strength. The different failure modes were predicted by a first-ply failure analysis using maximum stress and Hashin's criteria. We observed that the dispersions of the FE simulations were higher than the exprimental data assuming simple macroscopic criteria based on the failure bearing stress. Nevertheless, these simulations emphasized the significant effect of shear on the failure of composite bolted joints. In conclusion, the FE modeling must be improved and failure analysis of composite also. This should be done modeling the secondary bending effect that occurs in supported single-shear specimens (even if small) or introducing a non linear behavior of the composite material. The objective is to improve the simulation of strain and stress field distributions in supported single-shear specimens. The secondary bending illustrates the geometrical non linearities in joint behavior. It introduces out of plane stress in the laminate and thus modifies the stress field within the plate thickness. The non linear behavior of the composite material is related to the damage of its constituents such as fiber kinking, matrix cracks or delaminations. Our objective is to identify the role of each damage on the different failure modes and stresses.

Acknoledgement

The authors gratefully acknowledge the French DGAC for financially supporting this study

Figure 2 .

 2 Figure 2. Failure modes of composite bolted joints

Figure 3 .

 3 Figure 3. Bearing stresses at failure with respect to the edge distance

  illustrate these different failure modes showing the end edge of the failed specimens (L1 and L2 cleavage and L3 shearout) when P=1.5φ.The L1 and L2 specimens present a straight failure behind the bolt which is specific to cleavage. The 90° plies of the L3 specimen are strongly delaminated thus they are not within the laminates anymore.

Figure 4 .

 4 Figure 4. Post-failure specimens that failed in cleavage and shear-out (from left to right: L1, L2 and L3)

Figure 5 .

 5 Figure 5. Bearing stress at failure with respect to D (W=2D)

  th European Conference on Composite Materials Athens, Greece, 24-28 th June 2018 5 B. Montagne, F. Lachaud, E. Paroissien and D. Martini

Figure 7 .

 7 Figure 7. Boundary conditions of the 2D FE model

  ECCM18 -18 th European Conference on Composite Materials Athens, Greece, 24-28 th June 2018 6 B. Montagne, F. Lachaud, E. Paroissien and D. Martini

Figure 8 .

 8 Figure 8. ϵ 𝐹𝐹 in the different plies orientation (here P=3φ and D=2.5φ )

Figure 9 .

 9 Figure 9. Maximal stress criterion (A) with short and large end distance and Hashin's criterion (B)with short and large end distance for tension fiber failure in 90° plies

Table 1 .

 1 Mechanical properties of the T700GC/M21 UD ply.

	E 11	E 22 = E 33	G 12 = G 13	G 23	ν 12 = ν 13	ν 23 σ 11 RT	σ 11 RC	σ 22 σ 33 RT = RT	σ 22 RC = σ 33 RC	σ 12 R = σ 13 R	σ 23 R
	135000 8500 4800 3200 0.31 0.42 2300 1450	70	300	95	60


Valérie Pommier-Budinger

Marc Budinger
email: valerie.budinger@isae.fr

Pierrick Rouset

Fabien Dezitter

Florent Huet

Marc Wetterwald

Elmar Bonaccurso

Electro-mechanical Resonant Ice Protection Systems: Initiation of Fractures with Piezoelectric Actuators

Recent research is showing growing interest in low-power electromechanical de-icing systems and, in particular, de-icing systems based on piezoelectric actuators. These systems use the vibrations generated by piezoelectric actuators at resonance frequencies to produce shear stress at the interface between the ice and the support or to produce tensile stress in the ice. Many configurations of de-icing systems using piezoelectric actuators have been tested and showed that piezoelectric actuation may be a viable ice removal system. If the many experimental studies already achieved have the advantage to present tests in different configurations, they often lack analysis of the phenomena, which limits the optimization opportunities. This paper proposes a computational method for estimating voltages and currents of a piezoelectric de-icing system to initiate cohesive fractures in the ice or adhesive fractures at the ice/support interface. The computational method is validated by comparing numerical results with experimental results. Other contributions of this paper are the study of the types of mode (extensional or flexural) and of the frequency range with respect to de-icing performances and the proposal of some general rules for designing such systems while limiting their electric power consumption. 𝜎 𝑠ℎ𝑒𝑎𝑟 = transverse shear strength, MPa 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑖𝑐𝑒 = ice adhesive shear strength, MPa 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑋% = minimum shear stress among X% of the nodes with the highest shear stresses 𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒-𝑖𝑐𝑒 = ice tensile strength, MPa 𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒-𝑋% = minimum tensile stress among X% of the nodes with the highest tensile stresses 𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒-𝑖𝑐𝑒 = ice tensile strength, MPa 𝜎 𝑥𝑧 = transverse shear strength component on the face x and in the direction z, MPa 𝜎 𝑦𝑧 = transverse shear strength component on the face y and in the direction z, MPa  i = pulsation of mode i, rad/s

I. Introduction

CING occurs when an aircraft flies through clouds in which super cooled droplets are suspended in atmosphere with an ambient air temperature below the freezing point. The droplets impinge on the aircraft surfaces and freeze leading to ice accretion. The resulting change in the aircraft geometry can modify wing aerodynamic characteristics (loss of lift, rise of drag) or even damage the engine by ice ingestion. Regarding electrical de-icing systems, electro-thermal and electro-impulse technologies are already implemented on aircraft but are power-consuming. Thermal solutions require approximately from 2 up to 25 kW/m² for anti-icing and from 5 kW/m² up to 50 kW/m² for de-icing 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] Meier | A handbook method for the estimation of power requirements for electrical de-icing systems[END_REF]. Studies are currently ongoing to propose new solutions that consume less energy or have less bulky power supplies. Among these studies, the de-icing systems based on piezoelectric systems are a subject of growing interest. Note that one advantage of piezoelectric material is that it can also be used as sensors to determine ice characteristics [START_REF] Liu | Ultrasonic-attenuation-based technique for ice characterization pertinent to aircraft icing phenomena[END_REF]. The principle of de-icing systems based on piezoelectric actuators is to apply vibrations to the structure that create high-level stresses greater than those required to lead to delamination or cracking of the ice accumulated on the structure.

Delamination of the ice at the ice/support interface is linked to adhesive shear strength of the ice with the support and cracking is linked to the cohesive tensile strength of the ice.

I.I

The analysis of literature on de-icing systems using piezoelectric actuators shows that these vibrations are generated in different frequency ranges and by different types 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 piezoexcited 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][6] [START_REF] Venna | Mechatronic development of self-actuating in-flight de-icing structures[END_REF][8] 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 was more probably caused by thermal effects than by shear stress. S. Struggl et al. [START_REF] Struggl | A basic approach for wing leading de-icing 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 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 | Investigation of an ultrasonic ice protection system for helicopter rotor blades[END_REF] 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] Tarquini | Investigation of ice shedding properties of superhydrophobic coatings on helicopter blades[END_REF] he investigated the effect of hydrophobic coating combined with an ultrasonic de-icing system and showed that the ice adhesion depended on the support roughness. In [START_REF] Palacios | Ultrasonic de-icing of wind-tunnel impact icing[END_REF], he started using Finite Element models to predict the ultrasonic ice shedding transverse shear stresses responsible for ice shedding. He pursued this approach in [START_REF] Overmeyer | Ultrasonic de-icing bondline design and rotor ice testing[END_REF] to design a piezoelectric de-icing system able to promote shedding of ice layers ranging from 1.4 to 7.1 mm in thickness for varying icing conditions and in [START_REF] Diplacido | Enhancement of Ultrasonic De-Icing via Tone Burst Excitation[END_REF] to study the effect of tone burst excitation: enhanced de-icing ability is observed when employing multi-frequency bursts. These studies are interesting but were carried out without considering the tensile effects in the ice or the estimation of the electrical power consumption. More recently, Strobl [START_REF] Strobl | Feasibility study of a hybrid Ice protection system[END_REF] used multilayer piezoelectric patches at frequencies around 4 kHz and an icephobic coating to delaminate ice instantaneously on polished surfaces with a low supply voltage.

Villeneuve [START_REF] Villeneuve | Piezoelectric De-icing System for Rotorcraft[END_REF] used piezocomposite DuraAct to experiment de-icing on a leading edge sample. Using Abaqus software and introducing damping measurements in a harmonic analysis, a Finite Element model allowed him to compute the displacements for different modes but no link was made with the voltage or currents required to deice.

One conclusion of this review on de-icing systems using piezoelectric actuators is that different technologies of actuators and different frequency ranges have been tested and showed promising results. Another conclusion of this review is that many validations of the proposed architectures are experimental, the models proposed until now concerned mainly mechanical issues and not electrical issues. There is thus a lack of theoretical or numerical multiphysics analysis for assessing and designing electromechanical piezoelectric de-icing systems. In Budinger [START_REF] Budinger | Ultrasonic Ice Protection Systems: Analytical and Numerical Models for Architecture Tradeoff[END_REF], we compared different architectures of de-icing systems based on piezoelectric actuators and we focused on the frequency ranges and types of modes that are the most favorable for de-icing. We highlighted the fact that flexural modes and ultrasonic frequencies over 20 kHz allowed the shear stress at the ice/support interface to be maximized for a given displacement while minimizing fatigue in the structure and limiting breakdown of the actuators. We also promoted the use of pre-stressed piezoelectric actuators to avoid failure of the actuators.

The main contribution of the present paper is to present a numerical methodology to assess architectures of piezoelectric de-icing systems with respect to both mechanical stresses (in tensile and in shear mode) and electrical power consumption. It relies on a modal analysis that allows computation of voltages and currents (and thus power) for delaminating or cracking the ice on a structure at different resonance frequencies. This methodology is validated by comparing numerical and experimental results. The first experimental tests are carried out on a small rectangular plate. This very simple structure allows making comparison of numerical and experimental easier. Additional tests are performed with a leading edge sample to validate the method in a more realistic environment.

Another contribution of this article is to assess two kinds of frequency modes for ice delamination and cracking: extensional and flexural modes. To meet this objective, each test is performed on a narrow frequency range in order to excite only one mode and not several. For each mode, voltages and currents for delaminating or cracking the ice are measured in order to lay conclusions on the performance of each type of mode and help improve the design of piezoelectric de-icing systems.

Finally, based on outcomes from numerical analyses and experimental results, some design rules are proposed and some issues concerning piezoelectric de-icing systems are discussed.

The paper is organized as follows. Section 2 presents and justifies the choice of the test specimens used for the experiments. In section 3, the mechanical ice properties are selected according to the icing test conditions.

The proposed numerical methodology for assessment of piezoelectric de-icing systems method is detailed in section 4 and some first comments on the required inputs for the computation method are done. Section 5 shows the validation of the method on the simplest test specimen (plate-type specimen) for which numerical and experimental results are compared. The method is then assessed on a leading edge sample in section 6. Finally some rules for the design of piezoelectric de-icing systems are proposed in section 7.

II. De-icing Test Specimens

Our main objective in this paper is to propose and validate a numerical methodology for the preliminary design of piezoelectric de-icing systems. For the choice of the piezoelectric actuation system, the options currently available are pure ceramics, ceramic piezocomposites (made of pure ceramics encapsulated in a resin), Langevin transducers and piezoelectric stacks. Piezoelectric ceramics can be efficient to excite a large number of resonance frequencies on a large frequency range, in low and high frequencies. These ceramics exist in many shapes, can be cut to get the required size and are easy to implement. They must be bonded to the surface with a rigid glue to get a good coupling coefficient between the actuators and the structure. Piezocomposites have the advantage of being bendable and mechanically robust. Another advantage is that the piezoelectric ceramic used to manufacture the piezocomposite is encapsulated with prestressing in the resin and thus can withstand more stress than a pure piezoelectric ceramic [21][22]. These piezocomposites can also be efficient on a broad frequency range. They exist in a few sizes. The Langevin transducers are also made with pre-stressed ceramics and are thus resistant to stress. They can be screwed and thus disassembled, which is an advantage for the maintenance. Their drawback is that they can only be actuated at a single frequency. At last, piezoelectric stacks are piezoelectric linear actuators that provide high forces and high displacements compared to the other piezoelectric actuators, but that require special care for their mounting.

In our study, several resonance frequencies were used and their de-icing efficiency was compared. To carry on this study, we have used piezoelectric ceramics and piezocomposites for their ability to act on a broad frequency range and their ease of use. The objective of this article is not yet the optimization of the actuation system, thus fixed configurations were chosen to perform the analysis and the experiments. Two test specimens were used in this study.

Test specimen A consisted of an aluminum plate measuring 130x59x1.5 mm 3 with one piezoelectric ceramic (PIC 151size: 50x25x0.5mm 3 ) used as sensor and one piezocomposite DuraAct P-876.A15 from PI Ceramic as actuator. They were bonded at the clamped sides of the plate (Fig. 1). This simple structure allowed comparing simulation and experimental results while limiting uncertainties in the boundary conditions and in the ice profile. The test specimen A was tested in an icing wind tunnel. The ice was accreted at a temperature of -3

°C and with a flow velocity of 50 m/s. The conditions used for ice accretion were selected to obtain glaze ice. A liquid water content (LWC) of 0.3 g/m3 was measured and a droplet size (median volumetric diameter) of 22 µm is assumed at those settings. De-icing tests were then performed at a temperature of -10 °C and with a flow velocity of 25 m/s. A lower temperature and a lower speed were used in order to ensure that neither thermal nor aerodynamics effects would perturb the de-icing tests.

Fig. 1 -Test specimen A: plate with one piezocomposite and one piezoelectric ceramic

The test specimen B (Fig. 2) was closer to a more complex and realistic case. It consists in a leading edge sample in aluminum alloy that allowed testing the numerical method for a curved surface. The leading edge is a NACA 0012 profile of chord 350 mm and thickness 1mm. The actuation was performed by two piezocomposites DuraAct P-876.A15 bonded on one side of the leading edge. Another DuraAct P-876.A15 was used as sensor on the other side of the leading edge. This test specimen was also tested in the icing wind tunnel with the icing conditions above mentioned.

III. Mechanical Ice Properties

For a given architecture of piezoelectric de-icing system architecture, one of the required inputs is the type of ice and its mechanical properties. The work of Bennani [23] provides a good review of the studies carried out to determine the mechanical ice properties.

A. Young's modulus, Poisson's ratio and Density

Concerning the elastic behavior of ice, the Young's modulus, Poisson's ratio and density can be found in Gammon [START_REF] Gammon | Elastic constants of artificial and natural ice samples by Brillouin spectroscopy[END_REF]. Eskandarian [START_REF] Eskandarian | Ice shedding from overhead electrical lines by mechanical breaking[END_REF] has estimated the Young's modulus for atmospheric ice and Mohamed [START_REF] Mohamed | An experimental study on the tensile properties of atmospheric ice[END_REF] gave its density. All these results are given for glaze ice. It can also be interesting to have the elastic behavior for other types of ice and Gao [START_REF] Gao | Ice detection and classification on an aircraft wing with ultrasonic shear horizontal guided waves[END_REF] proposed some values for rime ice and mixed ice.

B. Tensile Strength and Adhesive Shear Strength

To be effective, piezoelectric de-icing systems must deform ice sufficiently to break and detach ice accumulations on surfaces. It implies two kinds of fractures: cracks in the ice layer or delamination of the ice at the ice/support interface (Fig. 3). The appearance of cracks is linked to the cohesive tensile strength of the ice and delamination to the adhesive shear strength of the ice with the support. Although many experiments have been carried out to measure these values, there is a relatively wide range of ice strengths and it is not easy to select a single value because ice strength depends on many factors: temperature [START_REF] Fortin | Ice adhesion models to predict shear stress at shedding[END_REF] [29] [START_REF] Petrovic | Review mechanical properties of ice and snow[END_REF][37], nature and roughness of the support [START_REF] Fortin | Ice adhesion models to predict shear stress at shedding[END_REF] [37], strain and strain rate [30][31], ice grain size [START_REF] Petrovic | Review mechanical properties of ice and snow[END_REF], flow speed [32][37]. Based The reduced model of a structure with k piezoelectric actuators with n degrees of freedom can be written as [START_REF] Pommier-Budinger | Sizing optimization of piezoelectric smart structures with meta-modeling techniques for dynamic applications[END_REF]:

{ [𝑀] 𝑛𝑥𝑛 {𝑞̈} 𝑛 + [𝐷 𝑠 ] 𝑛𝑥𝑛 {𝑞̇} 𝑛 + [𝐾] 𝑛𝑥𝑛 {𝑞} 𝑛 = [𝑁] 𝑛𝑥𝑘 {𝑉} 𝑘 {𝑞 𝑐 } 𝑘 = [𝑁] 𝑇 𝑘𝑥𝑛 {𝑞} 𝑛 + [𝐶 𝑜 ] 𝑘𝑥𝑘 {𝑉} 𝑘 (1)
This model consists in a mechanical equation and an electrical equation where q the vector of generalized coordinates, M the mass matrix (nxn), K the stiffness matrix (nxn), D s the damping matrix (nxn), N the electromechanical coupling factor matrix (nxk), q c the electrical charge vector (kx1), V the voltage vector (kx1), and C o the modal turned-off capacity matrix (kxk).

The parameters of the model can be computed analytically or using Finite Element software (such as COMSOL©, ANSYS© or ABAQUS©) that allows computations with piezoelectric elements. The method is described here for Finite Element software which enables handling easily complex geometry unlike analytical calculations.

The natural frequencies associated to the eigenmodes of this system are those that verify the relation:

| [𝐾] -𝜔 𝑖 2 [𝑀] | = 0. (2)
This is an eigenvalues problem and for each natural frequency there will be an associated eigenvector which verifies [𝐾]{Φ} 𝑖 = 𝜔 𝑖 2 [𝑀]{Φ } 𝑖 . These eigenvectors represent the mode shapes of each eigenmode. The eigenvectors can be normalized to the mass matrix and, thanks to the orthogonality of the eigenmodes, it is possible to find a base of eigenvectors that verifies the following properties:

{Φ} 𝑖 𝑇 [𝑀]{Φ } 𝑖 = 1 (3)
and

{Φ} 𝑖 𝑇 [𝐾]{Φ } 𝑖 = 𝜔 𝑖 2 (4)
This way we can obtain the modal matrix:

[Φ] = [{Φ} 1 , {Φ} 2 , … , {Φ} 𝑛 ]. (5)
This matrix represents a change of basis between the initial generalized coordinates {𝑞} and the coordinates {𝑢} that are the modal coordinates. This means:

{𝑞} = [Φ]{𝑢}. (6)
Let us consider the basis in which the mass matrix becomes the identity and the stiffness matrix is diagonal and composed by the squared natural frequencies.

[𝑀 * ] = [𝐼] (7)
[𝐾 * ] = 𝑑𝑖𝑎𝑔[𝜔 𝑖 2 ] (8)
[Φ] is chosen such that:

[Φ] = [𝑢 𝑖𝑗 ] (9)
where 𝑢 𝑖𝑗 is the modal displacement of the point 𝑖 (out-of-plane for flexural modes and in plane for extensional modes) for the j th natural frequency.

For our method, the analysis will be carried out frequency by frequency, for one degree of freedom and one natural frequency at each time. Consequently there is only one value for 𝑢 and the modal mass and modal stiffness for the considered mode are written as:

𝑀 * = 𝑢𝑀𝑢 = 𝑢 2 𝑀 = 1 (10)
𝐾 * = 𝑢𝐾𝑢 = 𝑢 2 𝐾 = 𝜔 2 , (11)
and thus the mass and stiffness of equation ( 1) are obtained by:

𝑀 = 1 𝑢 2 (12) 𝐾 = 𝜔 2 𝑢 2 . ( 13
)
In order to obtain the force factor, the modal analysis used for the computation of 𝑀 and 𝐾 is run with all the ceramic surfaces in short circuit, that is {𝑉} 𝑘 = 0 𝑘 applied on all the piezoelectric ceramic surfaces. This way the force factor is easily estimated by using the electrical equation of model [START_REF] Meier | A handbook method for the estimation of power requirements for electrical de-icing systems[END_REF] with {𝑉} 𝑘 = 0 𝑘 thus giving:

𝑁 = 𝑞 𝑐 𝑢 . ( 14
)
The last missing term in the reduced model ( 1) is the modal damping matrix. According to Hasselman [START_REF] Guerin | Analytical model based on experimental data of centrifuge ice adhesion tests with different substrates[END_REF], it is formulated for one mode as:

𝐷 𝑠 = 2𝜀√𝐾𝑀 ( 15
)
where 𝜀 is the damping ratio of the considered mode. The damping ratio is generally difficult to calculate because it depends on too many experimental parameters and it is preferable to measure it. It can also be defined according to the quality factor:

𝑄 𝑚 = 1 2𝜀 . (16)
Note: in our experiments, the damping ratios are low and the natural frequency and resonance frequency of each mode can be considered as identical. The next sections of the article refer to the resonance frequency which is the value that is experimentally measured.

B. Estimation of the Required Voltage for Ice Delamination or Cracking

Estimation of the Voltage, Current and Power for Ice Delamination Based on Adhesive Shear Strength

Ice delamination using piezoelectric transducers occurs when the transverse shear stress at the ice/support interface exceeds the adhesive shear strength of the ice accumulated onto the surface. The objective is to detach the ice accretion by the failure of the adhered interface. In this section, shear stresses at the ice/support interface are used to estimate the voltage, current and power required by the transducer to achieve ice delamination. With respect to the frame of Fig. 1, the transverse shear stresses are 𝜎 𝑥𝑧 and 𝜎 𝑦𝑧 and the transverse shear stress 𝜎 𝑠ℎ𝑒𝑎𝑟 is calculated as the resultant of 𝜎 𝑥𝑧 and 𝜎 𝑦𝑧 stress components:

𝜎 𝑠ℎ𝑒𝑎𝑟 = √(𝜎 𝑥𝑧 2 + 𝜎 𝑦𝑧 2 ) (17)
This shear stress 𝜎 𝑠ℎ𝑒𝑎𝑟 is computed through a modal analysis for the different resonance frequencies of the structure in a given frequency range.

It is possible to compute the maximal value of the shear stress, but it will correspond to the delamination of a single point of the ice/support interface. The objective is to achieve delamination in several points of the surface.

To reach this objective, 𝑋% of the nodes for which the shear stress is maximum are selected and the minimum shear stress of these nodes is considered. The nodes are selected on the ice surface in contact with the support.

We define 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑋% .as:

𝜎 𝑠ℎ𝑒𝑎𝑟-𝑋% ≡ minimum shear stress among X% of the nodes with the highest shear stresses

Using this X-percentile shear stress 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑋% computed by modal analysis and with the maximal modal displacement 𝑢 𝑚𝑎𝑥 at the ice/support interface also computed by modal analysis, it is possible to compute the required displacement, called 𝑞 𝑑𝑒𝑙 , that has to be reached on the node of 𝑢 𝑚𝑎𝑥 in order to overcome the ice adhesive shear strength 𝜎 𝑠ℎ𝑒𝑎𝑟_𝑖𝑐𝑒 . Note that 𝑞 𝑑𝑒𝑙 is a displacement calculated in the real basis (not the modal)

and is of the order of µm. 𝑞 𝑑𝑒𝑙 is expressed as:

𝑞 𝑑𝑒𝑙 = 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑖𝑐𝑒 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑋% 𝑢 𝑚𝑎𝑥 (19)
The estimation of the voltage, current and power requires the computation of the mass, stiffness and force factor of each mode. This is done by using the method previously described in section A and considering the maximum modal displacement 𝑢 𝑚𝑎𝑥 :

𝑁 = 𝑞 𝑐 𝑢 𝑚𝑎𝑥 ( 20
)
𝐾 = 𝜔 2 𝑢 𝑚𝑎𝑥 2 (21) 𝑀 = 1 𝑢 𝑚𝑎𝑥 2 (22)
Then as the stationary response of the forced mechanical equation of the reduced model (1) expressed for one resonance mode 𝜔 𝑖 leads to:

𝑗𝐷 𝑠 𝜔 𝑖 𝑞 = 𝑁𝑉 (23)
the value of the voltage 𝑉 𝑑𝑒𝑙 that generates the displacement 𝑞 𝑑𝑒𝑙 and thus leads to the delamination of Xpercentile of the nodes is then written:

𝑉 𝑑𝑒𝑙 = 𝐷 𝑠 𝑁 𝜔 𝑖 |𝑞 𝑑𝑒𝑙 | (24)
Using equations ( 15), ( 19), ( 20) and ( 21), |𝑉 𝑑𝑒𝑙 | can also be expressed as:

𝑉 𝑑𝑒𝑙 = 𝐾 𝑁𝑄 𝑚 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑖𝑐𝑒 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑋% |𝑢 𝑚𝑎𝑥 | (25)
For the estimation of the current, the electric equation of the reduced model ( 1) is now considered. The electric current is defined by 𝐼 = 𝑑𝑞 𝑐 𝑑𝑡 and is thus equal to:

𝐼 = 𝑁𝑞̇+ 𝐶 0 𝑉 ̇ ( 26
)
Current I is shared between motional current 𝐼 𝑚 and capacitive current 𝐼 𝑐 due to the dielectric behavior of the piezoelectric material:

𝐼 𝑚 = 𝑁𝑞̇= 𝑗𝜔 𝑖 𝑁𝑞 (27)
𝐼 𝑐 = 𝐶 0 𝑉 ̇= 𝑗𝜔 𝑖 𝐶 0 𝑉

Like for the voltages, the currents are estimated for the delamination of X-percentile of the nodes:

𝐼 𝑚-𝑑𝑒𝑙 = 𝑗 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑖𝑐𝑒 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑋% 𝜔 𝑖 𝑞 𝑐 (29)
𝐼 𝑐-𝑑𝑒𝑙 = 𝑗𝜔 𝑖 𝐶 0 𝑉 𝑑𝑒𝑙 = -𝜔² 𝑖 𝐶 0 𝐷 𝑠 𝑁 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑖𝑐𝑒 𝜎 𝑠ℎ𝑒𝑎𝑟-𝑋% 𝑢 𝑚𝑎𝑥 (30)

𝐼 𝑑𝑒𝑙 = √ |𝐼 𝑚-𝑑𝑒𝑙 | 2 + |𝐼 𝑐-𝑑𝑒𝑙 | 2 (31)
Once the required voltage and currents are known, it is simple to obtain the estimated apparent power 𝑆 𝑑𝑒𝑙 and active power 𝑃 𝑑𝑒𝑙 for the delamination of X-percentile of the nodes:

𝑆 𝑑𝑒𝑙 = 𝐼 𝑑𝑒𝑙 . 𝑉 𝑑𝑒𝑙 ( 32
)
𝑃 𝑑𝑒𝑙 = |𝐼𝑚_𝑑𝑒𝑙 |. 𝑉 𝑑𝑒𝑙 ( 33
)

Estimation of the Voltage, Current and Power for Ice Cracking based on Tensile Strength

Ice cracking using piezoelectric transducers occurs when the tensile stress in the ice exceeds the tensile strength of the ice accumulated onto the surface.

As for the shear stress, the tensile stress is computed for the different resonance frequencies thanks to a modal analysis. It is also computed to crack in several points of the surface.𝑋% of the nodes for which the tensile stress is maximum are selected and the minimum tensile stress of these nodes is considered. The nodes are selected on the upper ice surface (opposite to the ice/support interface). 𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒-𝑋% is defined as :

𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒-𝑋% ≡ minimum tensile stress among X% of the nodes with the highest tensile stresses [START_REF] Guerin | Analytical model based on experimental data of centrifuge ice adhesion tests with different substrates[END_REF] Using this X-percentile tensile stress 𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒-𝑋% and the maximal displacement 𝑢 𝑚𝑎𝑥 over the structure, it is possible to estimate the required displacement, called 𝑞 𝑐𝑟𝑎𝑐𝑘 that has to be reached on the node of 𝑢 𝑚𝑎𝑥 in order to overcome the ice tensile strength 𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒_𝑖𝑐𝑒 :

𝑞 𝑐𝑟𝑎𝑐𝑘 = 𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒-𝑖𝑐𝑒 𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒-𝑋% 𝑢 𝑚𝑎𝑥 (35)
The mass, stiffness and force factor of each mode are computed as previously, as well as the value of the voltage 𝑉 𝑐𝑟𝑎𝑐𝑘 that generates the displacement 𝑞 𝑐𝑟𝑎𝑐𝑘 and thus leads to the cracking of X-percentile of the nodes:

𝑉 𝑐𝑟𝑎𝑐𝑘 = 𝐷 𝑠 𝑁 𝜔 𝑖 |𝑞 𝑐𝑟𝑎𝑐𝑘 | = 𝐾 𝑁𝑄 𝑚 𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒-𝑖𝑐𝑒 𝜎 𝑡𝑒𝑛𝑠𝑖𝑙𝑒-𝑋% |𝑢 𝑚𝑎𝑥 | (36)
The formulas of the currents estimated for cracking of X-percentile of the nodes as well as the formulas of the powers are established in the same way as in the previous section.

C. Comments on the Numerical Methodology for Assessment of Piezoelectric De-icing Systems

Estimated Voltage, Current and Power for ice shedding

The estimated voltages, currents and powers depend strongly on two values which are subject to a large scatter:

 the adhesive shear strength of the ice σ shear-ice or the tensile strength of the ice σ tensile-ice  the quality factor Q m of the structure for the considered mode (note: the motional current is the sole value independent of the quality factor).

The voltages and current depend linearly of the shear or tensile strength of the ice and the power is proportional to the square of these values.

Quality Factor

Regarding the quality factor, it is has already been highlighted that this coefficient is difficult to estimate and that it will be measured. Measurement results will be given in the next section and will show that this coefficient depends on the displacement amplitude of the structure (and thus on the voltage amplitude applied to the actuator) and on the kind of modes (flexural or extensional).

Measured Resonance Modes and Classification of the Resonance Modes Given by the Method

In theory, on a given frequency range, all the resonance modes are computed. In practice, experimental results will show that the number of modes that can be observed is smaller than the number of modes that are computed. The modes that are the most likely to be measured are the modes for which the piezoelectric transducers are the best coupled with the structure, thus the modes for which the electromechanical coupling factors are the highest for both the piezoelectric actuators and sensors. Of course, these coupling factors depend on the number, size and location of the piezoelectric transducers. The optimization phase of the actuation system consists in providing a configuration with high coupling factors that can excite efficient modes for ice shedding i.e. that can completely remove ice with minimal consumed power and energy. This is a very important issue that is not addressed in this paper even if some recommendations will be given at the end of the paper for the design of piezoelectric actuation system.

Note: in order to read in a relevant way the results given by the numerical method, the modes computed by the Finite Element Analysis are classified by decreasing electromechanical coupling factors of the actuators. In this way, the first modes that appear in the results are the most likely to lead to ice shedding. Of course, these modes depend on the way the actuators are coupled with the structure. An optimized coupling would give the optimized modes.

Other Criteria for the Design of Robust Piezoelectric De-icing Systems

Additionally to the previous classification, two other filters are applied to the results in order to eliminate the modes that lead to solutions that could deteriorate the piezoelectric de-icing systems by non-respect of physical or mechanical properties of the material.

The first filter concerns the maximal voltage that can be applied to the actuators. For piezoelectric ceramics, this voltage is linked to the maximal electric field that can be applied without depolarizing the ceramics. Thus it is necessary to consider only cases where the voltages do not lead to the maximal electric field acceptable by the piezoelectric ceramics (2000V/mm for unipolar static and low-frequency applications and 400V/mm for highfrequency bipolar applications). For piezocomposites, the datasheet of DuraAct P-876.A15 from PI gives operating voltages from -250V to 1000V but without indicating if these values are for low-frequency or highfrequency applications. To avoid depolarization risks, we will consider that these values are only valid for lowfrequency applications.

The second filter concerns the mechanical rupture: the stresses in the piezoelectric actuators must be below the maximal stresses allowable for the piezoelectric ceramics or piezocomposites before ice fracture occurs. For piezoelectric ceramics, we consider that there is a risk of mechanical rupture if the stress in the ceramics exceeds 25 MPa. For piezocomposites, as the ceramics are pre-stressed, we consider a limit of 50 MPa [START_REF] Gall | Finite Element Formulation of a Piezoelectric Continuum and Performance Studies of Laminar PZT-Patch-Modules[END_REF]. The maximum stress in the bond must also be taken into account in order to avoid debonding of the actuators.

Note: other criteria can be taken into account to the methodology provided that values are available. It is possible for example to consider the cycle life of actuators submitted to cyclic loading.

V. Validation of the Numerical Methodology through Experiments on a Plate Sample

A. Test Procedure

In this article, our objective was not to deice a structure completely, but to validate the numerical methodology to estimate the voltages and currents required to initiate the delamination or the cracking of ice at one given frequency. Thus, for each test, the following values were measured:

 the quality factor to enter this value in the model for the analysis of the test  the current and voltage to compare the estimated and real values

The test procedure was established in order to obtain these values. One test corresponded to measurements in a frequency range in the vicinity of the frequency f r that needed to be excited. This frequency was measured on the structure with the ice layer. For each test, the voltage applied to the actuator was set as following:

 The frequency was swept on a narrow bandwidth, [f r -x Hz, f r + x Hz], x being selected for each test such that only one resonance mode is excited The frequency was swept in the frequency range with a step of 1Hz

for the tests in low frequencies and 5Hz for the tests in high frequencies. Each step lasted 0.1s approximatively. This value ensured that the steady state is reached.

 the amplitude was increased by step of 10V starting from 10V for the plates and 50V for the leading edge till the first crack or delamination was observed.

For each voltage applied to the actuator(s), the current consumed by the latter(s) was measured by a current probe. One piezoelectric transducer bonded on the structure was used as sensor to evaluate the quality factor by measuring the voltage at its terminals.

B. Comparison of Experimental and Numerical Results

Different tests at different frequencies were performed for the test specimen A. For each frequency and its measured quality factor, the numerical results gave the voltage and current expected to have the first delamination (computed based on the adhesive shear strength and a percentage of nodes of 1, 2 or 5%) and the first crack (computed based on the cohesive tensile strength and a percentage of nodes of 1, 2 or 5%). These theoretical values were then compared to the measured values. To obtain valid theoretical values, great attention was paid to the mesh of the structures. It was all the more true for values computed for the shear stress criterion.

A mesh convergence study allowed concluding that the mesh must be performed with 8 elements in the ice thickness and elements of 1 mm of side. The number of elements in the plate thickness is less critical since the stresses in the plate are not considered in our computations.

Results for a Flexural Mode around 5 kHz

A first test has been performed with the test specimen A tested around 5 kHz, which corresponded to a flexural mode (Fig. 4). The ice was formed as described in section II. The ice layer did not cover the whole plate since the plate was longer than the length of the wind tunnel test section (the length of the plate was the length of the wind tunnel test section plus the length of the two flanges). With respect to the pictures of the ice profile and measurements done after cracking, the ice layer was assumed to be 100 mm long and 50 mm width and its thickness was assumed to be constant and equal to 5 mm. These values were used for the simulations. The quality factor was measured during the test and is plotted in Fig. 5. The quality factor is low and tends to decrease with increasing voltages. This poor quality factor is justified by the loss in the metallic support and also in the Plexiglas flanges. Numerical computations (Table 1) were run with a quality factor of 25. Table 1 gives the estimated voltage and current to have a certain percentage of nodes cracked or delaminated (according respectively to cohesive tensile criteria and adhesive shear criteria) to see the influence of this coefficient in the results. This table also gives the computed stress in the piezoelectric actuators to estimate the solicitations in the actuator and the risks of damage. For information, the maximal displacement on the plate to obtain ice cracking or delamination is also reported in the table. The first fractures corresponded to cracks that appeared at 90V for a consumed current of 0.11A (Fig. 6).

Even if the ice profile is not perfectly modeled, numerical results are good to predict the voltage initiating the first cracks. Indeed the predicted voltage according to cohesive tensile criterion is 85 V with 2% of nodes used for the computation and the predicted current is 0.06 A. Thus, the prediction of the voltage is good. The prediction of the current is not as accurate but, as the current magnitude is small in this test, there may be uncertainties in the measurements.

Results for an Extensional Mode around 38 kHz

Extensional modes are present at higher frequencies than flexural modes. A test was performed around 38 kHz that corresponds to the extensional mode described in Fig. 7(a). Compared to the flexural mode at 5 kHz tested on the same test specimen, the quality factor of this extensional mode was slightly larger. This difference on the quality factor has been observed for different test specimens: the quality factors of extensional modes tend to be higher than the quality factors of flexural modes. According to the measurements performed during the tests, a value of 45 was considered to obtain the numerical results.

For this high frequency (compared to the one tested previously), the numerical results predicted that the initiation of cracking and delamination occurred for a same voltage of 145V and a current of 0.89A (Table 2, results with 1% of nodes used for the computation). This result is consistent with the picture of the test showing that ice cracking and delamination appeared simultaneously for a voltage of 140V with a current of 0.9A (Fig.

7(b)

). Note that more ice was removed from the plate compared to tests with flexural modes, but the current consumed during the excitation of this extensional mode was much higher (among other reasons because the frequency used for excitation was higher).

C. Conclusion on the Validity of the Method

The numerical results of the methodology for the assessment of piezoelectric de-icing systems depend strongly on the quality factor and the values of the tensile and shear strengths. If the uncertainties on these parameters are small, the method gives a good estimation of the voltages and currents required to cause the first crack in the ice or the first delamination at the ice/support interface. This estimation can be obtained for simulations run with 1% of the nodes where stress is higher than the limit strength, in tensile or in shear. Indeed, the computation results for 1 or 2% are very close. In a general way, performing the computations with 1% gives a good approximation.

The uncertainties on tensile and shear strengths can be reduced with measurements. The main issue for obtaining valid predictive results is to have a value for the quality factor. It varies for each resonance mode and with the type of mode. It tends to be larger for extensional modes or for coupled modes where extension is important than for flexural modes. The quality factor depends on the boundary conditions: it is larger for free conditions than for clamped conditions. For clamped conditions, it depends on the material of the support (the quality factor of a material is equal to the inverse of the loss factor that varies considerably with the material; for e.g., the loss factor of aluminum alloy is around 10 -4 while the loss factor of Plexiglas is 10 -2 ) and on the clamp systems (glue, screw, rivets,…). The quality factor also varies nonlinearly with the temperature and the magnitude and frequency of the vibrations. In the temperature range, frequency range and displacement range of our experiments, it decreases as the amplitudes of the vibrations increase. However, despite the influence of all these parameters, the quality factor seems to reach a minimal value. Experiments on the plates have shown that 25 was a minimum when the plate was fixed onto the flanges in Plexiglas. For metallic fixtures and support in aluminum alloy, a minimal value of 30 seems to be attainable if flexural modes are excited. If extensional modes are excited, this value can be higher (40 to 50).

Another conclusion that can be drawn from these first experiments is that, for bare support (not polished nor covered by a coating), the de-icing is mainly due to tensile effects, especially at low frequencies. It confirms the results of the analytical studies performed in [START_REF] Budinger | Ultrasonic Ice Protection Systems: Analytical and Numerical Models for Architecture Tradeoff[END_REF] for flexural modes. In this article, it was shown that, for bending plates, the shear stress is higher than the tensile stress for frequencies above 40 kHz approximately.

VI. Numerical and Experimental Results on a Leading Edge Sample

To further test the methodology, the leading edge sample described in section 2 was tested. In the icing wind tunnel, the ice accretion was on the top of the leading edge, as shown in Fig. 9 and Fig. 11. The two actuators were just below the iced surface and could be actuated in two ways: symmetric modes where the two actuators were supplied with the same voltage, and antisymmetric modes where the two actuators were supplied with opposite voltage. Tests were mainly performed in anti-symmetric modes because the simulation, which was confirmed by the experiments, showed that these modes required, for a given quality factor, lower voltages than symmetric modes.

First, the numerical method based on the FEM analysis was run with a quality factor of 30 and values of cohesive tensile strength of 1.5 MPa and adhesive shear strength of 0.5MPa. An ice profile as close as possible to the real ice profile was selected for the simulation. The computation highlighted frequencies of 31 kHz, 12 kHz, 6.5 kHz and 3 kHz as potential frequencies for de-icing (based on the force factor criterion).

Then, a spectral analysis on the bandwidth [100-50k] Hz was performed in order to check which modes were well-coupled and to compare these modes with the modes emphasized by the numerical method. Of course, the modes that could be observed were not only the well-coupled modes but also the modes that could be measured by the sensors. Among the resonance modes measured on the iced structure, two modes were tested: the mode around 3 kHz and the mode around 31 kHz because they were both among the first ones given by the numerical results and they could be easily observed.

A. Results around 3 kHz for the Leading Edge Sample

The mode at 3 kHz is mainly a flexural mode (Fig. 8). For this frequency and a quality factor of 30, the numerical results predict a first crack at 73 V with a current of 0.11 A according to the tensile strength criterion applied to 1 % of nodes for the computation. For this voltage, the numerical results estimated displacements of 140 µm and stresses in the piezoelectric actuators of around 19.5 MPa. The tests confirmed the value of the quality factor and the first ice crack occurred at 90 V with a current of 0.17 A (Fig. 9 (a)), which is not too far from the numerical results. The difference can come from imperfections in the ice profile model within the simulation. At 150 V with a current of 0.59 A, more ice shedding was observed (Fig. 9 (b)). Note that the quality factor was measured after the first crack and that it had decreased to 20. Although this mode was a flexural mode, compared to the experiment on the plates, the ice was well removed from the surface, especially at the tip of the leading edge. Further investigations are required to understand if this high de-icing efficiency with a flexural mode comes from the curved profile, the effect of the aerodynamic forces on the profile, or another effect. The mode at 31 kHz is shown in Fig. 10. It corresponds to an extensional mode coupled with flexion. First numerical computations were run with a quality factor of 30 but had to be reevaluated with the measured quality factor of 45. The numerical results then predicted a first crack at 110 V with a current of 0.95 A according to the tensile strength criterion applied to 1% of nodes for the computation.

Systems

In this last section, based on the analysis of the numerical and experimental results, some rules are proposed for assisting in the design of piezoelectric de-icing systems.

Low-power de-icing systems are systems that consume little power and that are lightweight, which is favorable for minimizing the overall consumption of the aircraft. As shown in equations ( 27) and ( 28), the current magnitudes increase with the frequency. To study the mass and energy consumption of piezoelectric deicing systems, it would be necessary to consider the electronic power supply of such systems. Two main technologies can be used: linear amplifiers and switching inverters. Linear amplifiers can generate perfectly sinusoidal waveforms but have very low energy efficiency, in particular with capacitive loads such as piezoelectric actuators. The losses to be dissipated, and thus the weight and volume of heatsinks, increase with frequency. Switching inverters have better energy efficiency but have to integrate heavy inductive components in order to compensate capacitive reactive energy of piezoelectric loads. Moreover, the nature of the piezoelectric loads at the resonance can vary quickly from capacitive to inductive. This can also be an issue for this type of power supply.

Another point that can be considered to reduce power consumption is the increase of the quality factor of the structures to be de-iced. If possible, for clamped structures, it is better to avoid materials with high losses and to favor rivets over screws, for example. For bonded structures, the bond layer must be minimized because it is a source of losses.

At last, to optimize the consumption of piezoelectric de-icing systems, the ice thickness should be taken into account. Firstly, thick ice layers are detrimental from an aerodynamic point of view and for the overall consumption of the aircraft. Secondly, the shedding of thick ice layers is more demanding in terms of power. On Fig. 13, the required power to induce the first cracks is plotted versus the ice thickness. This figure is plotted using computed results obtained on the clamped plate and considering the tensile criterion. It corresponds to a flexural mode whose frequency varies from 3 kHz to 4.6 kHz with ice thickness. Regarding the criteria of power consumption, Fig. 13 highlights that it is preferable to limit the thickness of the ice to be removed, but with some limits to reach the optimal thickness. Moreover, it is well-known that ice thickness should be also limited for aerodynamics issues. For all these reasons, there is clearly an interest in activating the piezoelectric de-icing system to avoid thick ice layers. The optimal value of ice thickness cannot be given here since more work on the deicing mechanisms and more results on deicing systems efficiency are required to give a clear answer to this question. Finally, another criterion that must be taken into account for the choice of the frequency range is the noise generated by the active de-icing system. It is obvious but it must not be neglected. Below 20 kHz, the piezoelectric de-icing system produces audible frequencies and the level of sound, especially in the range [START_REF] Meier | A handbook method for the estimation of power requirements for electrical de-icing systems[END_REF][START_REF] Liu | Ultrasonic-attenuation-based technique for ice characterization pertinent to aircraft icing phenomena[END_REF][START_REF] Ramanathan | De-icing of helicopter blades using piezoelectric actuators[END_REF][START_REF] Kalkowski | Removing surface accretions with piezoexcited high-frequency structural waves[END_REF][START_REF] Venna | In-Flight De-Icing Self-Actuating Wing Structures with Piezoelectric Actuators[END_REF][START_REF] Venna | Development of Self-Actuating In-Flight De-Icing Structures with Power Consumption Considerations[END_REF][START_REF] Venna | Mechatronic development of self-actuating in-flight de-icing structures[END_REF] kHz where the human ear is the more sensitive, can be seriously disturbing. If this noise can propagate through the structures till the cabin, it can become an issue.

VIII. Conclusions

This article provides a computational method to estimate voltages and currents required to induce the first ice cracks or first ice delamination with piezoelectric de-icing systems. The method is based on a modal analysis of the structure where the needs to be removed and was validated through experiments. In addition to the benefit of the proposed methodology, these experiments also allow drawing conclusions on piezoelectric deicing systems. The different experiments presented in this paper were conducted with structures in bare aluminum alloys.

These experiments and the associated numerical results allow showing that, with such material, the initiation of de-icing is more likely due to tensile stress in the ice than shear stress at the interface between the ice and the support, especially at low frequencies. Initiation of de-icing due to shear stress occurs for high frequencies, above around 40 kHz for our study cases. If ice-phobic treatments were applied, adhesive shear strength could be reduced, as well as the voltage required to have ice delamination by shear stress. Initiation of delamination could then occur before cracking, even at low frequencies. It would change the de-icing mechanisms.

This article also tried to compare flexural and extensional modes. The flexural modes come in low and high frequencies. For low frequencies, the required power to initiate the cracks is less than for high frequencies.

However, on plates, whatever the frequencies, for flexural modes, only cracks occurred and ice was still stuck on the plate at the end of the tests performed till 200V. For curved surfaces such as the leading edge, the results are different since both ice cracking and delamination were observed for flexural modes. These results are not yet fully understood and should be further investigated. It could be an important phenomenon for the design. For extensional modes, the de-icing efficiency was better for the plate and the leading edge, but as these modes exist at rather high frequency, the power consumption is much higher than for flexural modes.

For piezoelectric de-icing systems based on vibrations generation, the damping or quality factor of the structures is an important issue. All the design choices that permit to limit losses and as such increase the quality factor are to be favored.

At last, other parameters are to be taken into account in the design of low-power piezoelectric de-icing systems such as the electronic power supplies and acoustic noise issues. To find the good mode(s) for de-icing, more understanding of de-icing mechanisms is required to obtain complete ice shedding and not only cracks or partial delamination. Optimization is certainly required to take into account all the design drivers and compromise may be necessary at the final stage of the design.

Finally, concerning the perspectives to improve the computational method, the method proposed here is, in our point of view, the first step to the design of piezoelectric de-icing systems. To go further than the prediction of the de-icing initiation, more work is required to model and understand what happens after the generation of initial cracks or delamination, and how the cracks or delamination can propagate.

  the delamination of X-percentile of the nodes of the structure, A I c-del = capacitive current for the delamination of X-percentile of the nodes of the structure, A I m-del = motional current for the delamination of X-percentile of the nodes of the structure, A for the delamination of X-percentile of the nodes of the structure, W q = generalized coordinate q c = electrical charge vector q crack = required displacement for initiating cracks, m q del = required displacement for initiating delamination, m Q m = quality factor, -𝑆 𝑑𝑒𝑙 = apparent power for the delamination of X-percentile of the nodes of the structure, VA 𝑢 𝑖𝑗 = displacement of the point 𝑖, for the j th natural frequency, m 𝑢 𝑚𝑎𝑥 = maximal displacement over the structure under study, m V = voltage vector V crack = voltage generating the displacement q crack , V V del = voltage generating the displacement q del , V  = damping ratio, - = modal coordinate, m  = transformation matrix

Fig. 2 -

 2 Fig. 2 -Test specimen B: (a) Leading edge sample with three piezocomposites (two as actuators on the upper side and one as sensor on the lower side). (b) Positions of the piezocomposites on the upper side

  several articles[START_REF] Makkonen | Ice adhesion-theory, measurements and countermeasures[END_REF][36][START_REF] Laforte | De-icing strains and stresses of iced substrates[END_REF][34][START_REF] Jellinek | Adhesive properties of ice[END_REF][START_REF] Schulson | The structure and mechanical behavior of ice[END_REF], the range of the ice tensile strength for fresh water is between [0.6- 3] MPa and between [1.3-1.66] MPa for atmospheric water although the range of the shear strength of the ice with the support for fresh water is between [0.2-1] MPa and between [0.05-0.5] MPa for atmospheric water. The values of the adhesive shear strength are given for the supports used for the experiments, i.e. bare aluminum.

Fig. 3 - 33  3 

 3333 Fig. 3 -Cracks and delamination in piezoelectric de-icing system

Fig. 4 -

 4 Fig. 4 -Mode shape of the test specimen C with a 5mm-thick ice layer around 5 kHz.

Fig. 5 -

 5 Fig. 5 -Measured quality factor of the test specimen A around 5 kHz (flexural mode).

Fig. 6 -

 6 Fig. 6 -De-icing result for the test specimen C excited around 5 kHz for 90 V

Fig. 7 -

 7 Fig. 7 -Mode shape of the test specimen B around 38 kHz (a) and de-icing result at 140V (b)

Fig. 8 -Fig. 9 -

 89 Fig. 8 -Mode shape at 3 kHz for the leading edge profile

  Fig. 11 shows the ice cracking for different values of voltages. The first cracks (Fig. 11(a)) occurred for a voltage of 100 V and a current of 1.1 A once again in the range of the estimated values. The cracks were well positioned on the node lines of the structures. Tests were performed until 150 V, voltage for which almost all the ice has been removed except on the upper part of the leading edge. This high frequency was efficient for de-icing but at the cost of a high consumed current.

Fig. 10 -Fig. 11 -

 1011 Fig. 10 -Mode shape at 31 kHz for the leading edge profile

Fig. 12 -

 12 Fig. 12 -Test of the leading edge sample with thicker and rime ice excited around 31 kHz

Fig. 13 -

 13 Fig. 13 -Apparent power versus ice thickness (a) computed for the flexural mode (b)

Table 1 -Estimated numerical results for the test specimen A (computed for frequency around 5 kHz and a quality factor of 25)

 1

	Criteria	% node for the	Voltage,	Current,	Stress in the	Maximal
		computation	V	A	piezoelectric	Displacement,
					actuator, MPa	µm
	adhesive shear strength	1%	176	0.12	21	24
	cohesive tensile strength	1%	85	0.06	10	11.5
	adhesive shear strength	2%	187	0.13	22	25
	cohesive tensile strength	2%	87	0.06	10.5	12
	adhesive shear strength	5%	208	0.14	25	28
	cohesive tensile strength	5%	93	0.06	11	12.6

Table 2 -Estimated numerical results for the test specimen C (computed for frequency around 38 kHz and a quality factor of 45)

 2

	Criteria	% node for the	Voltage,	Current,	Stress in the	Maximal
		computation	V	A	piezoelectric	Displacement,
					actuator, MPa	µm
	adhesive shear strength	1%	145	0.89	6.7	3.8
	cohesive tensile strength	1%	145	0.89	6.7	3.8
	adhesive shear strength	2%	158	0.97	7.3	4.2
	cohesive tensile strength	2%	151	0.93	7	4
	adhesive shear strength	5%	184	1.12	8.5	4.4
	cohesive tensile strength	5%	164	1	7.6	4.3

Acknowledgements:

The present work is partially funded by Toulouse University and Occitanie region. The authors gratefully acknowledge the support of Airbus Central R&T in Munich, without which the present study could not have been completed. Many of the tests presented in this article were performed in the Airbus iCORE icing wind tunnel and were supported by Vittorio Vercillo, Javier Mayen, and Alexandre Laroche.