# Experimental Study of a Piezoelectric De-Icing System Implemented to Rotorcraft Blades
## 3 Numerical Model  
A similar ABAQUS FEA-based numerical model, as developed and validated in the previous investigations by [17; 19], was developed to replicate the structure and boundary conditions of the blade-based testing setup. The blade geometry was designed to replicate as accurately as possible the experimental setup. An extruded shell part was created using NACA0012 coordinates to generate the profile with a chord of 13.3 cm (5.25 in). The NACA profile was cut to 4.5 cm, which corresponds to the vibrating leading-edge part of the experimental setup, known to accumulate most of the ice. The resulting profile was extruded to a span of 22.9 cm (9 in). The thickness of the shell was set to 0.8 mm. Aluminum 2024-T4 material was used, with a mass density of 2780 kg/m\({}^{3}\), Young's modulus of 73 GPa and a Poisson's ratio of 0.33. The six actuator patches were modeled using piezoelectric elements with a curved profile to fit the shape of the blade, consistent with the experimental setup (Figure 4). The patches were positioned on the blade part to correspond to their actual position in the physical setup. The ice layer was created as a 3D homogenous deformable extrusion in the numerical model. The thickness was set to 6.35 mm (1/4 th of an inch)--the limit accepted on the rotorcraft blade in de-icing. Ice material properties were set to a density of 897 kg/m\({}^{3}\), Young's modulus of 9 GPa and a Poisson's Ratio of 0.31. The shape of the ice was based on the NACA 0012 airfoil profile as the ice tends to preserve the shape of the leading edge under this accumulation regime (Figure 5). Ice accumulation can be altered due to the rotation of the profile, which is something to consider for the numerical model of the rotating blade. With the increase in tangential velocity along the radius, more ice tends to accumulate at the tip of the blade than at the root. However, experimentation showed that for the low rotation speed and small radius of the rotating blade setup, the rotation had a negligible effect on the ice accumulation and the ice shape was similar to the one obtained during fixed tests (Figure 5b), also showing that the test dummy did not affect the ice accumulation. For this reason, the ice part was not modified for numerical analysis of the rotating setup. A convergence study was performed, and the number of elements was set to 26,400 C3D8R 8-node linear brick elements in the ice layer and 8664 linear S4R shell elements for the profile. To add the effect of rotation, a radial force load was applied to \begin{table}
\begin{tabular}{c c c c c c} \hline \hline \multirow{2}{*}{**Test**} & **Wind Speed** & **Air** & \multirow{2}{*}{**LWC** (**g/m\({}^{3}\)**)} & **MVD** & **Rotation** \\  & **(m/s)** & & **(\({}^{\circ}\)C)** & **(g/m\({}^{3}\))** & **(\(\mu\)m)** & **Speed** \\ \hline Static & 20 & \(-7\) & 0.8 & 20–30 & 0 \\ Rotating & 20 & \(-7\) & 0.8 & 20–30 & 600 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Test Conditions. Figure 4: Piezoelectric actuator patches attached to the curved blade shell with their respective identification numbers. the ice layer, in an additional general static step, determining the initial pre-strained state for the frequency and direct steady-state dynamic simulations. The force was calculated from the ice mass of 11 g, rotation speed of 600 RPM and average radius of 16 cm, or a total force of 7 N along the radius. Even if this method has not been validated experimentally on its own, it is a commonly used representation [25]. This method allowed a quick and reasonable estimation of the additional shearing provided by the centrifugal force, but a deeper investigation to validate this method experimentally should be completed in a more comprehensive future study. The ice and piezoelectric actuator patches were attached to the blade shell structure using Tie constraints.