Ultrasonic Deicing E ﬃ ciency Prediction and Validation for a Flat Deicing System

: An e ﬀ ective deicing system is needed to be designed to conveniently remove ice from the surfaces of structures. In this paper, an ultrasonic deicing system for di ﬀ erent conﬁgurations was estimated and veriﬁed based on ﬁnite element simulations. The research focused on deicing e ﬃ ciency factor (DEF) discussions, prediction, and validations. Firstly, seven di ﬀ erent conﬁgurations of Lead zirconate titanate (PZT) disk actuators with the same volume but di ﬀ erent radius and thickness were adopted to conduct harmonic analysis. The e ﬀ ects of PZT shape on shear stresses and optimal frequencies were obtained. Simultaneously, the average shear stresses at the ice / substrate interface and total energy density needed for deicing were calculated. Then, a coe ﬃ cient named deicing e ﬃ ciency factor (DEF) was proposed to estimate deicing e ﬃ ciency. Based on these results, the optimized conﬁguration and deicing frequency are given. Furthermore, four di ﬀ erent icing cases for the optimize conﬁguration were studied to further verify the rationality of DEF. The e ﬀ ects of shear stress distributions on deicing e ﬃ ciency were also analyzed. At same time, a cohesive zone model (CZM) was introduced to describe interface behavior of the plate and ice layer. Standard-explicit co-simulation was utilized to model the wave propagation and ice layer delamination process. Finally, the deicing experiments were carried out to validate the feasibility and correctness of the deicing system.


Introduction
Icing on the equipment is one of the problems needed to be solved in a cold region. Ice buildup on the structures may cause many adverse effects, including cracking, stiffness degradation, malfunctions, and low performance on the structural components, which remarkably endangers their reliability. Icing problems are widely existing in aircraft [1,2], wind turbine blades [3,4], marine structures [5,6] etc. So far, various deicing methods were proposed and developed, mainly including hot bleed air, electro-thermal, electro-impulse, electro-expulsive, vibratory, high frequency microwave, deicing fluids etc. [7][8][9]. Based on the theory of wave propagation in a multilayered structure, Lamb waves and shear horizontal (SH) waves can generate stress and displacement at the interface of the structure. These stresses can be utilized to detach the ice layer. Compared to the deicing techniques mentioned above, ultrasonic deicing is more advantageous in terms of power consumption and efficiency. Meanwhile, the damping coefficient of the Lamb wave is small and this property guarantees the long propagation distance of the Lamb wave. Thus, the Lamb wave has an advantage in deicing for large structures. system and ultrasonic waves propagating in the ice-plate layered structure for deicing is shown in Figure 1.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 21 [27]. The coordinate system and ultrasonic waves propagating in the ice-plate layered structure for deicing is shown in Figure 1. The governing equation for Lamb waves and SH waves can be written as: where ρ is the density, Cijkl is the stiffness tensor, and u is the displacement field. The waves are assumed to propagate along x1 direction. The solution of Equation (1) can be assumed as the following form: ( In Equation (2), Ui = Uαi, αi is the cosine component in i direction. Ui is the polarization vector presenting the displacement in each direction, k stands for the wave number along x1 direction, c represents the phase velocity along x1 direction, and α is the ratio of wave number in x3 with respect to the wave number in the x1 direction. By utilizing the inverse method [28], Christoffel's equation can be written as Equation (3):   2  11  12  13  1  2  12  22  23  2  2  13 23 33 3 0 0 0 The values of λim (i, m = 1, 2, 3) in Equation (3) can be calculated from Equation (4): where nk, nl are the direction cosines of the normal to the wave front. According to the boundary conditions, the displacements u1 and u3 are zero for SH waves and u2 is zero for Lamb waves. Thus, displacements and stresses can be expressed as Equation (5): The governing equation for Lamb waves and SH waves can be written as: where ρ is the density, C ijkl is the stiffness tensor, and u i is the displacement field. The waves are assumed to propagate along x 1 direction. The solution of Equation (1) can be assumed as the following form: In Equation (2), U i = Uα i , α i is the cosine component in i direction. U i is the polarization vector presenting the displacement in each direction, k stands for the wave number along x 1 direction, c represents the phase velocity along x 1 direction, and α is the ratio of wave number in x 3 with respect to the wave number in the x 1 direction. By utilizing the inverse method [28], Christoffel's equation can be written as Equation (3): The values of λ im (i, m = 1, 2, 3) in Equation (3) can be calculated from Equation (4): where n k , n l are the direction cosines of the normal to the wave front. According to the boundary conditions, the displacements u 1 and u 3 are zero for SH waves and u 2 is zero for Lamb waves. Thus, displacements and stresses can be expressed as Equation (5): B k e ik(x 1 +α k x 3 −ct) B k e ik(x 1 +α k x 3 −ct) B k U 3k e ik(x 1 +α k x 3 −ct) Appl. Sci. 2020, 10, 6640 4 of 21 B k α k µ(ik)e ik(x 1 +α k x 3 −ct) where B k is the weighting coefficient of the partial waves. U 3k is the relation of the polarization vectors between U 3 and U 1 . The other variables are the same as mentioned before.

FE Model and Simulation
In this section, harmonic analysis and dynamic analysis were performed to model deicing effect and deicing process, respectively. Harmonic analysis was implemented to calculate the stress tensor at expected frequency range. The finite element model equation of the PZT actuator bonded to the host structure can be written as Equation (6) [29]: where F is applied mechanical force, I is excitation current, U is elastic displacement field, Φ is electric potential, [K nn ] is elastic rigidity matrix, [K nΦ ] is piezoelectric stiffness matrix, [K ΦΦ ] is permittivity matrix, [M] is mass matrix, [R] is dissipation matrix, and w is the angular frequency. The equation of harmonic analysis used in FE software can be expressed as follows: The work started from the evaluation of the effectiveness of guided waves technology to produce the ice detaching from the structure. The 3D model of finite element simulation is shown in Figure 2.
where Bk is the weighting coefficient of the partial waves. U3 is the relation of the polarization vectors between U3 and U1. The other variables are the same as mentioned before.

FE Model and Simulation
In this section, harmonic analysis and dynamic analysis were performed to model deicing effect and deicing process, respectively. Harmonic analysis was implemented to calculate the stress tensor at expected frequency range. The finite element model equation of the PZT actuator bonded to the host structure can be written as Equation (6) where F is applied mechanical force, I is excitation current, U is elastic displacement field, Φ is electric potential, The work started from the evaluation of the effectiveness of guided waves technology to produce the ice detaching from the structure. The 3D model of finite element simulation is shown in Figure 2. The geometry size of the aluminum plate and ice layer was 200 mm × 150 mm × 2 mm and 140 mm × 80 mm × 3 mm. The material of the piezoelectric disc was selected to be PZT-4. The material properties of the aluminum and ice are shown in Table 1.  The geometry size of the aluminum plate and ice layer was 200 mm × 150 mm × 2 mm and 140 mm × 80 mm × 3 mm. The material of the piezoelectric disc was selected to be PZT-4. The material properties of the aluminum and ice are shown in Table 1. Material properties for PZT-4 are shown as follows: where ρ is density, c E is elasticity matrix, e is piezoelectric coupling constant, and ε is dielectric matrix. The FE simulations were performed with an Abaqus ® during processing. The host plate and ice layer were meshed by using 20-node quadratic brick elements with reduced integration (C3D20R) [30] and PZT actuator was meshed by using 20-node quadratic piezoelectric brick elements with reduced integration (C3D20RE). In the frequency range of research , the wavelengths of the Lamb wave and SH wave were taken into consideration. According to the theory of guided waves in layered structures in Ref. [27], the global matrix method can be used to solve the dispersion curves of the ice-plate layered structure when ice accretes on the plate. By solving dispersion curve equations for both Lamb wave and SH wave, the phase velocity of both Lamb wave and SH wave can be calculated. The results for the 2 mm aluminum plate and 2 mm aluminum plate with a 3 mm ice layer are shown in Figure 3.
where ρ is density, E c is elasticity matrix, e is piezoelectric coupling constant, and ε is dielectric matrix. The FE simulations were performed with an Abaqus ® during processing. The host plate and ice layer were meshed by using 20-node quadratic brick elements with reduced integration (C3D20R) [30] and PZT actuator was meshed by using 20-node quadratic piezoelectric brick elements with reduced integration (C3D20RE). In the frequency range of research (20-240 kHz), the wavelengths of the Lamb wave and SH wave were taken into consideration. According to the theory of guided waves in layered structures in Ref. [27], the global matrix method can be used to solve the dispersion curves of the ice-plate layered structure when ice accretes on the plate. By solving dispersion curve equations for both Lamb wave and SH wave, the phase velocity of both Lamb wave and SH wave can be calculated. The results for the 2 mm aluminum plate and 2 mm aluminum plate with a 3 mm ice layer are shown in Figure 3. According to phase velocity, the shortest wavelengths should be around 8.0 mm for the 2 mm aluminum plate with a 3 mm ice layer and around 7.8 mm for the 2 mm aluminum plate at frequency 240 kHz. After several times trial calculations, the global mesh seeding size was chosen as 0.74 mm, smaller than one over ten the smallest wavelength. Meanwhile, the mesh type adapted in simulations was 20-node quadratic brick elements with reduced integration, rather 8-node linear brick elements According to phase velocity, the shortest wavelengths should be around 8.0 mm for the 2 mm aluminum plate with a 3 mm ice layer and around 7.8 mm for the 2 mm aluminum plate at frequency 240 kHz. After several times trial calculations, the global mesh seeding size was chosen as 0.74 mm, smaller than one over ten the smallest wavelength. Meanwhile, the mesh type adapted in simulations was 20-node quadratic brick elements with reduced integration, rather 8-node linear brick elements with reduced integration to ensure accuracy. Thus, the amount of calculation would be several times larger than using mesh type 8-node linear brick elements with reduced integration. After meshing the model, the number of elements was more than 310,000. The FE mesh and reference points (A and B) for FE study is shown in Figure 4. In the model, the voltage applied to the positive and negative electrode of the actuator was 25 and 0 V, respectively. The left and right edge of the aluminum plate were fixed. To make out the effect of different actuators on deicing and stress distributions, the volume of the PZT actuator was kept as 100 π mm 3 , while their radius (r) and thickness (t) were changing as seven different configurations. The radius for seven configurations were 4, 5, 6, 7, 8, 9, and 10 mm, respectively.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 21 larger than using mesh type 8-node linear brick elements with reduced integration. After meshing the model, the number of elements was more than 310,000. The FE mesh and reference points (A and B) for FE study is shown in Figure 4. In the model, the voltage applied to the positive and negative electrode of the actuator was 25 and 0 V, respectively. The left and right edge of the aluminum plate were fixed. To make out the effect of different actuators on deicing and stress distributions, the volume of the PZT actuator was kept as 100 π mm 3 , while their radius (r) and thickness (t) were changing as seven different configurations. The radius for seven configurations were 4, 5, 6, 7, 8, 9, and 10 mm, respectively. Then, the deicing efficiency factor (DEF) was defined. The deicing efficiency for these seven configurations were compared. Based on simulation results in Section 4, the most efficient configuration was confirmed to be configuration 2. To verify the rationality of DEF, four different icing cases as shown in Figure 5 were studied. For these four cases, the radius of the PZT actuator was 5 mm. The icing case (a) in Figure 5 was the same as configuration 2 mentioned above. The size of the icing area and the centroid of ice layer for these four cases were the same. For the cases (b-d) in Figure 5, the distance between each ice block was 4 mm. Meanwhile, the applied voltage and frequency were the same for these four cases. Based on these, shear stresses and DEF were computed and verified. Finally, Abaqus standard-explicit co-simulation was conducted to model wave propagation and the deicing process for configuration 2. The aluminum plate and ice layer were defined as the parent model to conduct transient dynamic analysis in Abaqus explicit, while the PZT actuator was defined as the child model to conduct Abaqus implicit analysis. An interaction interface between the PZT Then, the deicing efficiency factor (DEF) was defined. The deicing efficiency for these seven configurations were compared. Based on simulation results in Section 4, the most efficient configuration was confirmed to be configuration 2. To verify the rationality of DEF, four different icing cases as shown in Figure 5 were studied. For these four cases, the radius of the PZT actuator was 5 mm. The icing case (a) in Figure 5 was the same as configuration 2 mentioned above. The size of the icing area and the centroid of ice layer for these four cases were the same. For the cases (b-d) in Figure 5, the distance between each ice block was 4 mm. Meanwhile, the applied voltage and frequency were the same for these four cases. Based on these, shear stresses and DEF were computed and verified.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 20 were fixed. To make out the effect of different actuators on deicing and stress distributions, the volume of the PZT actuator was kept as 100 π mm 3 , while their radius (r) and thickness (t) were changing as seven different configurations. The radius for seven configurations were 4, 5, 6, 7, 8, 9, and 10 mm, respectively. Then, the deicing efficiency factor (DEF) was defined. The deicing efficiency for these seven configurations were compared. Based on simulation results in Section 4, the most efficient configuration was confirmed to be configuration 2. To verify the rationality of DEF, four different icing cases as shown in Figure 5 were studied. For these four cases, the radius of the PZT actuator was 5 mm. The icing case (a) in Figure 5 was the same as configuration 2 mentioned above. The size of the icing area and the centroid of ice layer for these four cases were the same. For the cases (b)-(d) in Figure 5, the distance between each ice block was 4 mm. Meanwhile, the applied voltage and frequency were the same for these four cases. Based on these, shear stresses and DEF were computed and verified. Finally, Abaqus standard-explicit co-simulation was conducted to model wave propagation and the deicing process for configuration 2. The aluminum plate and ice layer were defined as the parent model to conduct transient dynamic analysis in Abaqus explicit, while the PZT actuator was defined as the child model to conduct Abaqus implicit analysis. An interaction interface between the PZT Finally, Abaqus standard-explicit co-simulation was conducted to model wave propagation and the deicing process for configuration 2. The aluminum plate and ice layer were defined as the parent model to conduct transient dynamic analysis in Abaqus explicit, while the PZT actuator was defined as the child model to conduct Abaqus implicit analysis. An interaction interface between the PZT actuator and aluminum plate was used to manage the data exchange between the two solvers. A bilinear cohesive zone model (CZM) [23] was considered to simulate delamination behavior at the interface of the ice layer and aluminum plate. Cohesive elements with zero thickness were inserted between the aluminum plate and ice layer. The constitutive behavior of the cohesive layer was defined as a traction-separation law to model the detachment at the interface of the ice layer and aluminum plate. It is assumed that there is a linear elastic traction separation law before failure, and the material stiffness will gradually degrade after failure. To describe cohesive material behavior, damage initiation and damage evolution should be defined. Initiation of interface damage is based on quadratic traction interaction failure criterion, written as Equation (8): where σ z , τ xz , τ yz represent the normal stress, shear stress in the x and y direction, separately. T and S are the tensile strength and shear strength, respectively. The damage evolution satisfies the mixed-mode fracture energy approach [31,32]. The fracture energy for the aluminum/ice interface was reported to be 1.0 N/m [33]. Therefore, the normal mode fracture energy G IC was considered to be 1.0 × 10 −3 N/mm in this study. The shear mode fracture energy G IIC and G IIIC was 2.0 × 10 −3 N/mm. The failure criterion can be expressed as Equation (9): where G I , G II , G III are energy produced in the normal, the x and y direction, respectively. G IC , G II C , G IIIC are the critical fracture energy to achieve debond failure in the normal, the x, and the y direction, respectively. The value of exponent α is usually 1 or 2. It was selected to be 2 in this study. The relevant parameters can refer to literatures reported before [31]. All parameters are summarized in Table 2.
To obtain the convergent result of the explicit dynamic procedure, the element size and time step should be considered. The element size was the same as mention above (0.74 mm). According to Courant criteria and Nyquist criteria, the time step was sited to be 2 × 10 −8 s. Time period was 6 × 10 −4 s. The voltage applied the actuator was a sinusoidal load of 25 V. The other boundary condition was the same as mentioned in the harmonic analysis.
In order to solve such a huge amount of calculation, the supercomputer of National Supercomputer Center in Guangzhou, which is currently the fourth fastest supercomputer in the world, with a measured 33.86 petaflop/s (quadrillions of calculations per second), was used to perform these simulations. In these processes, we changed the number of nodes which is related to the number of CPUs used for computation, threads, and processes by running scripts to adjust the efficiency of parallel operation of supercomputer Tianhe-2. In this way, simulation results can be obtained most efficiently.

Deicing Experimental Setup and Schemes
On the basis of FE simulations, the experimental platform was built. In this part, the most efficient configuration 2 (according to simulation results in part 4) was mainly considered to validate the feasibility and rationality of the ultrasonic deicing system.
Ice adhesion strength varies from icing temperature, humid, surface roughness, measurement error etc. In this study, the ice adhesion strength was measured using a low temperature centrifuge, which is similar to these reported [16,34]. After completing all settings, the temperature of the low temperature centrifuge was turned to −12 • C for another 40 min and the centrifuge apparatus started. The rotate speed of the centrifuge was recorded when the ice layer was detached from the sample. In the calculation process, the centrifugal force was written as F = mrω 2 . Then, ice adhesion strength can be calculated based on equation τ = F/A. F is centrifugal force (N), m is ice mass (0.00862 ± 0.0002 kg), r is radius of centrifuge apparatus (0.135 m), ω is speed of rotation (rad/s), τ is ice adhesion stress (Pa), A is icing area (0.00070686 m 2 ). Ice adhesion strength was calculated to be 0.36 ± 0.06 MPa after 8 experiments.
All the settings were the same as mentioned in Section 2. In the experiment, a piezoelectric disk with radius 5 mm and thickness 4 mm was used as an actuator to generate ultrasonic waves. The piezoelectric disk was bonded to the aluminum plate by epoxy resin. The setup for the deicing experiments is shown in Figure 6. To eliminate the effect of surrounding temperature on the results, the whole experiment process was conducted in the refrigerator at −12 • C.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 21 efficiency of parallel operation of supercomputer Tianhe-2. In this way, simulation results can be obtained most efficiently.

Deicing Experimental Setup and Schemes
On the basis of FE simulations, the experimental platform was built. In this part, the most efficient configuration 2 (according to simulation results in part 4) was mainly considered to validate the feasibility and rationality of the ultrasonic deicing system.
Ice adhesion strength varies from icing temperature, humid, surface roughness, measurement error etc. In this study, the ice adhesion strength was measured using a low temperature centrifuge, which is similar to these reported [16,34]. After completing all settings, the temperature of the low temperature centrifuge was turned to −12 °C for another 40 min and the centrifuge apparatus started. The rotate speed of the centrifuge was recorded when the ice layer was detached from the sample. In the calculation process, the centrifugal force was written as F = mrω 2 . Then, ice adhesion strength can be calculated based on equation τ = F/A. F is centrifugal force (N), m is ice mass (0.00862 ± 0.0002 kg), r is radius of centrifuge apparatus (0.135 m), ω is speed of rotation (rad/s), τ is ice adhesion stress (Pa), A is icing area (0.00070686 m 2 ). Ice adhesion strength was calculated to be 0.36 ± 0.06 MPa after 8 experiments.
All the settings were the same as mentioned in Section 2. In the experiment, a piezoelectric disk with radius 5 mm and thickness 4 mm was used as an actuator to generate ultrasonic waves. The piezoelectric disk was bonded to the aluminum plate by epoxy resin. The setup for the deicing experiments is shown in Figure 6. To eliminate the effect of surrounding temperature on the results, the whole experiment process was conducted in the refrigerator at −12 °C. Three different frequencies 190, 213, and 236 kHz, respectively, were adopted to verify optimal frequency for configuration 2. These three frequencies were optimal frequencies for configuration 7 (r = 10 mm), configuration 2 (r = 5 mm), and configuration 3 (r = 6 mm), respectively. All the boundary conditions were the same for these three deicing experiments (including the size and location of ice layer), except the frequency applied to the PZT actuator. After 90 s testing for these three different frequencies, the final deicing effects can be obtained. Moreover, deicing experiments were also carried out for the remaining six groups to compare their differences of deicing effect and validate the correctness of simulation results. Finally, deicing experiments of four different icing cases for configuration 2 were carried out to verify the rationality of DEF. Three different frequencies 190, 213, and 236 kHz, respectively, were adopted to verify optimal frequency for configuration 2. These three frequencies were optimal frequencies for configuration 7 (r = 10 mm), configuration 2 (r = 5 mm), and configuration 3 (r = 6 mm), respectively. All the boundary conditions were the same for these three deicing experiments (including the size and location of ice layer), except the frequency applied to the PZT actuator. After 90 s testing for these three different frequencies, the final deicing effects can be obtained. Moreover, deicing experiments were also carried out for the remaining six groups to compare their differences of deicing effect and validate the correctness of simulation results. Finally, deicing experiments of four different icing cases for configuration 2 were carried out to verify the rationality of DEF.

Deicing Energy Efficiency
After harmonic analysis, several nodes were selected as reference points to confirm optimal frequencies for each configuration. Due to the amount of data, only configuration 1 (r = 4 mm, t = 6.25 mm) was shown as an example to interpret the whole process. Here, only the two reference points A and B shown in Figure 4, which are located at the boundary and center of the interface between the host plate and ice layer, respectively, were presented as an example to output shear stress amplitude-frequency curves. The shear stress amplitude-frequency curves are shown in Figure 7. The optimal frequency results for the rest nodes were the same, so they were not presented here.

Deicing Energy Efficiency
After harmonic analysis, several nodes were selected as reference points to confirm optimal frequencies for each configuration. Due to the amount of data, only configuration 1 (r = 4 mm, t = 6.25 mm) was shown as an example to interpret the whole process. Here, only the two reference points A and B shown in Figure 4, which are located at the boundary and center of the interface between the host plate and ice layer, respectively, were presented as an example to output shear stress amplitudefrequency curves. The shear stress amplitude-frequency curves are shown in Figure 7. The optimal frequency results for the rest nodes were the same, so they were not presented here. From the results of Figure 7, it can be found that the optimal frequency for deicing is 221 kHz. Then, the same processes were performed for the remaining six configurations to obtain the optimal frequency of each model. For the purpose of estimating the deicing effect of each model, the average From the results of Figure 7, it can be found that the optimal frequency for deicing is 221 kHz. Then, the same processes were performed for the remaining six configurations to obtain the optimal frequency of each model. For the purpose of estimating the deicing effect of each model, the average shear stresses for delamination τ 13 and τ 23 at the interface of the ice layer were calculated. The results are shown in Table 3. From the results of Table 3, the total shear stress τ total for delamination can be calculated from Formula (10): The total shear stresses for seven different configurations are shown in Figure 8. These shear stresses were enough to overcome ice adhesion strength measured in lab condition, which was 0.36 ± 0.06 MPa at temperature −12 • C (detailed information is in part 3).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 21 shear stresses for delamination τ13 and τ23 at the interface of the ice layer were calculated. The results are shown in Table 3.
The total shear stresses for seven different configurations are shown in Figure 8. These shear stresses were enough to overcome ice adhesion strength measured in lab condition, which was 0.36 ± 0.06 MPa at temperature −12 °C (detailed information is in part 3). Besides deicing shear stresses analysis, energy consumed for the deicing system is another important aspect needed to take into consideration. For simplicity, the energy dissipation and heat loss were ignored, only the kinetic energy and strain energy for the whole system were considered to estimate the energy consumed for each model at optimal frequency. Therefore, the energy density for deicing can be written as Formula (11): where w, Ek, U, and V represent energy density, kinetic energy, strain energy, and volume of system, respectively, for deicing. From these data, the energy density needed to detach ice from the host plate for each model can be estimated. From computation results, the conclusion can be drawn that the energy density needed for deicing varies from the radius of the piezoelectric wafer with the system at optimal frequency for each configuration. Meanwhile, the deicing effect also differs from one to Besides deicing shear stresses analysis, energy consumed for the deicing system is another important aspect needed to take into consideration. For simplicity, the energy dissipation and heat loss were ignored, only the kinetic energy and strain energy for the whole system were considered to estimate the energy consumed for each model at optimal frequency. Therefore, the energy density for deicing can be written as Formula (11): where w, E k , U, and V represent energy density, kinetic energy, strain energy, and volume of system, respectively, for deicing. From these data, the energy density needed to detach ice from the host plate for each model can be estimated. From computation results, the conclusion can be drawn that the energy density needed for deicing varies from the radius of the piezoelectric wafer with the system at optimal frequency for each configuration. Meanwhile, the deicing effect also differs from one to another. Therefore, it is difficult to distinguish which model is better. The ideal goal is to achieve deicing with minimum energy consumption. Therefore, a coefficient called deicing efficiency factor (DEF) was proposed to describe the relationship between deicing effect and deicing energy density. DEF was defined as Equation (12): where τ total , τ 13 , τ 23 , w, E k , U, and V have the same meaning as mentioned prior. This coefficient can describe deicing efficiency for each model. The greater this coefficient is, the greater total shear stress obtained at the same energy density, indicating that more energy is applied to detach ice, which represents the system as more efficient. From above data and Equation (12), DEF for each configuration can be calculated. The energy density needed for deicing and DEF of seven different configuration are shown in Figure 9. From these simulation results, it can be found that configuration 2 (r = 5 mm, t = 4 mm) is the most efficient one.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 21 another. Therefore, it is difficult to distinguish which model is better. The ideal goal is to achieve deicing with minimum energy consumption. Therefore, a coefficient called deicing efficiency factor (DEF) was proposed to describe the relationship between deicing effect and deicing energy density. DEF was defined as Equation (12): where τtotal, τ13, τ23, w, Ek, U, and V have the same meaning as mentioned prior. This coefficient can describe deicing efficiency for each model. The greater this coefficient is, the greater total shear stress obtained at the same energy density, indicating that more energy is applied to detach ice, which represents the system as more efficient. From above data and Equation (12), DEF for each configuration can be calculated. The energy density needed for deicing and DEF of seven different configuration are shown in Figure 9. From these simulation results, it can be found that configuration 2 (r = 5 mm, t = 4 mm) is the most efficient one. Then four icing cases based on configuration 2 mentioned in Section 2.2 were studied and compared. The shear stresses and DEF for these four cases were computed and shown in Figure 10. For these four icing cases, the deicing energy density was the same for the voltage and frequency applied to the PZT actuator. Therefore, the trend of DEF is consistent with that of the total shear stress for delamination (as shown in Figure 10b). Based on the simulation results, deicing experiments can be carried out to verify the rationality of DEF.  Then four icing cases based on configuration 2 mentioned in Section 2.2 were studied and compared. The shear stresses and DEF for these four cases were computed and shown in Figure 10. For these four icing cases, the deicing energy density was the same for the voltage and frequency applied to the PZT actuator. Therefore, the trend of DEF is consistent with that of the total shear stress for delamination (as shown in Figure 10b). Based on the simulation results, deicing experiments can be carried out to verify the rationality of DEF.

Stress Distribution and Deicing Affect Analysis
The following study revolved around configuration 2. Stress distribution of ice layer was mainly investigated in this section. The ice/substrate interface shear stress nephogram at XZ plane and YZ plane at frequency 213 kHz were studied. In order to show these results more intuitively, 13 paths at the interface of the ice layer including left and right boundaries were established along the long side of the ice layer, which were parallel to the Y axis. Similarly, 22 paths at the interface of the ice layer including two boundaries were established along the short side of ice layer, which were parallel to the X axis. Based on these, the variations of shear stress with paths can be obtained. Then, the shear stress amplitude of 13 paths parallel to the Y axis and 22 paths parallel to the X axis were averaged separately. Average the shear stress amplitude of 13 paths parallel to the Y axis, the average shear stress τ 13 and τ 23 vary along the Y axis can be calculated. Similarly, the results of average shear stress amplitude vary along the X axis can also be estimated. The results are shown in Figure 11. Then four icing cases based on configuration 2 mentioned in Section 2.2 were studied and compared. The shear stresses and DEF for these four cases were computed and shown in Figure 10. For these four icing cases, the deicing energy density was the same for the voltage and frequency applied to the PZT actuator. Therefore, the trend of DEF is consistent with that of the total shear stress for delamination (as shown in Figure 10b). Based on the simulation results, deicing experiments can be carried out to verify the rationality of DEF.

Stress Distribution and Deicing Affect Analysis
The following study revolved around configuration 2. Stress distribution of ice layer was mainly investigated in this section. The ice/substrate interface shear stress nephogram at XZ plane and YZ plane at frequency 213 kHz were studied. In order to show these results more intuitively, 13 paths at the interface of the ice layer including left and right boundaries were established along the long side of the ice layer, which were parallel to the Y axis. Similarly, 22 paths at the interface of the ice layer including two boundaries were established along the short side of ice layer, which were parallel to the X axis. Based on these, the variations of shear stress with paths can be obtained. Then, the shear stress amplitude of 13 paths parallel to the Y axis and 22 paths parallel to the X axis were averaged separately. Average the shear stress amplitude of 13 paths parallel to the Y axis, the average shear stress τ13 and τ23 vary along the Y axis can be calculated. Similarly, the results of average shear stress amplitude vary along the X axis can also be estimated. The results are shown in Figure 11. From the simulation results of Figure 11a,b, the conclusion can be drawn that the average shear stress τ13 and τ23 amplitudes are symmetric about the neutral axis. Combining these two figures results with shear stress nephogram, it can be found that τ13 is symmetric about the neutral axis, while τ23 is antisymmetric about the neutral axis. From the results of Figure 11, it can be seen that the average shear stresses are higher at edges and corners of the ice layer which means ice debonding starts from the edges and corners of the ice layer. Meanwhile, the shear stress distributions of deicing efficiency top three, configuration 2 (r = 5 mm), configuration 6 (r = 9 mm), and configuration 7 (r = 10 mm) were compared. In order to find the principle of stress distributions of these three configurations, the shear stress τ13 and τ23 were normalized. The results are shown in Figures 12 and From the simulation results of Figure 11a,b, the conclusion can be drawn that the average shear stress τ 13 and τ 23 amplitudes are symmetric about the neutral axis. Combining these two figures results with shear stress nephogram, it can be found that τ 13 is symmetric about the neutral axis, while τ 23 is antisymmetric about the neutral axis. From the results of Figure 11, it can be seen that the average shear stresses are higher at edges and corners of the ice layer which means ice debonding starts from the edges and corners of the ice layer. Meanwhile, the shear stress distributions of deicing efficiency top three, configuration 2 (r = 5 mm), configuration 6 (r = 9 mm), and configuration 7 (r = 10 mm) were compared. In order to find the principle of stress distributions of these three configurations, the shear stress τ 13 and τ 23 were normalized. The results are shown in Figures 12 and 13.   Figures 12 and 13 reveal the relationship between shear stress distribution and deicing efficiency. The normalized shear stress distribution can reflect the stress gradient from edges to the center of the ice layer. In order to describe this type of stress gradient, a coefficient by referring to the definition of variance was defined. The coefficient named variance of maximum shear stress (VMSS) represents the degree of dispersion between shear stress at the edges and middle of the ice layer. For the average shear stress at the edges of the ice layer at maximum, the average shear stress at the edges would be equal to one after normalization. Thus, this coefficient can be described as Formula (13): where N, τi represent the number of node and shear stress for each node, respectively, τmax means shears stress at edges, which equals to one after normalization. The larger the value of VMSS, the larger the degree of dispersion. That presents the shear stresses at the edges of the ice layer reach the deicing stress threshold, the shear stresses at the rest of the part of the ice layer are still relatively small. It means shear stresses are more concentrated on the edges, which is more efficiency and power saving. The results for these three configurations are shown in Figure 14.   Figures 12 and 13 reveal the relationship between shear stress distribution and deicing efficiency. The normalized shear stress distribution can reflect the stress gradient from edges to the center of the ice layer. In order to describe this type of stress gradient, a coefficient by referring to the definition of variance was defined. The coefficient named variance of maximum shear stress (VMSS) represents the degree of dispersion between shear stress at the edges and middle of the ice layer. For the average shear stress at the edges of the ice layer at maximum, the average shear stress at the edges would be equal to one after normalization. Thus, this coefficient can be described as Formula (13): where N, τi represent the number of node and shear stress for each node, respectively, τmax means shears stress at edges, which equals to one after normalization. The larger the value of VMSS, the larger the degree of dispersion. That presents the shear stresses at the edges of the ice layer reach the deicing stress threshold, the shear stresses at the rest of the part of the ice layer are still relatively small. It means shear stresses are more concentrated on the edges, which is more efficiency and power saving. The results for these three configurations are shown in Figure 14.  Figures 12 and 13 reveal the relationship between shear stress distribution and deicing efficiency. The normalized shear stress distribution can reflect the stress gradient from edges to the center of the ice layer. In order to describe this type of stress gradient, a coefficient by referring to the definition of variance was defined. The coefficient named variance of maximum shear stress (VMSS) represents the degree of dispersion between shear stress at the edges and middle of the ice layer. For the average shear stress at the edges of the ice layer at maximum, the average shear stress at the edges would be equal to one after normalization. Thus, this coefficient can be described as Formula (13): where N, τ i represent the number of node and shear stress for each node, respectively, τ max means shears stress at edges, which equals to one after normalization. The larger the value of VMSS, the larger the degree of dispersion. That presents the shear stresses at the edges of the ice layer reach the deicing stress threshold, the shear stresses at the rest of the part of the ice layer are still relatively small.
It means shear stresses are more concentrated on the edges, which is more efficiency and power saving. The results for these three configurations are shown in Figure 14.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 21 Figure 14. Variance of maximum shear stress along Y axis (a) and X axis (b).
In Figure 14, coefficient VMSS shows a downward trend from configuration 2 (r = 5 mm), configuration 6 (r = 9 mm) to configuration 7 (r = 10 mm). This conclusion is consistent with the results of coefficient DEF, which proves the correctness of DEF from another perspective. In addition, 15 paths along the thickness direction of the ice layer were established to analyze the shear stress distribution along the Z axis direction. The results are shown in Figure 15. The result of Figure 15 shows that amplitude of shear stress τ13 and τ23 decrease from interface to free surface while shear stress τ12 exists along the entire thickness direction of ice layer, which is corresponding with theory and analysis.

Standard-Explicit Co-Simulation Analysis
With the help of transient dynamic analysis by using standard-explicit co-simulation, the whole process of wave propagation and deicing can be revealed more clearly. Similarly, configuration 2 was studied in the process. The distribution of displacement for the aluminum plate and ice layer at different times was calculated and shown in Figure 16. After the PZT actuator was turned on 6 μs, waves were generated and started to spread around the actuator. After 12 μs, waves reached and reflected from the left edge of the aluminum plate. Then, the reflected waves were superimposed with waves generated by the actuator and propagated toward the ice layer. Around 27 μs later, waves reached the left edge of the ice layer. In the process, the waves continued spreading to the right edge of the aluminum plate and reflected from the edges around 180 μs. Gradually, displacement started In Figure 14, coefficient VMSS shows a downward trend from configuration 2 (r = 5 mm), configuration 6 (r = 9 mm) to configuration 7 (r = 10 mm). This conclusion is consistent with the results of coefficient DEF, which proves the correctness of DEF from another perspective. In addition, 15 paths along the thickness direction of the ice layer were established to analyze the shear stress distribution along the Z axis direction. The results are shown in Figure 15. In Figure 14, coefficient VMSS shows a downward trend from configuration 2 (r = 5 mm), configuration 6 (r = 9 mm) to configuration 7 (r = 10 mm). This conclusion is consistent with the results of coefficient DEF, which proves the correctness of DEF from another perspective. In addition, 15 paths along the thickness direction of the ice layer were established to analyze the shear stress distribution along the Z axis direction. The results are shown in Figure 15. The result of Figure 15 shows that amplitude of shear stress τ13 and τ23 decrease from interface to free surface while shear stress τ12 exists along the entire thickness direction of ice layer, which is corresponding with theory and analysis.

Standard-Explicit Co-Simulation Analysis
With the help of transient dynamic analysis by using standard-explicit co-simulation, the whole process of wave propagation and deicing can be revealed more clearly. Similarly, configuration 2 was studied in the process. The distribution of displacement for the aluminum plate and ice layer at different times was calculated and shown in Figure 16. After the PZT actuator was turned on 6 μs, waves were generated and started to spread around the actuator. After 12 μs, waves reached and reflected from the left edge of the aluminum plate. Then, the reflected waves were superimposed with waves generated by the actuator and propagated toward the ice layer. Around 27 μs later, waves reached the left edge of the ice layer. In the process, the waves continued spreading to the right edge of the aluminum plate and reflected from the edges around 180 μs. Gradually, displacement started The result of Figure 15 shows that amplitude of shear stress τ 13 and τ 23 decrease from interface to free surface while shear stress τ 12 exists along the entire thickness direction of ice layer, which is corresponding with theory and analysis.

Standard-Explicit Co-Simulation Analysis
With the help of transient dynamic analysis by using standard-explicit co-simulation, the whole process of wave propagation and deicing can be revealed more clearly. Similarly, configuration 2 was studied in the process. The distribution of displacement for the aluminum plate and ice layer at different times was calculated and shown in Figure 16. After the PZT actuator was turned on 6 µs, waves were generated and started to spread around the actuator. After 12 µs, waves reached and reflected from the left edge of the aluminum plate. Then, the reflected waves were superimposed with waves generated by the actuator and propagated toward the ice layer. Around 27 µs later, waves reached the left edge of the ice layer. In the process, the waves continued spreading to the right edge of the aluminum plate and reflected from the edges around 180 µs. Gradually, displacement started to concentrate on the ice layer after 429 µs. Then, 498 µs later, the left edge of ice layer started to peel off symmetrically along the neutral axis of the aluminum plate. Delamination developed towards the inside of the ice layer within 519 µs. After that, delamination started at the bottom and top of the ice layer about 546 µs. Then, detachment developed towards the bottom and top of the ice layer after 555 µs. After 564 µs, detachment developed to the center and right edge of the ice layer. Finally, delamination finished at 570 µs.  To reveal the deicing process more clearly, an index for damage degree (SDEG) was utilized to elaborate the interface damage evolution process. Figure 17 shows the evolution process of SDEG for the cohesive layer. When SDEG was equal to 1, the area was detached and deleted. After the actuator was turned on 450 μs, damage initiation occurred at the left, bottom, and top edge of the cohesive layer. After 498 μs, the left edge of the cohesive layer started to damage and deleted. The ice layer stared to detach from the aluminum plate. Then, damage developed along the wave propagation direction after 531 μs. With the evolution of damage, the bottom and top of the ice layer was peeled off at 549 μs. Then, the detached area continued to increase. The detachment developed to the center and right edge of the ice layer around 564 μs. Finally, the delamination finished within 570 μs. The detachment area was around 88.2% of the icing area. Although the ice layer was not totally detached, the remaining area of the cohesive layer also showed varying degrees of damage, which meant that the ice layer was likely to be removed completely under a small external load. To reveal the deicing process more clearly, an index for damage degree (SDEG) was utilized to elaborate the interface damage evolution process. Figure 17 shows the evolution process of SDEG for the cohesive layer. When SDEG was equal to 1, the area was detached and deleted. After the actuator was turned on 450 µs, damage initiation occurred at the left, bottom, and top edge of the cohesive layer. After 498 µs, the left edge of the cohesive layer started to damage and deleted. The ice layer stared to detach from the aluminum plate. Then, damage developed along the wave propagation direction after 531 µs. With the evolution of damage, the bottom and top of the ice layer was peeled off at 549 µs. Then, the detached area continued to increase. The detachment developed to the center and right edge of the ice layer around 564 µs. Finally, the delamination finished within 570 µs. The detachment area was around 88.2% of the icing area. Although the ice layer was not totally detached, the remaining area of the cohesive layer also showed varying degrees of damage, which meant that the ice layer was likely to be removed completely under a small external load. Based on results of dynamic analysis, the processes of wave propagation, reflection, and superposition can be better observed and understood. Besides, the entire peeling process of the ice layer can be studied and presented. From the results of the above simulation, it can be found that the detachment started at the edge of the ice layer, which coincided with the stress distribution analysis in Section 4.2. In the harmonic analysis for configuration 2, the detached area was 83.3%. In the above dynamic analysis, the detached area was 88.2%. It showed that good agreements between these two analyses were obtained.

Experimental Verification
The deicing results for configuration 2 are shown in Figure 18. Three deicing experiments under different frequencies were compared. When frequency was 190 kHz, there was no visible change observed. When frequency was 236 kHz, only some cracks but without any delamination were observed. When frequency was 213 kHz, the ice layer was cracked and partially detached from the aluminum plate. After the actuator was turned on 24 s, delamination started from the edge of the ice layer. The delamination area increased along the wave propagation direction with some cracks after 47 s. Then, the detached area continued to increase with more cracks generated. Finally, delamination finished within 71 s. The detached area could be estimated from the nephogram in an open source image processing software named ImageJ [35]. Firstly, the experimental results image was converted into grayscale images. The detached area can be identified by comparing the images before and after the experiment. Then, the color threshold was set according to the color of the detached area. Based on these, the detached area was calculated. In the ice layer nephogram, the detached area was in white. The percent of shear stress exceeded the adhesion stress at 83.3% and 88.2% in harmonic Based on results of dynamic analysis, the processes of wave propagation, reflection, and superposition can be better observed and understood. Besides, the entire peeling process of the ice layer can be studied and presented. From the results of the above simulation, it can be found that the detachment started at the edge of the ice layer, which coincided with the stress distribution analysis in Section 4.2. In the harmonic analysis for configuration 2, the detached area was 83.3%. In the above dynamic analysis, the detached area was 88.2%. It showed that good agreements between these two analyses were obtained.

Experimental Verification
The deicing results for configuration 2 are shown in Figure 18. Three deicing experiments under different frequencies were compared. When frequency was 190 kHz, there was no visible change observed. When frequency was 236 kHz, only some cracks but without any delamination were observed. When frequency was 213 kHz, the ice layer was cracked and partially detached from the aluminum plate. After the actuator was turned on 24 s, delamination started from the edge of the ice layer. The delamination area increased along the wave propagation direction with some cracks after 47 s. Then, the detached area continued to increase with more cracks generated. Finally, delamination finished within 71 s. The detached area could be estimated from the nephogram in an open source image processing software named ImageJ [35]. Firstly, the experimental results image was converted into grayscale images. The detached area can be identified by comparing the images before and after the experiment. Then, the color threshold was set according to the color of the detached area. Based on these, the detached area was calculated. In the ice layer nephogram, the detached area was in white. The percent of shear stress exceeded the adhesion stress at 83.3% and 88.2% in harmonic analysis and dynamic analysis while it was 76.4% in the experiment. The main error came from the difference in ice properties between the simulation and experiment. Simultaneously, the different performances at frequency 190 and 236 kHz were also corresponding with simulation results. The final results are shown in Table 4. Neither frequency, 190 nor frequency 236 kHz could detach the ice layer, since both shear stresses of them were smaller than the adhesion strength of the ice. The crack existed at frequency 236 kHz rather frequency 190 kHz due to shear stress τ12 which contributed to ice cracking. Shear stress τ12 was higher at frequency 236 than 190 kHz. Thus, the crack could observe at frequency 236 kHz rather than 190 kHz. Simultaneously, the different performances at frequency 190 and 236 kHz were also corresponding with simulation results. The final results are shown in Table 4. Neither frequency, 190 nor frequency 236 kHz could detach the ice layer, since both shear stresses of them were smaller than the adhesion strength of the ice. The crack existed at frequency 236 kHz rather frequency 190 kHz due to shear stress τ 12 which contributed to ice cracking. Shear stress τ 12 was higher at frequency 236 than 190 kHz. Thus, the crack could observe at frequency 236 kHz rather than 190 kHz. Meanwhile, deicing experiments were also performed for the remaining six groups at their optimal frequency, respectively. The trend of detached areas after deicing experiments and total shear stress τ total for delamination from FE simulation results were compared. The trend of total shear stress τ total and the detached area are shown in Figure 19.  Meanwhile, deicing experiments were also performed for the remaining six groups at their optimal frequency, respectively. The trend of detached areas after deicing experiments and total shear stress τtotal for delamination from FE simulation results were compared. The trend of total shear stress τtotal and the detached area are shown in Figure 19. The result of Figure 19 can prove the optimal frequency for configuration 2 was 213 kHz. From the results of Figure 19, it can be found that the trend of detached areas in experiments is consistent with total shear stress in FE simulations, which can verify the correctness of FE simulation results and feasible of the deicing method. Finally, deicing experiments for four different icing cases mentioned before were conducted. The deicing results for these four cases are shown in Figure 20. The trend of the detached area for these four cases is consistent with the trend of DEF predicted in the FE simulation. The result of Figure 19 can prove the optimal frequency for configuration 2 was 213 kHz. From the results of Figure 19, it can be found that the trend of detached areas in experiments is consistent with total shear stress in FE simulations, which can verify the correctness of FE simulation results and feasible of the deicing method. Finally, deicing experiments for four different icing cases mentioned before were conducted. The deicing results for these four cases are shown in Figure 20. The trend of the detached area for these four cases is consistent with the trend of DEF predicted in the FE simulation.
the results of Figure 19, it can be found that the trend of detached areas in experiments is consistent with total shear stress in FE simulations, which can verify the correctness of FE simulation results and feasible of the deicing method. Finally, deicing experiments for four different icing cases mentioned before were conducted. The deicing results for these four cases are shown in Figure 20. The trend of the detached area for these four cases is consistent with the trend of DEF predicted in the FE simulation.

Conclusions
In this paper, a flat deicing system using ultrasonic waves to detach the ice layer from the host structure was studied. The feasibility of the deicing system was validated in both FE simulations and experiments. In FE simulations, both harmonic analysis and dynamic analysis were conducted. In the harmonic analysis, deicing parameters and deicing frequencies were given. For different configuration, a factor (DEF) was proposed to estimate and predict deicing efficiency. Based on optimal configuration, four icing cases were discussed for further validation of DEF. Simultaneously, the relationship between shear stress distributions and deicing was explored. In dynamic analysis, wave propagation and the deicing process were investigated. Finally, deicing experiments were performed to verify the feasibility of the deicing system and the reasonableness of predicted results in FE simulations.
Based on the results above, the following conclusions were drawn:

1.
A parameter design method was proposed according to simulation results on the supercomputer. The optimal configuration (the radius and the thickness of the PZT actuator) and deicing frequency were also given.

2.
The coefficient deicing efficiency factor (DEF) was proposed based on FE simulation, which can estimate and predict the deicing efficiency. With the help of this factor, the effect of the different model and deicing parameters on deicing efficiency can be calculated and compared. For shear stresses distributions, the more concentrated the shear stresses are on the edges of ice/substrate interface, the better for deicing. This can be one of the potential ways to achieve deicing more effectively.

3.
The calculation results of harmonic analysis and dynamic analysis are consistent. Harmonic analysis can obtain deicing stress and deicing results more concise and clearer. Wave propagation, reflection, superposition, and other processes are better displayed in dynamic analysis, which can help to understand the deicing process and provide guidance for further optimization. For different models and problems, these two analyses should be used rationally to achieve the purpose of obtaining the expected results with lower calculation costs.