Numerical and Experimental Investigation of the Design of a Piezoelectric De-Icing System for Small Rotorcraft Part 1/3: Development of a Flat Plate Numerical Model with Experimental Validation

Abstract

The objective of this research project is divided in four parts: (1) to design a piezoelectric actuator-based de-icing system integrated to a flat plate experimental setup and develop a numerical model of the system with experimental validation, (2) use the experimental setup to investigate actuator activation with frequency sweeps and transient vibration analysis, (3) add ice layer to the numerical model and predict numerically stresses for different ice breaking with experimental validation, and (4) bring the concept to a blade structure for wind tunnel testing. This paper presents the first objective of this study. First, preliminary numerical analysis was performed to gain basic guidelines for the integration of piezoelectric actuators in a simple flat plate experimental setup for vibration-based de-icing investigation. The results of these simulations allowed to optimize the positioning of the actuators on the structure and the optimal phasing of the actuators for mode activation. A numerical model of the final setup was elaborated with the piezoelectric actuators optimally positioned on the plate and meshed with piezoelectric elements. A frequency analysis was performed to predict resonant frequencies and mode shapes, and multiple direct steady-state dynamic analyses were performed to predict displacements of the flat plate when excited with the actuators. In those steady-state dynamic analysis, electrical boundary conditions were applied to the actuators to excite the vibration of the plate. The setup was fabricated faithful to the numerical model at the laboratory with piezoelectric actuator patches bonded to a steel flat plate and large solid blocks used to mimic perfect clamped boundary condition. The experimental setup was brought at the National Research Council Canada (NRC) for testing with a laser vibrometer to validate the numerical results. The experimental results validated the model when the plate is optimally excited with an average of error of 20% and a maximal error obtained of 43%. However, when the plate was not efficiently excited for a mode, the prediction of the numerical data was less accurate. This was not a concern since the numerical model was developed to design and predict optimal excitation of structures for de-icing purpose. This study allowed to develop a numerical model of a simple flat plate and understand optimal phasing of the actuators. The experimental setup designed is used in the next phase of the project to study transient vibration and frequency sweeps. The numerical model is used in the third phase of the project by adding ice layers for investigation of vibration-based de-icing, with the final objective of developing and integrating a piezoelectric actuator de-icing system to a rotorcraft blade structure.

1. Introduction

In-flight icing is an important problematic for helicopters. Ice accretion leads to, amongst other consequences, decrease in lift and increase in drag, torque, vibration, and imbalance [1]. All helicopters currently habilitated to fly into icing conditions are equipped with an electrothermal ice protection system for the main rotor. Those systems require large amounts of power and add a lot of weight to the helicopter. They are not suitable for small rotorcraft, which cannot provide the excess power required and the weight capacity to support those systems [2]. Rotorcraft are often used for their great flexibility and easy access to terrain and locations that can’t be reached by other means of transportation. Not being able to fly under certain conditions such as icing greatly reduces this flexibility and accessibility advantage, which is a major drawback for those vehicles. The community is in dire need of finding an alternative for small rotorcraft and is in active research to develop and integrate it. Different alternative systems have been investigated by the community. The weeping wings system is a method based upon the freezing point depressant concept [3]. Fluid with low freezing point is pumped to the leading edge of the wings and other critical surfaces of aircraft, depressing the freezing point of the water impacting. For the rotorcraft industry, Bell Helicopter investigated this type of system on UH-1 in the 1960s. Experiments showed that for a LWC of 0.8 g/m³, a standard value for atmospheric icing, the system used 0.41 kg/min of fluid down to −20 °C, which was considered effective [2]. However, the added weight for the fluid container as well as the room it requires makes it not practical for small rotorcraft. Pneumatic boots are another proven system for the aircraft industry. An inflatable boot is installed on the wings of the aircraft and when a sufficient amount of ice has accreted on the surface, air is pumped and fills the boot. By inflating, the boot breaks the ice accumulation and is then deflated. The weight added and energy consumed by this system can be considered negligible. However, its application on helicopters is impossible mainly due to the damage sustained by the boots caused by the high erosion affecting the rotating blades. Also, the cyclic nature of the system allows small amounts of residual ice to accumulate, whose roughness can cause significant aerodynamic losses [4].

As a solution to this problem, the current work investigates the use of piezoelectric actuators as a low-energy de-icing system that could be implemented on small rotorcraft. Piezoelectric actuators are devices that generate strains when subjected to an applied voltage. Different approaches using piezoelectric actuators have been explored in the past. First investigated by Ramanathan [5] and later by Palacios [6], resonating piezoelectric elements were used to generate ultrasonic surface waves to produce shear stress at the ice interface. Ramanthan used an aluminum plate on which a piezoelectric actuator patch was bonded. The actuator was excited at a frequency of 1 MHz, the resonant frequency of the actuator, to remove an ice layer accreted a few centimeters ahead of the actuator. Delamination of the ice occurred after more than a hundred seconds and thermal energy released from the actuator played an important role in it.

The method later provided interesting results for flat and shell structures. Venna et al. [7] used piezoelectric actuators located directly below the iced zones of a NACA 0012 leading edge to generate both local shear strains and normal impulse forces to de-ice, however with limited success. Unlike Ramanathan, they used the resonant frequencies of the structure instead of the resonant frequency of the actuator. The leading edge structure was fixed spanwise at its extremities and actuators positioned at the stagnation line. Ice was accumulated in a cold room. Instantaneous de-icing wasn’t observed as they forecasted. Instead de-icing was obtained after more than a minute of excitation. Similar to the de-icing mechanism explored by Venna and al., Kandagal and Venkatraman [8] partially succeeded in de-icing a cantilevered flat plate with piezoelectric actuators by exciting its resonant frequencies. Excitation of the fourth mode, at 805 Hz, resulted in instantaneous delamination of the ice at several locations on the plate, without any trace of heating from the actuators. Analysis of the results showed that shear stress generated at the interface seemed to be responsible for breaking the ice adhesion to the substrate.

Palacios [6] studied ultrasonic waves and their dispersion on a flat plate. It was showed that as the excitation frequency approached the resonant frequency of the actuator, the adhesion of the ice diminished and the plate temperature increased. It was not found what percentage of adhesion reduction was due to the heating of the plate. The author successfully delaminated an ice layer of 2.54 mm on a steel plate with a disk actuator at 28.5 kHz and 50 W, satisfying results from numerical simulations. The author tried designing actuators to excite specific mode at frequencies of 200 kHz and 500 kHz without success. The study concluded that the excitation frequency should be below 100 kHz since strains generated diminished with a frequency increase. Zhu et al. [9] followed up on Palacios research by adding discontinuities on the surface of setup similar to what Palacios used. Using the same activation methods for the actuators, they hoped that the high frequency waves would be guided by the discontinuities. Numerical results were in agreement with this, predicting increases in the strain generated by up to 400% for some places on the plate. Experimental results showed diminution of power required to de-iced by more than a factor of two.

Quilan et al. [10] performed experimental testing with a thin curved aluminum surface equipped with piezo-electric actuators applied over a rotor blade-shaped section. The rotor blade-shaped section was built out of wood and a two millimeter rubber band was installed between the section and the thin curved aluminum surface. A total of seven actuators were bonded to the aluminum surface. Wind tunnel testing showed that the most effective method of de-icing was by cascade activation. In cascade activation, two actuators were operated at once, alternating the activated piezo actuators along the span. This led to the emission of a patent for vibrating sleeves. The small thickness of the sleeve was tested and its weak strength makes it impractical for rotorcraft application. No numerical simulations were conducted during this study.

In addition to experimental testing, numerical modelling of piezoelectric de-icing systems has been investigated by the scientific community. Palacios [11,12] implemented a 3D finite element model and used it to design de-icing system experimental setups. With the help of this tool, they investigated in particular the impedance matching for excitation, which lead to promising experimental results. Boutros [13] proposed a de-icing solution using lightweight Macro Fiber Composite (MFC) actuators on the leading edge of a wing. The developed numerical model for MFC actuators was validated experimentally on cantilever beams. The resonant frequencies of the structure excited by the actuators were calculated to predict if the adhesion of ice could be overcome by the shear stress generated. The study concluded that by exciting the third mode of the structure with three MFC actuators at optimal locations, a more than sufficient shear stress could be generated to delaminate the ice. A fixed value of 1.6 MPa was used for the adhesion value of ice. Habibi et al. [14] investigated the use of piezoelectric actuators to excite shear stresses on Carbon-Fibre-Reinfoced Polymer (CFRP) as a de-icing system. They used FEM to study the wave propagation through this material. With the results they established an optimal excitation of the system by comparing the resulting shear stress with the adhesive shear strength of ice. The obtained results will be tested experimentally in a later study.

Budinger et al. [15] developed an analytical model for piezoelectric de-icing at ultrasonic frequencies. The model showed that flexural modes generate more shear stress at the ice/substrate interface than extensional modes. It also showed that Langevin type actuators, piston-like, were more suitable than patches mainly due to the resistance to mechanical stress and lower energy consumption. The model was validated experimentally and showed promising results, enough to convince the authors to pursue their study. They also showed that using this system with the addition of an icephobic coating could significantly reduce the power consumed by the actuators and at the same time lengthen their life duration [16]. Kalkowski [17] used numerical simulations to estimate the shear stress at the interface generated by piezoelectric actuators and the matching power consumption of the actuators. The simulations were conducted using a wave finite element method. The experimental investigation demonstrated that the model predicted reasonably well the power requirements for ice delamination. Yongjiu [18] performed ANSYS numerical simulations on a wing edge skin model, which revealed the feasibility to remove ice by piezoelectric ceramic actuators. Their experimental investigation showed that successful de-icing can be achieved with a voltage of 650 V and a forced frequency at 1530 Hz. This was in agreement with the results obtained with the ANSYS model. Tian et al. [19] developed with ANSYS a simplified 2D model of a vibrating flat plate. A validation of the model was done for the frequency values of the resonant modes, but not for the vibration amplitudes. While they obtained good accordance between the model and the experiments for most of the modes investigated, some discrepancies were found for some resonant modes, which were explained by the 2D simplification as well as the inaccuracy of some model parameters. They studied the vibration response to the length of the piezoelectric actuator with their ANSYS model and concluded that maximum excitation is acquired when the length of the actuator is an odd integer multiple of the half wavelength.

Harvey [20] followed by Villeneuve et al. [21] demonstrated a proof of concept for de-icing with actuator patches at structural resonant frequencies matching the first modes of the structures on which they were bonded. Total or partial de-icing was successfully achieved for a flat plate, a thinned Bell 206 main rotor blade, and a Bell 206 tail rotor blade sections. Numerical simulations were used to predict the resonant frequencies and the mode shapes to optimally position and drive the actuators.

In the present study, a numerical model of a simple flat plate was elaborated to understand phenomena involved in vibration-based de-icing and assist with the design of the piezoelectric de-icing system. The symmetry and simplicity of this structure facilitated the study and the elaboration of useful conclusions to design an optimal system and adapt it to a blade structure. A Modal Frequency Analysis was performed to compute the natural frequencies and the modal shapes of the structure. Direct-Solution Steady-State Dynamic Frequency Response Analysis was performed to define the optimal actuator positioning and phasing. An experimental setup was built with actuators located in accordance with the numerical model result analysis. Experimental testing was performed with the setup to validate the numerical model with a scanning laser vibrometer at the NRC in Ottawa. The validated numerical model and the experimental setup designed will serve as the basis for further investigations of the vibrations of the flat plate to frequency sweep excitation, and when covered with an ice layer. Numerous challenges have been observed and addressed while implementing, numerically validating, and characterizing this set-up. The numerical modal analysis implemented in this study faced some challenges related to the number and the dimensions of the required piezoelectric actuators. While it was evident that the success of the method was depending on the capability of the actuators to sustain and transfer the required power of excitation at a structural natural resonant frequency of interest, the dimensions of the actuator relative to the wavelength of vibration, their positioning, and their number, as well as the phasing of excitation, were expected to affect the integrity of some of the high frequency modes and consequently the ability to fully excite the targeted modes. This comprehensive analysis comprises the evaluation of the numerical and experimental errors due to the sometime inefficient mode of excitation.

2. Materials and Methods

2.1. Numerical Investigation for Flat Plate Experimental Setup Design

A study on a simple flat plate structure was performed in order to assess phenomena involved in structure resonant vibration-based de-icing methods and to develop design guidelines for piezoelectric-based de-icing systems. The flat plate structure was selected for its simplicity to facilitate numerical modeling and experimental testing. The comprehension and conclusions obtained from this study will be used in future studies to include ice layers and eventually consider a more complex structure with an aerodynamic profile, for the design of a piezoelectric de-icing system integrated to a helicopter blade prototype. The first step was to model a flat plate structure in order to draw basic guidelines for the integration of piezoelectric actuators into a flat plate experimental test setup that will generate cracking and delamination of ice layers. ABAQUS CAE software was selected because it handles piezoelectric type finite elements and to build on the work performed in a previous study [21]. To obtain the structure’s vibration amplitude response to different excitation scenarios, the direct steady-state dynamic analysis was used with ABAQUS CAE. In this procedure, the solution of the perturbed system is obtained by linearization from the current base state. The formulation is based on the dynamic virtual work equation shown under its discretized form at Equation (1), with $u$ and its associated derivative the displacement, velocity, and acceleration; and $M^{NM}$ the mass matrix, $C_{\left(m\right)}^{NM}$ the mass proportional damping matrix, $I^N$ the internal load vector, and $P^N$ the external load vector. It is assumed that the structure undergoes small harmonic vibration to obtain the steady-state harmonic response and since it is a perturbation procedure, the change from the base state is defined by the step’s load and response.

$$\delta u^n\left\{M^{NM}{\ddot{u}}^M+C_{\left(m\right)}^{NM}{\dot{u}}^M+I^N-P^N\right\}=0$$

(1)

In this step, the structural damping was used, which is defined in the model by Equation (2) with $F_D^N$ the damping forces, s the structural damping, and $I^N$ the forces caused by stressing (excitation) of the structure.

$$F_D^N=is{I}^N$$

(2)

A flat rectangular plate with the dimensions of 0.5 m × 0.2 m and 0.0016 m thick was modeled in ABAQUS as shown in Figure 1. The plate is made of Stainless Steel 304 with a mass density of 8000 kg/m³, a Young’s modulus of 193 GPa, and a Poisson’s ratio of 0.27. An inherent structural damping coefficient of 0.01 was used as a standard value for the preliminary simulations. This coefficient was readjusted for the assembled system when piezoelectric actuator testing was performed to obtain a more accurate value. As can be observed in Figure 1, a fully clamped boundary condition (fixed translational and rotational degrees of freedom) was applied to the longest edges of the plate. This configuration was used to fit an existing support available in the laboratory. Linear S4R shell elements were used to mesh the plate. A convergence study was performed in order to ensure that the numerical implementation accurately simulates the dynamic behavior of the plate.

The structure’s frequency response was estimated using the direct steady-state dynamic analysis for an optimal positioning of the actuators to generate maximum vibration levels. The dynamic analysis succeeded a modal frequency analysis performed from 0 to 2500 Hz, which calculated the resonant frequencies of the structure. This range was selected due to de-icing obtained at those frequencies during a previous study [21]. Forty-six resonant modes were found within this frequency range, with the first mode at 177 Hz. Examples of mode shapes are shown in Figure 2.

A load of 1 N, defined as a nodal force in the z-axis direction (direction of thickness of the plate), was applied first at the center of the plate. It is well known that the 1 N value is arbitrary and facilitates comparison for each potential actuator configuration. Since the direct steady-state dynamic analysis is a linear perturbation procedure, this load value has no significance and the results can, as expected, be extrapolated linearly to any load values. The direct steady-state dynamic analysis was run in the same frequency range as the modal analysis, from 0 to 2500 Hz. As expected, the results obtained show that all of the resonant modes with an anti-node located at the center of the plate were properly excited and resulted in significant modal displacements and accelerations (Figure 3). On the other hand, as expected, the resonant modes without an anti-node at the center of the plate were not properly excited and the resulting vibration shape was not consistent with the resonant mode shapes, and levels were negligible in terms of displacements and accelerations compared to other modes (Figure 4).

Three different simulations were run with the 1 N load applied at three different positions lengthwise on the center anti-node of mode 3 (Figure 3a): at the edge of the anti-node, at one third of the anti-node, and on the center of the anti-node. The load was applied at the center of the plate in the width direction for the three cases. Displacements at those three positions are presented in Figure 5 for the three load cases. The results show as expected that the optimal load positioning is at the center of an anti-node. Simulations were run for mode 9 (Figure 6) for three load configurations: 1 N load at center of anti-node 1, 1 N load at center of anti-node 3, and 0.33 N loads on anti-nodes 1, 3, and 5. Results, presented in

Table 1. Displacement for three load position configurations at 640 Hz (Mode 9).

Mode 9	Displacement at Center of Anti-Node 1 (µm)	Displacement at Center of Anti-Node 3 (µm)	Displacement at Center of Anti-Node 5 (µm)
1 N load at anti-node 1	67	67	67
1 N load at anti-node 3	66	66	66
0.33 N load at anti-nodes 1, 3 and 5	66	66	66

, demonstrate as expected that the load can be positioned on any anti-node for the same results and that it can even be divided between different anti-nodes for equivalent results. This was important to validate since the load on each piezoelectric actuators can be reduced by distributing it on multiple actuators to prevent their damaging or solicitations close to their operational limit, which can greatly reduce their lifespan. On the other hand, this also indicated that the number of actuators can be reduced to limit the cost and complexity of their integration.

As it can be observed in Figure 6, each antinode is out of phase with its neighboring antinodes, forming a sine wave. To evaluate and validate the impact of load phasing on vibration, two load configurations were applied under a direct-solution steady-state dynamic numerical simulation analysis for mode 9 at 640 Hz (Figure 6). First, a load of 0.2 N was applied at the center of each anti-node, for a total of 1 N. The second configuration consisted of loads of 0.2 N applied at the center of anti-nodes 1, 3, and 5 and loads of −0.2 N applied at the center of anti-nodes 2 and 4, so that each successive load is out of phase compared to the ones next to it.

Table 2. Displacement for three load position configurations at 640 Hz (Mode 9).

Mode 9	Displacement at Center of Anti-Node 1 (µm)	Displacement at Center of Anti-Node 2 (µm)	Displacement at Center of Anti-Node 3 (µm)	Displacement at Center of Anti-Node 4 (µm)	Displacement at Center of Anti-Node 5 (µm)
0.2 N loads on anti-nodes	12	12	12	12	12
0.2 N and −0.2 N loads on anti-nodes	66	66	66	66	66

presents the displacements at the center of the anti-nodes for the two configurations. When the loads were applied in opposite phasing, the computed displacements were similar to the three configurations of Table 1, which is expected since the same total force load is applied. However, when all the force loads were in the same direction, the computed displacements decreased by a factor of 5.5, meaning that some forces were acting, as expected, against the optimal deployment of the mode.

Those simulations allowed to obtain useful information for actuators positioning and also to confirm some of the results obtained in a previous study [21]. These results were used to design the actuators integration to the flat plate setup. Testing performed in these previous studies has demonstrated that de-icing for the flat plate was mainly obtained for frequency sweeps between 0 and 2500 Hz. In those experiments, an ice layer of 40 × 450 mm was accreted at the center of the flat plate, simulating an ice layer accreted at the leading edge of an airfoil. To generate ice breaking in the ice layer, maximum vibration and stress must be generated on the centerline of the flat plate. Simulations have shown that piezoelectric actuators must be placed at the center of the plate in the width direction in order to be positioned at the center of the anti-nodes most susceptible to generate maximum displacement and stress in the ice layer.

To position the actuators lengthwise, positions of anti-nodes in the length direction were studied for the modes obtained from 0 to 2500 Hz. Modes with two lines of anti-nodes were neglected. Every mode has half anti-nodes located on each extremity of the plate (Figure 3 and Figure 4), which was selected to position actuators. Every mode with an odd number of anti-nodes has an anti-node located at the center of the plate, meaning an actuator at that position will contribute to the optimal deployment of half of the modes investigated. Anti-node positions lengthwise is different for each mode, meaning that no other position is ideal for a majority of the modes investigated. It was therefore decided to position an actuator at both edges of the plate and at the center of the plate. It was also decided to install two additional actuators at anti-node positions of mode 4 (16 cm from each edges), which has no anti-node on the center of the plate. Those two actuators will allow to activate mode 4 and can contribute to excite other modes. Final actuator positioning is shown in Figure 7 as yellow arrows.

Physik Instrumente P-876.A15 patches actuator were used for the flat plate setup due to their availability at the laboratory and for their demonstrated performances in a previous study [21]. The patches were about double in length compared to their width (61 mm in length by 35 mm wide with a total thickness of 0.8 mm). Simulations were conducted to define in which direction should the patches be installed for maximum mode deployment. A solid deformable 3D part was created to model the P-876.A15 patches. The part was installed at the center of the plate in the width direction for a first analysis (Figure 8a) and in the length direction for a second one (Figure 8b). A pressure force of 1 N was applied to both width faces of the actuator, which simulated the contraction/elongation of the actuator when vibrating, and the displacement results were compared for both cases. Results at Figure 9 show that the optimal actuator direction varies between modes. The same process was repeated with five patches installed at positioned defined for the experimental setup (Figure 7) in both directions, and simulations were run for the first five modes with the optimal phasing for each mode. Displacement amplitudes at the edge of the plate for each mode showed that the optimal placement for the actuators is in the width direction (Figure 10).

2.2. Numerical Model of Experimental Setup with Forced Vibration Generated by Piezoelectric Actuator Patches

The results of the numerical investigation presented in the previous section has allowed to confirm the conclusions obtained in the previous study [21] as well as to obtain new information for the integration of the actuators to a flat plate structure to optimize de-icing of an ice layer. A numerical model integrating the simplified numerical modeling study conclusions was developed for the experimental setup incorporating piezoelectric actuators to generate forced vibration on the plate. The five actuator patches were positioned as follows: at both extremities of the plate, at its center and at positions of anti-nodes of mode 4, which correspond to 16 cm from each edge. The actuators were all installed in the width direction with their centerline centered with the plate’s width (Figure 11). The actuators were modeled as 3D extruded deformable solid and tied to the flat plate with a tie constraint.

PIC255 material was created in the software. ABAQUS uses the e-form constitutive equation to couple stress and electrical field, and to resolve numerical piezoelectric analysis with the finite element method. The material properties required for these equations are the stiffness, the electrical permittivity, and the piezoelectric coefficient. Those values are dependent on the piezoelectric material characteristics and are provided by the manufacturer. The elastic component for stiffness is anisotropic and the values input are shown in this stiffness matrix, as follows:

$${\mathrm{c}}^E=\left[\begin{array}{cc}\begin{array}{ccc}13.27& 8.667& 8.563\\ 8.667& 13.27& 8.563\\ 8.563& 8.563& 11.92\end{array}& \begin{array}{ccc}0 & 0& 0\\ 0 & 0& 0\\ 0 & 0& 0 \end{array}\\ \begin{array}{ccc}0 & 0& 0\\ 0 & 0& 0\\ 0 & 0& 0 \end{array}& \begin{array}{ccc}2.128& 0& 0\\ 0& 2.128& 0\\ 0& 0& 2.299\end{array}\end{array}\right]\ast {10}^{10} \left(\mathrm{N}/{\mathrm{m}}^2\right)$$

(3)

The electrical permittivity (dielectric component), which is orthotropic, is:

$${\mathsf{\epsilon }}^S=\left[\begin{array}{ccc}8.28& 0& 0\\ 0& 8.28& 0\\ 0& 0& 6.72\end{array}\right]\ast {10}^{-9} \left(\mathrm{F}/\mathrm{m}\right)$$

(4)

The piezoelectric coefficient, an anisotropic constant, can be described by the stress related coefficient or the strain related coefficient. The strain related coefficient matrix was input in the model and presented below:

$$e=\left[\begin{array}{cc}\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ -6.73& -6.73& 15.68\end{array}& \begin{array}{ccc}0& 13.14& 0\\ 13.14& 0& 0\\ 0& 0& 0\end{array}\end{array}\right] \left(\mathrm{N}/\mathrm{Vm}\right)$$

(5)

The coupling relationship between the piezoelectric element properties can be expressed under the matrix form by Equation (6), giving the stress (T), and (7) giving the electric induction (D). In Equation (6), the first term represents the elastic relationship, with $c^E$ the constant electric field stiffness matrix and S the strain, while the second term is the coupling for the direct piezoelectric effect with $e$ the constant strain piezoelectric coefficient and E the electric field.

$$T=c^{E}S-e^{t}E$$

(6)

In Equation (7), the first term is the reverse piezoelectric effect and the second term is the dielectric relationship, with ${\epsilon }^S$ the constant strain permittivity.

$$D=eS+{\epsilon }^{S}E$$

(7)

The mass density of 7800 kg/m³ and the material orientation of the non-isotropic components were also defined. The actuators were meshed with standard linear C3D8E 8-node piezoelectric brick, which allowed coupling between mechanical and electrical properties by including displacement and electrical potential degrees of freedom. Electrical boundary condition was applied to the top and bottom face of each actuator to generate the deformation of the patches. On the bottom face, the electrical potential was set constant to zero and on the top face the potential was set to the value of the voltages applied during the different experiment testing. This simulated the two electrodes on the top and bottom of the actuator patch.

2.3. Flat Plate Experimental Setup

The setup consists of a 304 stainless steel flat plate, 1.6 mm thick, with five actuators bonded to it to excite at its resonant frequencies. The plate was cut in a rectangle shape of 500 mm by 504.8 mm. Ten holes were drilled lengthwise in each side of the plate. The two sides were placed between two massive steel 44 w 50.8 mm (2”) thick blocks. Ten screws and bolts were used on each side to tighten the plate with the blocks. The screws were tightened to 165 lb-in with a torque wrench. This was designed and fabricated to recreate the fully clamped boundary conditions (fixed translational and rotational degrees of freedom) applied on the longest edges of the plate in the numerical model (Figure 1). It is necessary to reproduce the boundary conditions of the model as accurately as possible and maximize vibration of the plate. The middle section of the plate, which is 500 m × 200 mm long, is free to vibrate and represents the plate in the numerical model (Figure 12).

Five Physik Instrumente (PI) P-876.A15 actuator patches were installed on the lower surface of the plate. The patches, which are 61 mm in length by 35 mm wide with a total thickness of 0.8 mm, were installed at the positions shown on the experimental setup in Figure 13, as defined from the results of the numerical simulations in Section 2.1. The actuators were bonded to the surface following a similar process as used to install strain gages. The surface was cleaned with a CSM-2 degreaser and was sandpapered at 45° in both direction with a 400 grit. The surface was cleaned two more times, first with conditioner and then with neutralizer. The surface of each actuator was also cleaned with acetone. Epo-tek 353ND glue, as recommended by the piezoelectric actuators manufacturer, was applied both on the surface of the actuators and the plate. The actuators were carefully applied on the plate at the selected positions and pressure was applied on the actuators by depositing a heavy mass on each actuator for 72 h.

To drive and monitor the piezoelectric actuator system, an electrical system was required. The system used the laser vibrometer described in the next paragraph to provide the electrical signal of excitation to the actuators. The driving signal for the actuators was amplified using an Amp-line AL-1000-HF-A amplifier with a range of 50 to 1000 V to generate the required vibration level. Since the operational voltage of the actuators is −250 to 1000 V, an Amp-line AL-100DC power source was used to offset the voltage to allow testing close to the limit of the actuators. By offsetting the voltage to values around 400 or 500 V, an alternative voltage of 1000 to 1250 Vpp could be applied without compromising the integrity of the actuators. The voltage applied to the actuators was measured with a Fluke 105B oscilloscope (Figure 14).

A Polytec PSV-300 Laser Scanning Doppler vibrometer (sensor head OFV 056) available at the CNRC in Ottawa was used to measure the structure eigen modes and deflection shapes. The scanning head of laser vibrometer was positioned on a perpendicular direction at a predetermined distance from the flat plate (Figure 15). The excitation signal was generated using the Polytec signal generator card (PCI-6711) as stated in the previous paragraph.

3. Results

3.1. Numerical Model Validation

Numerical simulations and experimental tests were performed using the flat plate setup to validate the numerical model. The displacement, velocities, and accelerations of the plate were measured using the Polytec Scanning Laser Vibrometer PSV-300 (data acquisition board PCI-4452) at a sampling rate of 2.56 times the maximum frequency of interest. A grid of measurement points was defined over the surface of the flat plate to measure the frequency response and the modal shapes of the panel. The experimental tests were performed under the piezoelectric actuators excitation using sinusoidal signals at fixed frequencies to obtain steady-state mode deployment to compare with the numerical analysis results. The experimental results, given by the laser vibrometer software, present the upwards displacements in different shades of red, with lighter shades for higher magnitudes, and downwards displacements in different shades of green, also with lighter shades for higher magnitudes (see Figures 16, 17, 19 and 22). The operating voltage of the piezoelectric actuator was set to 200 Vpp, both experimentally and in the numerical model. Testing and numerical simulations were performed for the structural modes in the frequency range of interest of 0 to 1000 Hz that were proven efficient for de-icing in a previous study [21].

3.1.1. Frequency Analysis Validation

A comparison between the numerically estimated and experimentally measured modal frequencies of the structure is presented in

Table 3. Comparison of resonant frequencies obtained with numerical model and experimental results.

Resonant Mode (#)	Numerical Model (Hz)	Experimental Results (Hz)	Difference (%)
1	200	177	13
2	206	192	7
3	227	216	5
4	270	258	5
5	344	327	5
6	453	424	7
7	570	501	14
8	570	557	2
9	581	527	10
10	613	-	-
11	671	604	11
12	752	713	5
13	763	-	-
14	899	799	13
15	942	913	3

. As expected, the modal frequencies are independent of which actuators are excited and at what voltage. However, it is worth mentioning that two modes, mode 10 and 13, were not activated or measured during the experimental testing. For the other modes, the maximum observed difference between experimental and numerical results was 14%.

It is possible to divide the resonant modes in two distinct groups within this frequency range: the modes with one line of anti-nodes (Figure 16) and the modes with two lines of anti-nodes (Figure 17). The modes with one line of anti-nodes are the modes in bold in Table 3. The actuators have been positioned to optimize displacement at the center of the plate. Since the system has not been designed to excite the modes with two lines of anti-nodes, those modes are barely observable experimentally. This explains the higher difference between the numerical simulation and experimental results for those modes as well as the two missing modes. When excluding those modes (mode 7, 9, 10, 11, 13, and 14), the difference between experimental and numerical results falls under 7% for all but the first mode, which was expected since it is well known that the first mode of vibration is highly influenced by the boundary conditions. Those modes are no longer considered in this study since the system has been designed to maximize displacements at the center of the width of the plate, thus optimizing the modes with one line of anti-nodes. The results have also shown that there is no or negligible coupling between the modes investigated and that correlation techniques such as the Modal Assurance Criterion or other similar techniques were not required. This could be expected due to the simplicity of the structure.

3.1.2. Direct-Solution Steady-State Dynamic Analysis Validation

Damping characterizes the dissipation of elastic energy in a material or structure under cyclic solicitation. Theoretical damping values tend to be difficult to determine. With the experimental setup, it is possible to determine experimentally the value of the damping for the structure. In order to do so, the half-power method, also called −3 dB method, is used. In a frequency response function, the damping is proportional to the width of the resonant peak of a structural mode. To obtain the damping experimentally, sweeping from 160 to 1000 Hz was performed with each actuator individually at 200 Vpp, except with actuator 1 because it got debonded from the plate during transportation of the test setup between AMIL and NRC laboratories. The frequency response function was measured using the Polytec laser vibrometer and the loss factor estimated using the half-power method. The loss factor (η), which is the inverse of the damping factor (Q), was calculated at each resonant frequency for each actuator activation. Loss factor results are presented in Figure 18. Mode 15 was not considered for the calculation of the damping loss factor. The experimentally measured mode shape shown in Figure 19 seems to be in disagreement with the numerically simulated mode shape. The experimental mode shape obtained at 931 Hz seems to be the combination of two different modes, mode 15 and mode 16, which were predicted numerically at 942 Hz and 1002 Hz, respectively. Those two modes are close in frequency and the −3 dB frequencies for mode 15 cannot be distinguished from the vibration of mode 16. Since these two modes seem to be highly coupled, the half-power method would have not been appropriate to be used for this case. The remaining modes in the frequency range of interest, measured and identified using the laser vibrometer, were sufficiently decoupled (no double peaks observed), which allowed the use of the half-power bandwidth method to estimate the damping loss factor. For the piezoelectric actuator 3, which is positioned at the middle of the plate, only modes with odd number of anti-nodes are activated correctly. The absence of an anti-node at the center of the plate for the other modes resulted in very small displacements at those frequencies when the actuator 3 was driven and are therefore not considered here. The loss factor obtained for all actuator excitations were of the same order of magnitude. However, it slightly varied from mode to mode, with values ranging from 0.004 to 0.008. The average loss factor obtained from all actuators activation and all resonant modes was 0.006 ± 0.002. This value was used as a structural damping input for the direct-solution steady-state dynamic analysis.

Single Piezoelectric Actuator Excitation

Fast-scan measurements were performed using the laser vibrometer for each of the nine structural modal frequencies selected in Section 3.1.1. The fast-scan operation is based on a steady-state vibration excitation at each selected structural modal frequency that corresponds to the Direct Solution Steady-state Dynamic Analysis that were performed numerically. Fast-scan was performed with one actuator activated at a time at each structural resonant mode and for each actuator, except actuator 1. A direct-solution steady-state dynamic analysis was performed for each test case with the numerical model. Examples of experimental and numerical displacement results at mid-width along the length of the plate are shown in Figure 20. Due to the symmetry of the plate and position of actuators 2 and 4, numerical results were identical when excited by one or the other actuator. Four test conditions were repeated twice and the repeatability of the experimental results was within an approximate range of 20%.

Table 4. Average displacements at edge of plate for activation of a single actuator.

Mode	Actuator 2		Actuator 3		Actuator 4		Actuator 5
(#)	Num (µm)	Exp (µm)	Num (µm)	Exp (µm)	Num (µm)	Exp (µm)	Num (µm)	Exp (µm)
1	-	-	-	-	-	-	-	-
2	99	160	-	-	99	143	204	253
3	105	72	168	226	105	134	152	190
4	142	178	-	-	142	189	82	113
5	65	65	113	187	65	67	41	55
6	52	82	-	-	52	83	17	28
8	77	72	76	72	77	74	2	2
12	32	17	-	-	32	17	6	3
15	19	11	44	22	19	11	10	3

and

Table 5. Average displacements at anti-nodes for activation of a single actuator.

Mode	Actuator 2		Actuator 3		Actuator 4		Actuator 5
(#)	Num (µm)	Exp (µm)	Num (µm)	Exp (µm)	Num (µm)	Exp (µm)	Num (µm)	Exp (µm)
1	72	121	71	112	72	90	66	54
2	-	-	-	-	-	-	-	-
3	86	45	137	158	86	88	113	129
4	121	117	-	-	121	129	71	77
5	56	46	98	133	56	47	35	41
6	43	63	-	-	43	66	14	22
8	59	54	59	53	59	54	2	2
12	26	15	-	-	26	13	5	2
15	14	8	33	16	14	8	7	3

present the numerical simulation and the experimental testing average displacements on the edges of the plate and at center of anti-nodes for each mode. The results obtained for excitations with the actuator 2 and 4 confirm the symmetry of the vibration predicted by the numerical model. Except for mode 3, the experimental results obtained with these actuators are within the experimental repeatability of the test. For mode 3, the experimental results obtained with actuator 2 are abnormally different from experimental results obtained with actuator 4, and when compared to numerical results. More details are provided in the discussion (see Section 4.1).

Activation of Multiple Piezoelectric Actuators

This section presents the results obtained when multiple actuators were driven at the same time. As observed during the preliminary simulations, the piezoelectric actuators must be driven with the proper phasing for optimal mode deployment. Each anti-node is out of phase with its neighboring anti-nodes, and actuators must be activated with the same phasing as the anti-node it is positioned onto. For these tests, piezoelectric actuator 1 was bonded in its place directly at NRC with instant glue in the best way possible so it has as few impacts on the results as possible.

Table 6. Average displacements at edge of plate for activation of multiple actuators.

Mode	Actuator Activation	Displacement
(#)	Phase/Out of Phase	Num (µm)	Exp (µm)
2	1-2/4-5	491	452
3	1-5/2-3-4	510	464
4	1-4/2-5	368	345
5	1-3-5	199	186
5	1-3-5/2-4	293	293
6	1-2/4-5	140	118
6	1-4/2-5	75	113
8	1-3-5/2-4	159	110
15	1-3-5	28	19
15	1-3-5/2-4	62	34

and

Table 7. Average displacements at anti-nodes for activation of multiple actuators.

Mode	Actuator Activation	Displacement
(#)	Phase/Out of Phase	Num (µm)	Exp (µm)
1	1-3-5	230	213
1	1-2-3-4-5	335	282
3	1-5/2-3-4	405	358
4	1-4/2-5	296	249
5	1-3-5	166	194
5	1-3-5/2-4	241	225
6	1-2/4-5	112	142
6	1-4/2-5	61	83
8	1-3-5/2-4	119	90
15	1-3-5	20	13
15	1-3-5/2-4	46	23

present the average displacement at the edges of the plate and at anti-nodes for each mode for the different actuator activation tested. Actuators that are out of phase from the other actuators were separated by a backslash when marked in the table. For example, in activation 1-2/4-5, the actuators 1 and 2 are activated with a 200 Vpp loading and actuators 4 and 5 are activated out of phase with a −200 Vpp loading. During experimental testing, the connections of the actuators out of phase were inversed with the positive wires connected to the negative connectors and negative wires in positive connectors. Only optimal activation was tested, in addition to three tests with actuators 1, 3, and 5 activated in phase.

4. Discussion

4.1. Numerical Model Validation

Table 8. Absolute relative discrepancy of the numerical and experimental displacements at the edge of the plate for activation of a single actuator.

Mode (#)	Actuator 2 (%)	Actuator 3 (%)	Actuator 4 (%)	Actuator 5 (%)
1	-	-	-	-
2	38	-	31	19
3	45	26	22	20
4	20	-	25	27
5	0	40	3	25
6	37	-	37	39
8	7	6	4	0
12	88	-	88	100
15	73	100	73	233

and

Table 9. Absolute relative discrepancy of the numerical and experimental displacements at the center of anti-nodes for activation of a single actuator.

Mode (#)	Actuator 2 (%)	Actuator 3 (%)	Actuator 4 (%)	Actuator 5 (%)
1	40	37	20	22
2	-	-	-	-
3	91	13	2	12
4	3	-	6	8
5	22	26	19	15
6	32	-	35	36
8	9	11	9	0
12	73	-	100	150
15	75	106	75	133

present the absolute value of the relative discrepancy in the predicted displacements obtained with the numerical model when compared to the experimental results at the edges of the plate and at the anti-nodes, respectively. The results predicted by the numerical model are within 40% of the results obtained experimentally for the first seven modes with an average of 20%, except for mode 3 when excited with the piezoelectric actuator 2. This deviation can be explained by the experimental repeatability of the test, as well as by the boundary conditions and the bonding of the actuators. In the numerical simulation, they are both considered ideal, which is not the case for the experimental setup despite all considerations taken to build the set-up. The material parameters input for the stainless steel are generic values for stainless steel obtained in literature, which could slightly influence the results. Also, PIC255 material values given by the manufacturer are general values for the actuator model and can slightly vary for each patch individually. For mode 3 with piezoelectric actuator 2, the measured displacement is unusual, with a difference of 91% at anti-nodes. When compared with actuator 4 for the same mode, the results are outside the experimental variability, which is also not in agreement with the symmetry predicted and observed for all other modes. For this reason, this result is deemed questionable and should be retested.

For modes 12 and 15, the results predicted by the numerical model were higher than those measured by the laser vibrometer. At those frequencies, the anti-nodes become much smaller compared to the size of the actuator patches. This creates a situation where the actuators act on the plate outside their anti-nodes, and for those modes, also act on neighboring anti-nodes. As was observed during the preliminary simulations, each anti-node is in counter-phase, which lead to a less ideal activation of the mode by the actuator due to their size (Figure 21). A slight loss of symmetry is also observed for those two modes that can result from the non-negligible stiffness added, which affect the modal shape. Figure 22 shows very narrow anti-nodes and even narrower half anti-nodes at edges. It also shows a big imperfection at one of the edges. This demonstrates that the size of the actuators is a very important criterion for design. When the dimensions of the actuator are equal or larger than half the wavelength of vibration of the structure, it is expected that (locally) the high order mode shapes and (globally) the behavior of the structure could be greatly affected. The vibration obtained experimentally is less symmetric and ideal for these reasons. These results are in accordance with the numerical predictions obtained by Tian et al. [19]. They concluded that the ideal length of an actuator patch is an odd integer multiple of the half wavelength of a mode, and an odd integer of one is optimal. The results in this paper extrapolates those conclusions by demonstrating the negative effect of exceeding the size of half wavelength of the mode. The numerical model does not seem to completely account for this, which is why the numerical results are all overpredicted compared to the experimental results. The displacements are also much smaller compared to other modes, below 22 µm, and the macroscopic equation used by the numerical model could be less representative of the phenomenon at this scale. With those weak displacements, small differences also result in much higher error percentages. Finally, some coupling effect between mode 15 and mode 16 was observed for some of the tests, which could have impacted the comparison between numerical model and experimental results for this mode. The model gives less accurate prediction when the actuator positioning and activation is not optimal and the resulting vibration shape are slightly less symmetric than the shape of the resonant mode predicted by the frequency analysis, leading to overprediction of the numerical model. The displacement results are still satisfying for the numerical model, with a difference in the order of the micrometers between the model and experimental results for those two modes. Moreover, the numerical model is created to design an optimal de-icing system and is not intended to be used for those cases.

The absolute values of the discrepancy obtained between the numerical and experimental results when multiple actuators are activated (

Table 10. Absolute relative discrepancy of the numerical and experimental displacements at the edge of the plate for activation of multiple actuators.

Mode	Actuator Activation	Discrepancy
(#)	Phase/Out of Phase	(%)
2	1-2/4-5	9
3	1-5/2-3-4	10
4	1-4/2-5	7
5	1-3-5	7
5	1-3-5/2-4	0
6	1-2/4-5	19
6	1-4/2-5	34
8	1-3-5/2-4	43
15	1-3-5	47
15	1-3-5/2-4	82

and

Table 11. Absolute relative discrepancy of the numerical and experimental displacements at the center of anti-nodes for activation of multiple actuators.

Mode	Actuator Activation	Discrepancy
(#)	Phase/Out of Phase	(%)
1	1-3-5	8
1	1-2-3-4-5	19
3	1-5/2-3-4	13
4	1-4/2-5	19
5	1-3-5	14
5	1-3-5/2-4	7
6	1-2/4-5	21
6	1-4/2-5	27
8	1-3-5/2-4	32
15	1-3-5	54
15	1-3-5/2-4	100

) are similar to a value when a single actuator is activated. Maximum difference for the first seven modes is 43% with an average difference of 17%. For modes 12 and 15, differences are also similar with the same phenomenon occurring. The maximum difference of the numerical model is established at 43% for modes efficiently activated with piezoelectric actuators properly positioned and with proper sizing.

4.2. Piezoelectric Actuator Activation Efficiency

The actuators are positioned in order to activate optimally different modes. The actuators on both edges are always positioned on a half anti-node, which means they can activate all of the 9 resonant modes investigated (Figure 16). However, since it is a half anti-node, the center of the actuator never matches the position of maximum vibration like for the other anti-nodes. For modes 3 to 15, the displacements, both predicted numerically and obtained experimentally, are smaller than for the actuators that are positioned directly at the center of the anti-nodes (Table 4 and Table 5). Moreover, starting from mode 5, both the displacements predicted and measured diminishes greatly compared to the displacements obtained with the other actuators, even when actuators 2 and 4 are not positioned directly at the center of anti-nodes of a mode. As explained in a previous section, the behavior of the structure could be greatly affected when the dimensions of the actuator are equal or larger than half the wavelength of vibration of the structure. The half anti-nodes on each side, corresponding to a quarter of the wavelength of vibration, is affected at even lower frequencies, which explains these smaller displacements for actuator 5 at those modes. For mode 1, the mode is composed of one anti-node over the complete length of the plate, which explains why the model predicts similar displacements for all actuators. For mode 2 (Figure 16a,b), which consists of only a half anti-node at each extremity of the plate, actuator 5 is positioned as close as possible to the most efficient position, which explains the highest numerical and experimental displacements for this mode.

For mode 3 (Figure 16c,d), actuator 2 and 4 are not positioned optimally but can still activate the resonant mode. When activated on their own, the numerical prediction for the anti-node displacement is 86 µm compared to 137 µm for actuator 3 (Table 5). When the five actuators are activated with the optimal phasing (1-5/2-3-4), the prediction increases to 405 µm (Table 7). For five times the power, the displacement at anti-nodes is only increased by three times compared to actuator 3. For mode 4, the actuator 2 and 4 produced a maximum displacement of 121 µm with the numerical model. With activation of 1-4/2-5, this increases to 296 µm, 2.5 times higher for four times the injected power. The actuators are not all optimally positioned for these modes, which explains efficiency lost. An actuator that is positioned to activate a mode can also be used during other mode activation and help increase displacement, but with a loss of power efficiency depending on its position compared to anti-nodes center.

5. Conclusions

A numerical model has been developed for a simple flat plate structure. Preliminary simulations have allowed the design and integration of piezoelectric actuators for de-icing research to a flat plate experimental setup. They have provided the necessary information for actuator positioning and phasing of the actuator patches for optimal mode deployment. A numerical model of the final setup was elaborated with the optimal positions of the piezoelectric actuator patches obtained from the preliminary simulations and meshed with piezoelectric elements. A physical experimental setup was also built in accordance with the numerical model for experimental validation using a scanning laser vibrometer. The numerical model was used for resonant frequency and mode shape calculations with a frequency analysis. The experimental setup confirmed the results with a maximum variation of 14%.

Following this investigation, direct steady-state dynamic analysis was performed to predict vibration amplitudes of the plate. Electric boundary conditions were applied to the actuators, and vibration of the plate resulted from the piezoelectric actuators excitation. The experimental results validated the model when the plate was optimally excited with an average of error of 20% and a maximal error obtained of 43%. This deviation was explained by the experimental repeatability of the test, as well as by the assumption of the numerical model for boundary conditions and bonding of the actuators. However, when the plate was not correctly excited for a mode, the prediction of the numerical data were not as accurate. In this study, piezoelectric actuators were too large for two of the resonant modes investigated, which resulted in a less symmetric and optimal vibration of the plate experimentally. This was not totally accounted for by the numerical model and explains the overprediction of vibration amplitudes of the model for those two modes. This is not a concern since the numerical model is developed to design and predict optimal excitation of structures for de-icing purposes. The results also confirmed the conclusions obtained with the numerical model on optimal phasing of the piezoelectric actuators in function of each resonant mode shape.

This comprehensive study allowed to develop a numerical model for a simple flat plate and a better understanding and optimization of the actuators driving strategies in preparation for the design and instrumentation of an experimental piezoelectric-based de-icing system for small rotorcraft. The experimental setup designed is used in the next phase of the project to study transient vibration and frequency sweeps [22]. The numerical model is used in the third phase of the project by adding ice layers for numerical and experimental investigation of vibration-based de-icing [23], with the final objective of developing and integrating a piezoelectric actuator de-icing system to a rotorcraft blade structure.