Numerical Analysis of the Wing Leading Edge Electro-Impulse De-Icing Process Based on Cohesive Zone Model

Abstract

Aircraft icing has historically been a critical cause of airplane crashes. The electro-impulse de-icing system has a wide range of applications in aircraft de-icing due to its lightweight design, low energy consumption, high efficiency, and other advantages. However, there has been little study into accurate wing electric-impulse de-icing simulation methods and the parameters impacting de-icing efficacy. Based on the damage mechanics principle and considering the influence mechanisms of interface debonding and ice fracture on ice shedding, this paper establishes a more accurate numerical model of wing electric-impulse de-icing using the Cohesive Zone Model (CZM). It simulates the process of electric-impulse de-icing at the leading edge of the NACA 0012 wing. The numerical results are compared to the experimental results, revealing that the constructed wing electro-impulse de-icing numerical model is superior. Lastly, the effects of varying ice–skin interface shear adhesion strengths, doubler loading positions, and impulse sequences on de-icing effectiveness were studied. The de-icing rate is a quantitative description of the electro-impulse’s de-icing action, defined in the numerical model as the ratio of cohesive element deletions to the total elements at the ice–skin interface. The findings reveal that varying shear adhesion strengths at the ice–skin interface significantly impact the de-icing effect. The de-icing rate steadily falls with increasing shear adhesion strength, from 66% to 56%. When two, four, and seven impulses were applied to doubler two, the de-icing rates were 59%, 71%, and 71%, respectively, significantly increasing the de-icing efficiency compared to when impulses were applied to doubler one. Doubler one and two impulse responses are overlaid differently depending on the impulse sequences, resulting in varying de-icing rates. When the impulse sequence is 20 ms, the superposition results are optimal, and the de-icing rate reaches 100%. These studies can guide the development and implementation of a wing electric-impulse de-icing system.

1. Introduction

Aircraft icing poses a significant threat to flight safety [1,2]. When the leading edge of an aircraft wing icing, it can increase weight, decrease lift, and degrade aerodynamic performance, ultimately affecting the aircraft’s safety and maneuverability [3]. There are two approaches to addressing aircraft icing: anti-icing and de-icing. Passive anti-icing technology works primarily through surface coating to change surface properties and prevent aircraft surface ice adhesion [4,5,6,7]. However, a single coating technology cannot be used for aircraft de-icing. It must be used in conjunction with active de-icing technology to solve the problem of aircraft icing through a hybrid anti-icing technology [8,9]. As a result, research into active de-icing equipment remains the primary focus of the present aircraft icing problem. Thermal [10,11,12], chemical [13,14], and mechanical de-icing [15,16,17,18,19,20,21,22,23] are examples of active de-icing technologies that can eliminate ice formation before it has visible harmful consequences. Thermal de-icing affects engine efficiency and heat utilization when the engine generates heat [24]. However, it costs much electrical energy and has a restricted application when heat is generated by electric heaters placed under the wing surfaces [25]. Chemical de-icing has been utilized for ground-based anti-icing, although its environmental impact is a source of worry [26,27,28]. One of the oldest mechanical de-icing solutions is airbag de-icing, which uses airbag tubes to expand and protrude outward, shattering the ice covering. This de-icing technology has been employed in aircraft models MV-22 and ATR42; the de-icing effect is poor, and leftover ice increases drag, altering the aircraft’s original aerodynamic shape [29]. Ultrasonic and piezoelectric de-icing are being developed and have not yet been deployed to aircraft [30,31].

Dr. Levine’s publication of the electro-impulse system design sparked a rush of study into electro-impulse de-icing technologies [32]. The basic principle of electro-impulse de-icing technology is that the capacitor, through the switch to the impulse coil discharge, will form a rapid change of the magnetic field; the time-varying magnetic field will be induced in the metal skin of the eddy currents, which produces high amplitude and a short duration of the impulse force; the skin of the vibration in the action of the impulse force, which will enable the ice layer and the skin to be separated; and finally, to achieve the purpose of de-icing. Compared to other de-icing techniques, the electro-impulse de-icing system has been widely used for aircraft de-icing because of its low energy consumption, ease of maintenance, high dependability, and lack of substantial negative impact on engine characteristics. Although the experimental research of the electro-impulse system can view the de-icing process, it requires a particular icing environment, such as a wind tunnel, and the entire experiment must take place in an insulated box. Calculating experimental data for complicated structures can be time-consuming. With the advancement of electrical technology, the application of numerical simulation methods to study electro-impulse deicing technology can avoid these strict research settings, making the calculation of experimental data easier. Furthermore, simulating the ice-shedding process and predicting the de-icing effect is critical for developing electro-impulse de-icing devices. The finite element method can obtain the stress and strain responses of the skin and ice layer during electro-impulse de-icing by discretizing a complex set of differential equations into a linear set of equations. It can calculate the ice-shedding process using damage mechanics principles. The accurate de-icing calculation model enables people to swiftly finish the deicing system design, analysis, and optimization. Khatkhate et al. [33] observed that the shear force between the skin and the ice layer under electro-impulse loading dominated the ice-shedding process and developed a preliminary prediction of the ice-shedding region on the skin surface by comparing the shear stress at the interface between the ice layer and the skin. Zhang et al. [34] suggested a better de-icing criterion that emphasizes the de-icing influence of shear stress, and they employed an explicit central difference algorithm to describe the leading edge structure’s electro-impulse de-icing. These studies highlight ice shedding when the shear stress between the ice and skin reaches a threshold value, but they do not look at how ice shedding is affected by normal skin–ice interaction. Labeas et al. [35] investigated the effect of interface normal tensile stress on ice shedding using shear stress analysis. They presented an ice interface shedding failure criterion based on coupled shear stress and normal tensile stress at the ice interface. Wang et al. [36] also investigated the application of the coupled tensile-shear shedding failure criterion, the single tensile shedding failure criterion, and the single shear shedding failure criterion for electro-impulse de-icing simulations. The tensile shear linked failure criterion is more accurate than the single shedding failure criterion. It has also been argued that one of the reasons for the low accuracy of ice-shedding prediction is that the current model only considers ice interface failure [37].

This paper first introduces the theory of CZM and the Embedding Method of cohesive elements, with an aim to improving the current ice shedding model, which primarily focuses on the effect of ice–skin interface debonding on ice shedding and generally ignores the role of ice fracture on the de-icing effect. Then, it is proven that the cohesive element model can simulate the ice fracture and also characterize the debonding of the ice and the base material using the ice beam three-point bending experiment, single cantilever beam (SCB) experiment, and zero-degree cone angle push-out (push-out) experiment. Ultimately, an accurate finite element model of wing electro-impulse de-icing is established, taking ice fracture and the ice-skin interface debonding into complete account. The modeling’s reasonableness and accuracy are evaluated by comparing skin deflection curves and de-icing rates once the doubler is loaded, both experimentally and numerically. This paper examines the impacts of various ice–skin interface shear adhesion strengths, doubler loading positions, and impulse sequences on the de-icing effect to support the design and implementation of the current wing electric-impulse de-icing system.

2. CZM Theory and Cohesive Element Embedding Methods

2.1. Initial Linear Elastic Traction-Separation Behaviour

CZM, first presented by Dugdale et al. [38] and Barenblatt et al. [39], is a popular method for simulating fracture failure in various materials. For the first time, Hillerborg et al. [40] employed CZM and Finite Element Method (FEM) to simulate fracture in brittle materials. CZM is explained using the traction–separation law. The traction–separation model assumes initial linear elastic behavior, followed by initiation and evolution of damage [41,42]. The elastic behavior is written in an elastic constitutive matrix that relates the normal and shear stresses to the normal and shear separations across the interface. The initial linear elastic behavior of the traction–separation law is represented as follows:

$$t=\left\{\begin{array}{l}{t}_s\\ t_n\\ t_t\end{array}\right\}=\left[\begin{array}{ccc}{E}_{nn}& E_{ns}& E_{nt}\\ E_{ns}& E_{ss}& E_{st}\\ E_{nt}& E_{st}& E_{tt}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_n\\ {\epsilon }_s\\ {\epsilon }_t\end{array}\right\}=\frac{1}{T_0}\left[\begin{array}{ccc}{E}_{nn}& E_{ns}& E_{nt}\\ E_{ns}& E_{ss}& E_{st}\\ E_{nt}& E_{st}& E_{tt}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_n\\ {\epsilon }_s\\ {\epsilon }_t\end{array}\right\}$$

(1)

$$k_{ij}=\frac{E_{ij}}{T_0}$$

(2)

where $t$ is the nominal traction stress vector, which consists of three components in 3D problems ($t_n$, $t_s$, and $t_t$). $E_{ij}$ is Young’s modulus. ${\epsilon }_n$, ${\epsilon }_s$ and ${\epsilon }_t$ are the three components of the nominal strain. ${\delta }_n$, ${\delta }_s$, and ${\delta }_t$ are the corresponding separations. $T_0$ is the original thickness of the cohesive element, and the default value is 1, ensuring that the nominal strain equals the separation displacement. $k_{ij}$ is the material stiffness. The elasticity matrix provides fully coupled behavior between the traction vector and the separation vector components.

2.2. CZM Damage

2.2.1. Damage Initiation

Damage initiation is the initial indication that a material’s stiffness is deteriorating. The degradation process begins when the stress value meets a specific damage beginning criterion, which several factors, including the maximum stress, maximum separation, quadratic stress, and quadratic separation criterion, can determine. The specific characteristics of the material will determine which criterion should be used. However, the quadratic stress criterion is usually applied more often. The following equation demonstrates this criterion.

$${\left\{\frac{\langle t_n\rangle }{t_n^0}\right\}}^2+{\left\{\frac{t_s}{t_s^0}\right\}}^2+{\left\{\frac{t_t}{t_t^0}\right\}}^2=1$$

(3)

In the discussion below, $t_n^0$, $t_s^0$, and $t_t^0$ represent the peak values when the separation is either purely normal ($\mathrm{{\rm I}}$) to the interface or purely in the first ($\mathrm{{\rm I}}\mathrm{{\rm I}}$) or the second ($\mathrm{{\rm I}}\mathrm{{\rm I}}\mathrm{{\rm I}}$) shear direction, respectively. Where $<>$ is the Macaulay bracket, which is adapted to distinguish whether the cohesive element is under compression or tension:

$$\langle t_n\rangle =\{\begin{array}{l}0\Leftarrow t_n\le 0\left(\mathrm{tension}\right)\\ t_n\Leftarrow t_n>0\left(\mathrm{compresion}\right)\end{array}$$

(4)

2.2.2. Damage Evolution

The term “damage evolution” refers to the state of the material after damage has appeared. This stage of the stiffness weakening process can be used to describe the material’s subsequent mechanical properties during the degeneration process. A scalar damage variable, D, represents the overall damage. It initially has a value of 0. If damage evolution is modeled, D monotonically evolves from 0 to 1 upon further loading after damage initiation. The nominal traction stress vector and stiffness are affected by the damage according to

$$k_{ij}=k_{ij}^0(1-D)$$

(5)

$$t_n=\{\begin{array}{l}\left(1-D\right)t_{no},t_{no}\ge 0\\ t_{no},t_{no}<0(\mathrm{no} \mathrm{damage} \mathrm{to} \mathrm{compressive} \mathrm{stiffness})\end{array}$$

(6)

$$t_s=(1-D)t_{so}$$

(7)

$$t_t=(1-D)t_{to}$$

(8)

where $t_{no}$, $t_{so}$, and $t_{to}$ are the contact stress components predicted by the elastic traction–separation behavior for the current separations without damage. $k_{ij}^0$ is the stiffness of the material without damage.

It is evident from the preceding equation that the computation of D holds the key to explaining the damage evolution process of the material. Several intrinsic models are typically used to compute the damage variable D. Many different intrinsic models can often be used to calculate the damage variable D. Giulio et al. [43] compared them and came to the conclusion that the bilinear model could strike a balance between the needs of computing efficiency and accuracy. Thus, in this paper, modeling is completed using the bilinear model. In the model, D is expressed as follows:

$$D=\frac{{\delta }_m^f\left({\delta }_m^{\max }-{\delta }_m^o\right)}{{\delta }_m^{\max }\left({\delta }_m^f-{\delta }_m^o\right)}$$

(9)

where ${\delta }_m^o$ represents the element node’s opening displacement at the beginning of damage (i.e., D > 0), ${\delta }_m^f$ represents the node’s opening displacement at D = 1, and ${\delta }_m^{\max }$ represents the node’s maximum opening displacement during the loading history.

Rather than tensile damage or pure shear in one direction, materials and structures are usually destroyed by a combination of factors. The default CZM in this investigation is the same for $\mathrm{{\rm I}}\mathrm{{\rm I}}$ and $\mathrm{{\rm I}}\mathrm{{\rm I}}\mathrm{{\rm I}}$. According to the cohesive crack model, the crack surface has two tangential traction forces, $t_s$ and $t_t$, and one normal traction force, $t_n$. These traction forces are assumed to decrease monotonically with the corresponding displacement of the crack surface. Figure 1 displays a schematic of the mixed-mode cohesive traction response. Where ${\delta }_n^o$, ${\delta }_s^o$, and ${\delta }_t^o$ are the relative $\mathrm{{\rm I}}$, relative $\mathrm{{\rm I}}\mathrm{{\rm I}}$, and relative $\mathrm{{\rm I}}\mathrm{{\rm I}}\mathrm{{\rm I}}$ at the onset of damage, respectively. As the traction force decreases, ${\delta }_{nf}$, ${\delta }_{sf}$, ${\delta }_{tf}$, and ${\delta }_{mf}$ are the pure $\mathrm{{\rm I}}$ displacement, pure $\mathrm{{\rm I}}\mathrm{{\rm I}}$ displacement, pure $\mathrm{{\rm I}}\mathrm{{\rm I}}\mathrm{{\rm I}}$ displacement, and effective displacement, respectively. The increasing linear upward branch in each softening curve describes the initial intact material. The $\mathrm{{\rm I}}$, $\mathrm{{\rm I}}\mathrm{{\rm I}}$, and $\mathrm{{\rm I}}\mathrm{{\rm I}}\mathrm{{\rm I}}$ fracture energies are represented, respectively, by the areas $G_n^c$, $G_s^c$, $G_t^c$ and bounded by the traction-detachment curves and the axes in Figure 1. The Benzeggagh–Kenane fracture criterion is applied to characterize the outcomes of various damage mechanisms efficiently. This criterion is most useful when the pure $\mathrm{{\rm I}}\mathrm{{\rm I}}$ and $\mathrm{{\rm I}}\mathrm{{\rm I}}\mathrm{{\rm I}}$ fracture energies are equal and is expressed as follows:

$$G_n^c+(G_s^c-G_n^c){\left\{\frac{G_S}{G_T}\right\}}^{\eta }=G^C$$

(10)

Denoting by $G_n$, $G_s$, and $G_t$ the work completed by the tractions and their conjugate separations in the normal, first, and second shear directions, respectively, and defining $G_T=G_n+G_s+G_t$. It is also useful to define the quantity $G_S=G_s+G_t$ to denote the portion of the total work completed by the shear traction and the corresponding separation components. $G^C$ is the effective fracture energy after mixing. $\eta$ is a cohesive property parameter [44].

2.3. Cohesive Element Embedding Methods

The cohesive element methods, which combine the benefits of continuum-based and discontinuity-based techniques, have been developed and applied extensively in recent years to simulate the fracture of brittle and quasi-brittle materials [45,46,47]. The initial finite mesh is embedded with zero-thickness cohesive elements to simulate ice fracture precisely without a predetermined crack path. The main steps are: (a) Set up the finite element model and extract the information on the nodes and solid elements; (b) copy and renumber the specific nodes in the embedded zone according to the number of shared node elements; and (c) set up the zero-thickness cohesive elements according to the coordinate information of the renumbered nodes. A shared face connected with two eight-node linear bricks, reduced integration, and hourglass control (C3D8R) in the initial finite mesh is illustrated in Figure 2a, which is taken as the case for revealing the approach for embedding a 3D cohesive element. Each parent node in the shared interface is replaced with two separate child nodes with the exact coordinates. For example, shared node 1 is changed to node 1 and node $1^{\prime }$. The shared side is replaced by an eight-node cohesive element with zero thickness, as shown in Figure 2b, and the detailed structure of the cohesive element is illustrated in Figure 2c.

3. Validation of Numerical Models

This section begins with numerical simulations of ice fragmentation, followed by type $\mathrm{{\rm I}}$ and type $\mathrm{{\rm I}}\mathrm{{\rm I}}$ numerical simulations of aluminum–ice interface debonding. The validity of the cohesive element characterization of ice fragmentation and aluminum–ice interface debonding is then verified by comparing them to ice beam three-point bending experiments, SCB experiments, and pout-out experiments, respectively. A numerical simulation of wing electric-impulse de-icing is established, and one impulse load is applied to the doubler one without ice conditions. The calculated center deflection curve of doubler one can be compared to the experimental curve to evaluate the effectiveness of the uniform pressure instead of the impulse force. Finally, doubler one received two, four, and seven impulse loads. The accuracy of the numerical model of wing electric-impulse de-icing was confirmed by comparing the simulated and experimental doubler one-center deflection curves and de-icing rates.

3.1. Numerical Simulation of Three-Point Bending of Ice Beams

Based on Zhang et al.’s experiment [48], a numerical simulation of the three-point bending of ice beams is carried out. Before we begin modeling, we must make two assumptions. On the one hand, the fixed support and indenter of the ice beam three-point bending experiment are metal materials with a substantially higher stiffness than ice. Therefore, they are modeled as rigid bodies. This treatment ensures that the fixed support and indenter are not deformed during the simulation, and this assumption has a negligible impact on simulation accuracy. On the other hand, there will be initial defects within the actual ice body, such as microcracks and air bubbles, and the presence of these defects is ignored during the modeling process with the assumption that the ice body has no initial defects. However, studies such as wing de-icing in the later stages of this paper are concerned with the macroscopic fracture characteristics of the ice body; the initial defects of the ice body have been avoided to the most significant degree possible when preparing the ice beam samples. Hence, this assumption has little effect on the accuracy of the simulation results.

Figure 3 displays the sample size, numerical model, and experimental setup. Two fixed rigid bodies support the ice beam, while the other rigid body is loaded downward at a speed of 0.08 mm/min at the center of the upper surface of the beam. Spagni et al. [49] conducted three friction studies between stainless steel pin planes and ice of varying roughness. The results showed that the friction coefficients of the ice and the stainless steel were 0.02–0.11. The friction coefficient between the ice beam and the stiff body, where both are in universal contact, is set to 0.045 in this research. C3D8R is used to mesh the ice beam in the modeling process, and zero-thickness cohesive elements are introduced in the same manner described in Section 2.2.1.

Table 1. Lists an ice beam’s cohesive element and C3D8R material properties [48,50,51].

Element Type	Parameter	Symbols	Value
	Young’s modulus	$E$	9100 MPa
C3D8R	Poisson’s ratio	$\nu$	0.3
	Density	$\rho$	900 kg/m3
	Fracture energy	$G_n^c,G_s^c,G_t^c$	0.0008 N/mm
COH3D8	Peak tractive force	$t_{no},t_{so},t_{to}$	1 MPa
	Stiffness	$k_{no},k_{so},k_{to}$	320 MPa/mm

displays the material properties of the zero-thickness cohesive element and the C3D8R.

This work models ice beams using three different element sizes: 4 mm, 4.5 mm, and 5 mm. It is advised that the number of elements in the cohesive zone is between 2 and 10, depending on the sample size [52,53]. NI Baoyu et al. [54] claim that to characterize the material’s fracture process more precisely, there must be at least four cohesive elements inside a cohesive zone length. In order to satisfy the accuracy requirements, the cohesive zone in the ice beam modeling has to contain at least six cohesive elements. Figure 4 compares the experimental and numerical results of the ice beam’s three-point bending load-deflection relationship, demonstrating that the numerical load-deflection curve for the ice beam with an element size of 4 mm is mainly consistent with the experimental results. The numerical peak load is 141 N, with a compression deflection of 0.17 mm. The relative inaccuracy is 5.53%. As the element size decreases, the slope of the load-deflection curve and the peak load decrease. Similar to the results of Hauke Herrnring et al. [55], in Figure 5, the stress cloud is displayed. As seen in Figure 5a,b, compressive stresses were created in the specimen’s upper layer and tensile stresses in its lower layer during the three-point bending loading, according to Figure 5c, which corresponds with the experiment’s center fracture and further supports the model’s validity, when loading increases, the cohesive element at the ice beam’s middle bottom first experiences damage before being destroyed.

3.2. Numerical Simulation of SCB

According to Bishoy and Yavas et al. [56,57,58,59], the SCB sample comprises a narrow beam at the bottom, an ice layer, and a wide beam at the top. Figure 6 displays the sample dimensions and numerical model. The bottom, somewhat longer, narrow beam was subjected to a central loading auger velocity of 0.02 mm/s, while the top wide beam was fixed in an appropriate location on the test apparatus. The cohesive element is embedded between the bottom narrow beam and the ice layer. In contrast, the top wide beam and the ice layer are considered to adhere perfectly without failure, which is modeled this way because, according to the experimental results, the debonding of the ice from the aluminum plate occurs between the bottom narrow beam and the ice layer around the negative fracture phase angle. The ice material characteristics are the same as in Table 1, and the cohesive element parameters are displayed in

Table 2. SCB and Push-out numerical simulation of material parameters [56,57,58,59].

Material Name	Parameter	Symbols	Value
	Young’s modulus	$E$	73,100 MPa
Al	Poisson’s ratio	$\nu$	0.33
	Density	$\rho$	2750 kg/m3
	Fracture energy	$G_n^c,G_s^c,G_t^c$	0.0005 N/mm
Cohesive behavior (SCB)	Peak tractive force	$t_{no},t_{so},t_{to}$	0.8 MPa
	Stiffness	$k_{no},k_{so},k_{to}$	100,000 MPa/mm
	Fracture energy	$G_n^c,G_s^c,G_t^c$	0.004 N/mm
Cohesive behavior (Push-out)	Peak tractive force	$t_{no},t_{so},t_{to}$	1.7 MPa
	Stiffness	$k_{no},k_{so},k_{to}$	250 MPa/mm

. Ice and aluminum are modeled using linear elastic materials and C3D8R. The SCB debonding process begins when the peak stress between the interfaces exceeds a threshold value and continues until debonding is completed.

In this experiment, the SCB samples are modeled using four element sizes: 2.5 mm, 3 mm, 3.5 mm, and 4 mm. Figure 7 compares the SCB experiments’ load-displacement relationship to the numerical results. The results reveal that the numerical load-displacement curves are primarily compatible with the actual results, and element size has little effect on simulation results. When the element size is 4 mm, the simulated peak load is 18 N, with a displacement of 0.0288 mm and a relative error of 8.9%. Figure 8 shows the SCB ice–aluminum debonding process. Figure 8a shows that the cohesive element begins to suffer damage, and the debonding between the ice and the aluminum plate begins. Figure 8b shows that the cohesive element has been destroyed, and the ice and aluminum plate have been wholly debonded. Comparing the experimental phenomena [56,57,58,59], it is clear that the simulation fits with the experimental results, confirming the model’s reliability. Validated model parameters will be used for wing electric-impulse de-icing simulations.

3.3. Numerical Simulation of Push-Out

According to Bishoy [56] and Yavas et al. [59], the push-out experiment uses an aluminum rod with an ice cylinder adhered around it. Figure 9 displays the sample size and numerical model. During the experiment, a velocity of 0.02 mm/s was provided to the top of the aluminum rod, gradually pushing it out of the ice cylinder until it was entirely out. The lower surface of the ice cylinder was fixed. A cohesive element is embedded between the aluminum bar and the ice, and a C3D8R is used to simulate the aluminum bar and ice. Table 2 displays the material properties of the cohesive element, ice, and aluminum bar.

The push-out experiment is modeled using four element sizes: 0.6 mm, 0.8 mm, 1 mm, and 1.2 mm. Figure 10 shows a comparison of numerical and experimental load-displacement relationships. The findings indicate that element size has essentially little effect on calculation results. The simulated and experimental load-displacement curves are often in good agreement. The simulated peak load was 165 N. The load-displacement at peak load was 0.0086 mm, with a relative error of 1.9%. The stress cloud in Figure 11 shows the process of releasing ice from an aluminum bar. Figure 11a,b depicts the cohesive elements from the start of destruction to ultimate disappearance. It simulates the start and completion of debonding in ice cylinders and aluminum rods.

3.4. Numerical Simulation of Electro-Impulse De-Icing on Wing Leading Edge Structure

This study uses a model NACA 0012 wing with a chord length of 1000 mm. In Figure 12, the wing profile line is displayed. Only the de-icing of the wing’s leading edge is studied because it is primarily covered with ice. The electric current distribution on the skin’s surface near the coil is dense because of the skin effect of the electric current, and the skin force reduces as the electric current falls along the skin’s thickness. Only the skin thickness is added at the coil’s relevant point because adding the entire wing thickness would lead to further issues; this additional skin layer is also referred to as the doubler. In a subsequent study, the impulse force acted on the surface of the doubler. In both the span and chord directions, the ice size is close to what was seen in the experiment. The ice thickness in modeling is fixed at a constant value. The error was acceptable in the Sommerwerk et al. [51] study comparing steady-thickness ice modeling with two other modeling techniques. By offsetting the wing profile lines, the ice and leading edge of the wing were modeled in three dimensions. The ice thickness is 4 mm, the skin thickness is 0.5 mm, and the span length of the leading-edge portion of the wing is 500 mm. Figure 13 illustrates the doubler’s position at the wing’s leading edge. A total of 60% of the chord length of the leading edge is where doubler one’s center is situated, and 80% of the chord length is where doubler two’s center is situated. A total of 24,892 C3D8R with an element size of 4 mm are used in the doubler’s finite element modeling and the wing’s leading edge. Globally embedded cohesive elements simulate the ice layer that adheres to the wing’s leading edge; there are 19,402 cohesive elements, and 9800 solid elements, which have a dimension of 4 mm. Cohesive elements are again embedded between the ice layer and the wing contact surface to simulate the shedding of the ice layer on the wing surface after the ice layer and wing are combined. A total of 9800 cohesive elements are embedded between the ice layer and the wing contact surface. The uniform pressure on the surface of the doubler simulates the impulse force generated by the coil. Compared to the electromagnetic coupling analysis, this approach reduced analyzing time and demonstrated efficacy in the Möhle et al. [60] investigation. The wing skin was attached to the ribs and spar caps with appropriate boundary conditions determined by the experimental structure. Figure 14 displays the established finite element model of the wing icing.

3.5. Results and Discussions

Sommerwerk et al. [61] conducted wing electric-impulse de-icing experiments in an icing wind tunnel under actual icing conditions, obtaining the deflection curves of the doubler one center skin and de-icing results with impulse action time and comparing the simulation and experimental results for the subsequent verification of the numerical model of wing electric-impulse de-icing. In order to adjust the uniform pressure operating on the doubler region, numerical simulations were first performed in a without-ice environment. For the without-ice situation, an initial voltage is chosen. Figure 15a shows experimental and numerical simulation comparisons of the skin deflection curves at the center of doubler one when doubler one is pulsed once. The displacement cloud is shown in Figure 16. The maximum deflection distance of the experimental clean skin is 2.05 mm, while the maximum deflection distance of the simulated clean skin is 2.04 mm, both with a relative error of 0.5%. Next, the de-icing process is simulated using the established finite element model of wing icing. Under the wing ice-covering situation, one, two, and three pulses are applied to the doubler one, and the experimental and numerical simulations of the skin deflection curves at the center of the doubler one are compared in Figure 15b–d. The ice layer increases the structure’s stiffness, reducing the skin deflection distance under impulse action. The comparison of experimental and numerical simulation results shows that the loads acting on the doubler are appropriate.

Finally, based on the verified load amplitudes operating on the doubler’s surface, two, four, and seven impulses were applied to the doubler one’s surface, with a comparison of experimental and simulated de-icing outcomes given in Figure 17. The fourth impulse load amplitude is equal to the third impulse load amplitude in the ice-covered condition, and the fifth, sixth, and seventh impulse amplitudes are equal to the first impulse amplitude in the without ice state. The de-icing effect of the electric impulse is quantitatively described by the de-icing rate, which is defined as the ratio of the area of the ice-shedding region to the original total area. In the numerical simulation, this is equivalent to the ratio of the number of cohesive elements deleted between the ice-wing skin to the total number. A higher de-icing rate indicates better de-icing performance. Under two impulses, skin deformation is small. However, de-icing may be discovered at the ice edge downstream of the airfoil, with an experimental de-icing rate of 20% and a simulation de-icing rate of 25% compared to the initially iced-over area. The experimental skin was deformed with four impulses, and a substantial amount of ice was removed, leaving only 10% of the ice on the wing, resulting in a 64% simulated de-icing rate. The experiment with seven impulses resulted in additional de-icing, eventually leaving 5% of the ice on the wing and a simulated de-icing rate of 67%. Compared to Sommerwerk et al.’s simulation work [51], the number of de-icing impulses in this study is consistent with the experiment, and the skin deflection curve and de-icing effect are more accurate.

4. Analysis of Influential Factors

Variations in the ice–skin interface shear adhesion strengths, doubler loading positions, and impulse sequences significantly impact the de-icing effect in the wing electro-impulse de-icing process simulation. Extensive research on these parameters is required to shed light on how various parameters affect the de-icing impact and serve as a guide for designing the anti-de-icing system that comes after.

4.1. Impact of Varying Shear Adhesion Strengths on the Rate of De-Icing

The shear adhesion strength at the ice–skin contact is estimated in this paper to be between 1.1 and 2.3 MP. The de-icing rate of the electric impulse at the ice–skin interface is shown in Figure 18 for various shear adhesion strengths. According to the simulation results, there is a significant connection between the shear adhesion strength fluctuation at the ice–skin interface and the de-icing impact. Specifically, increased shear adhesion strength causes the de-icing rate to decline from 66% to 56%. As a result, the parameter index of the shear adhesion strength of the ice–skin contact requires special consideration during the design and simulation calculations of the anti-deicing system. On the one hand, precise measurement of the shear adhesion strength at the ice-skin interface under various icing meteorological conditions can aid in increasing the effectiveness of electric-impulse de-icing, lowering the system’s energy consumption, and enhancing the precision of numerical simulations used to prevent de-icing. On the other hand, the effectiveness of electric-impulse de-icing can be significantly enhanced by applying surface treatment technologies, superhydrophobic material preparation, and other methods to reduce the ice interface shear adhesion strength artificially. More thorough testing and research in related fields is essential because the broad application of electric-impulse de-icing technology and other technologies to lower ice adhesion strength to create a new kind of composite de-icing system has the potential to improve de-icing efficiency further and reduce energy consumption.

4.2. Effect of Varying Doubler Loading Positions on the Rate of De-Icing

Since there is more than one impulse coil in the actual electric-impulse de-icing system installation on the aircraft, it is essential to consider how other positional coils may affect the de-icing effect. Doing so will aid in designing and optimizing the multi-coil electric-impulse de-icing system. Figure 19 shows the de-icing rate when doubler two receives the same impulse load as doubler one. When doubler two receives two, four, and seven impulses, the de-icing rates are 59%, 71%, and 71%, respectively. Significant increase in de-icing rate compared to doubler one. When an airplane is in flight, the thickness and area of the icing on its lower surface are always more significant than those on its top surface. Additional coils should be placed on the wing’s lower surface or used more for de-icing to save energy and improve efficiency.

4.3. Effect of Varying Impulse Sequences on the Rate of De-Icing

An electric-impulse de-icing system changes the sequence of impulses acting on doublers one and two. The optimal de-icing effect can be achieved by modifying the impulse sequence given to doublers one and two. The impulse interval parameter α represents the impulse sequence acting on doublers one and two. The doubler one starts to be loaded α ms after the electric impulse is applied to the doubler two. The larger α is, the longer the impulse interval between doublers one and two. In this work, A ranges from 0 to 80 ms, and the de-icing rates at various impulse sequences are presented in Figure 20. Depending on the impulse sequence α, the residual response of doubler one to doubler two and the superposition of doubler two’s impulse response fluctuate, resulting in varying de-icing rates. When α < 20 ms, the superposition effect is better due to the shorter impulse interval, and as α increases, the superposition effect worsens, and the deicing rate decreases; when α = 20 ms, the superposition effect is the best, and the deicing rate reaches a peak of 100%; and when α > 20 ms, the superposition effect worsens due to the longer impulse interval, so the deicing rate decreases as α increases. An impulse interval of 20 ms is recommended when deicing an airplane. At this point, the impulse response of the two coils can provide a better superposition effect and a better deicing effect.

5. Conclusions

This paper first demonstrates how the cohesive element can simulate ice fracture and characterize ice debonding with the matrix material using ice beam three-point bending experiments, SCB, and push-out experiments. The effects of different element sizes on the computational results were compared in the numerical calculations. The cohesive element dramatically affects the computational results when simulating ice beam fracture. In contrast, the cohesive element has a negligible effect on the computational results when characterizing the debonding of the ice from the substrate material and has almost no effect in the push-out experiments. Then, the finite element model of wing electric-impulse de-icing was established, and its reasonableness was verified by comparing the skin deflection curves and de-icing rate after the doubler was loaded through experiments and simulations, as well as the feasibility of using uniform pressure to simulate the impulse force.

Finally, the effects of various ice–skin interface shear adhesion strengths, doubler loading positions, and impulse sequences on de-icing were studied. It has been discovered that (1) varied ice–skin interface shear adhesion strengths significantly impact the de-icing effect, and the de-icing rate steadily declines from 66% to 56% as shear adhesion strength increases. To achieve the goal of reducing energy consumption and improving de-icing rates, it is recommended that the anti-de-icing system and the actual de-icing process be designed to accurately obtain the real-time ice–skin interfacial shear adhesion strength and reduce the interfacial shear adhesion strength of the ice layer using surface treatment technology or superhydrophobic material preparation. (2) When two, four, and seven impulses were delivered to doubler two, the de-icing rates were 59%, 71%, and 71%, respectively, showing a significant increase in de-icing efficiency over when the impulses were delivered to doubler one. It is recommended to place more coils on the wing’s lower surface or use more lower surface coils for de-icing during the de-icing process. (3) Depending on the impulse sequence α, doublers one and two impulse responses are overlaid differently, resulting in varying de-icing rates when α = 20 ms, the superposition effect is the best, and the de-icing rate reaches 100%. If multiple coils must function simultaneously during the de-icing process, the impulse sequence should be shortened suitably to increase the de-icing rate and efficiency.

This paper uses the cohesive zone model in material fracture and delamination damage, establishing a finite element model of wing electric-impulse de-icing that considers ice–skin interface debonding and ice fracture and the factors that influence the de-icing rate. In the future, more emphasis should be placed on developing electromagnetic-force-thermal multi-field coupling models in simulation research and developing hybrid anti-icing technology in de-icing system design, in which active and passive anti-icing technologies work together.