# Dynamic model of an electromechanical de-icing system coupling fracture and resonance

Jason Pothin, Aurélien Reysset, Marc Budinger, Valérie Pommier-Budinger, G Boscoletto

► To cite this version:

Jason Pothin, Aurélien Reysset, Marc Budinger, Valérie Pommier-Budinger, G Boscoletto. Dynamic model of an electromechanical de-icing system coupling fracture and resonance. Smart Materials and Structures, 2025, 34 (11), pp.115043. 10.1088/1361-665X/ae1a24. hal-05488919

HAL Id: hal-05488919

https://ut3-toulouseinp.hal.science/hal-05488919v1

Submitted on 2 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

PAPER • OPEN ACCESS

# Dynamic model of an electromechanical de-icing system coupling fracture and resonance

To cite this article: J Pothin et al 2025 *Smart Mater. Struct.* **34** 115043

View the [article online](#) for updates and enhancements.

## You may also like

- ICRH modelling of DTT in full power and reduced-field plasma scenarios using full wave codes

A Cardinali, C Castaldo, F Napoli et al.

- Thermal conductivity of suspended MBE-grown PtSe₂

Juliette Jolivet, Arkadiusz P Gertych, Eva Desgué et al.

- Global evidence that cold rocky landforms support icy springs in warming mountains

Stefano Brighenti, Constance I Millar, Scott Hotaling et al.
# Dynamic model of an electromechanical de-icing system coupling fracture and resonance

J Pothin<sup>1,2,*</sup> ID, A Reysset<sup>1</sup>, M Budinger<sup>1</sup>, V Pommier-Budinger<sup>2 ID</sup> and G Boscoletto<sup>2</sup>

E-mail: jason.pothin@isae-supaero.fr

Received 18 June 2025, revised 7 October 2025

Accepted for publication 31 October 2025

Published 24 November 2025

<figure><img src="image_3.png" /></figure>

CrossMark

## Abstract

Preventing ice accumulation on aircraft surfaces is a real challenge in industrial aviation. Notably, ice buildup significantly impacts the functionality of wings and control surfaces, thereby diminishing overall aircraft performance. Existing in-flight ice protection solutions, such as the use of hot air from engines or electrical resistance for surface heating, as well as the use of inflatable pneumatic boots, are characterized by either high energy consumption or maintenance requirements. Moreover, it is anticipated that the aerospace sector will undergo a substantial transition toward electric technology in the coming decade. In response to this shift, ongoing research and development efforts explore electromechanical systems to mitigate ice accumulation. This study concerns resonant electromechanical systems that reduce the power required for de-icing by exploiting the phenomenon of resonance, which naturally amplifies the vibrations used to break up and expel ice, but whose frequency varies with ice shedding, which makes it difficult to control. The formulation of a representative dynamic model of the actuation system and the ice fracture mechanism would enable control laws to be developed to operate the de-icing system at its optimum operating point, even if the latter is variable. Consequently, this study aims to develop a numerical model capable of accurately simulating the dynamics of resonant electromechanical de-icing systems. At the end of the study, this model will be used to tune the excitation signal controlling the de-icing system in open-loop.

Keywords: dynamic model, de-icing, electromechanical system, fracture, model-in-the-loop (MIL), resonance, HIL (hardware-in-the-loop)

## 1. Introduction

Under certain atmospheric conditions, ice can accumulate on structures and compromise their integrity.

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Furthermore, ice accumulation negatively affects the aero-dynamic performance of various components such as rotor blades and airfoils, leading to higher energy consumption and associated costs. If not prevented, the resulting changes in an aircraft's lift characteristics, particularly stall speed, would pose considerable safety risks.

Literature highlights two categories of systems devised to mitigate these issues: anti-icing and de-icing technologies. Anti-icing systems are engineered to activate before encountering icing conditions, thereby preventing ice formation continuously. Conversely, de-icing systems allow a certain

<footnote n="1"><sup>1</sup> Institut Clément Ader (ICA), Université de Toulouse, CNRS-INSA-ISAE-Mines Albi-UPS, Toulouse, France</footnote>

<footnote n="2"><sup>2</sup> Fédération ENAC ISAE-SUPAERO ONERA, Université de Toulouse, Toulouse, France</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
thickness (a few millimeters) of ice accumulation before its removal.

Recent research has indicated that resonant electromechanical systems could become a viable solution for de-icing systems [1–5]. Specifically, the method investigated in this paper involves piezoelectric ceramic actuators that induce vibrations in the structure near its natural structural frequencies. By exploiting resonance, this approach achieves the displacement amplitude necessary to remove ice while minimizing power consumption. However, the resonant frequency can change rapidly and significantly with the breaking and shedding of the ice, although it is of utmost importance to operate the system at its resonant frequency. This change in resonant frequency is due to the fracture, as several properties, such as the modal stiffness and modal mass, change in a non-linear way. Studies on the variation of modal parameters in non-linear systems include those by Wang and Wang [6], who examined nonlinear vibrations in matrix-cracked composite beams where effective stiffness and mass depend on crack evolution, by Yang *et al* [7], who proposed efficient methods for solving vibration equations with time-varying parameters, and by González *et al* [8] who presented a simplified approach for analyzing vibrations in plates with variable stiffness. While these studies do not directly address ice fracture, they illustrate how systems can vary with internal physical changes, like how fractures in an electromechanical de-icing system can locally alter mass and stiffness.

The objective of this paper is to propose a dynamic model that could be used to test different controller architectures to drive a resonant piezoelectric de-icing system at its resonant frequency. The study does not address the problem of controller synthesis in closed-loop, which is another scientific issue, but aims to provide a complete dynamic model of the de-icing system. The model includes not only the piezoelectric actuation system bonded to the structure but also the ice-fracture mechanism within the ice or at the ice-structure interface. The dynamics of the structure and the actuation system are represented by lumped parameters, while the ice fracture dynamics is modeled by its fracture propagation velocity. This complete lightweight transient model is capable of representing the main physics and the dynamic behavior of the resonant de-icing system in a framework particularly well-suited for controller synthesis. The novelty of the study lies in the combination of a dynamic model of the actuation system with a dynamic model of the ice fracture mechanism, which has never been proposed before. A further contribution of the article is the utilization of the model for the open-loop control of the ice protection system, employing a straightforward strategy. The system is controlled by means of a frequency sweep signal, the range of which extends from the frequency of the iced structure to the frequency of the de-iced structure. For this purpose, the response of the ice protection system to a frequency sweep will be studied.

Following a review of the state of the art in modeling piezoelectric structures and ice fracture propagation, the article begins with a description of the system studied and used to experimentally validate the dynamic model developed in the article. The hypotheses considered when formulating the

model, in particular those related to fracture mechanisms, are also described. The dynamic models of the structure and fracture mechanisms are then developed. The implementation in Simulink is detailed, and the proposed model is validated by comparing the propagation estimated by the model with that obtained from experiments. Finally, the model is used to calculate a law for determining the duration of a frequency sweep signal, which can then be used to control the ice protection systems with variable resonance in open-loop.

## 2. State of the art in modeling piezoelectric structures and fracture mechanisms

Piezoelectric structures have been extensively studied and modeled for actuation [9] and vibration control [10] purposes. While some studies present couplings with multiple physics, such as hydraulics [11], coupling with fracture phenomena is much less common. Several analytical or numerical models of the electromechanical impedance of structures incorporating piezoelectric transducers (PZTs) have been proposed by [12, 13] for detecting structural damage and by [14] for detecting ice on stay cables. However, these models are not relevant for describing the fracture dynamics as required in the de-icing application control.

In the field of fracture mechanics, several numerical methods can be used, such as the extended finite element (FE) method [15], the virtual crack closure technique (VCCT) [16] and cohesive zone models (CZMs) [17]. Gastaldo performed an analysis of these methods in her thesis [18]. The VCCT method has been used to predict mechanical de-icing of structures [19, 20] or to determine ice adhesion properties [21]. CZM has emerged as a viable method for modeling the propagation of ice adhesion failure. Chen *et al* Milanese [22] used CZM to predict adhesion failure of an ice/aluminum interface, considering interfacial roughness and initial damage, with results demonstrating CZM's efficiency when supported by precise experimental data. Sommerwerk and Horst [23] investigated crack propagation in ice-adherent structures using a FE CZM framework in quasi-static and dynamic conditions. All these works are of the out-of-loop types, i.e. with no possible interaction with a control law. In [24], Mohanty *et al* proposed a phase-field model for crack growth in electromechanically coupled functionally graded piezoelectric ceramics. It allows simulating crack initiation and propagation without predefined crack paths, offering a powerful tool for studying fracture in smart materials. However, this approach remains computationally expensive and is designed for quasi-static problems, making it unsuitable for testing control strategies in fast dynamic systems such as resonant piezoelectric de-icing applications. In [25], the authors presented a study on the implementation of a hybrid test control for a crack. The model used is an inverse model designed to develop a control strategy for crack propagation on a concrete sample with a multi-axial testing machine. For security and set-up design reasons, the experiment is first replaced by a FE simulation of crack propagation. The approach proposed in this work allows coupling with control loops but is restricted to
quasi-static problems due to the computational cost of FE sim-
ulations. There is therefore a need for efficient methodologies
to represent fast dynamics, such as resonant structures com-
bined with fracture, as observed in the de-icing application.

In this study, the modeling methodology is based on
the coupled stress-energy-based approach, introduced by
Leguillon [26] and already used to analyze ice fracture
mechanisms [25], and a dynamic model of the fracture
propagation, coupled with a lumped model of the resonant
structure and actuation system. The main features of this mod-
eling methodology are:

* Incorporation of a dynamic fracture evolution, enabling a more accurate prediction of de-icing performance under operational conditions. The paper will demonstrate the application of the coupled criterion (CC) in a state flow to represent multiple fracture paths with initiation or propagation phases.

* A link between fracture propagation, structural dynamics and system actuation, providing a framework for control synthesis. It should be noted that fractures introduce strong nonlinearities that can induce complex temporal responses.

* A reduction in the computational time, facilitating experi-
mental validation (more simulations can be carried out in
less time, reducing the number of validation tests). The
proposed approach is based on off-line FE calculations,
enabling real-time Hardware-in-the-Loop (HIL) tests.

**3. System and ice fracture mechanism under study**

3.1. *System under study*

Although the presented methodology may be extended to
complex structures, a simple use case will be studied as an
example in this study: a titanium plate equipped with two
PIC 181 piezoelectric ceramics (figure 1). The titanium plate
(130 × 50 × 1 mm³) is in free boundary conditions on its four
edges. Its surface is covered with a 3 mm-thick ice layer. The
ice is modeled as an isotropic material with a Young's mod-
ulus of 9 GPa (according to [27]), a density of 900 kg m⁻³
and a Poisson's ratio of 0.3 (according to [28]), which cor-
responds to the properties of the freezer ice used in the tests.
The material of titanium chosen is the Ti-6Al-4 V, which has
a Young's modulus of 110 GPa, a density of 4510 kg m⁻³, and
a Poisson's ratio of 0.33 (according to [29]).

3.2. *Ice fracture and resonant mode*

When examining the de-icing mechanism using a resonant electromechanical de-icing system, two distinct types of fractures can occur: cohesive fractures (or bulk fractures) within the ice block and adhesive fractures at the interface between the ice and the substrate. Cohesive fractures are not always present during de-icing and, if they occur without adhesive fractures, they prevent ice shedding as the ice can remain adhered to the surface. Adhesive fractures are crucial for ice de-bonding.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Geometry of the plate.</figcaption></figure>

Two main types of modes can be used for resonant deicing: flexion and extension. A performance comparison between these two types of modes is provided in [2]. Extensional modes permit complete surface de-icing on a single mode, with a typical resonance frequency variation of up to 30%, a stable energy release rate, and generally operate at higher frequencies (10 kHz and above). In contrast, flexural modes provide a more limited de-icing capability when relying on a single mode, unless specific constructive precautions are taken [30]. They generally exhibit greater resonance frequency variations (approximately 50%) and a decreasing energy release rate with fracture length, making them more difficult to control. They operate in a lower frequency range (1–10 kHz). This work focuses on flexural modes, as they require lower deicing power than extensional modes and may be compatible with a wider range of actuator technologies. Moreover, as operating frequencies are lower, implementation in apid Control Prototyping (RCP) and HIL environments will be easier. This makes flexural modes a suitable and relevant preliminary case study to demonstrate the relevance of the proposed approach. In this study, the first flexural mode is considered the sole source of fracture initiation and propagation. For the structure under study, the modal shape is given in figure 2.

3.3. *Criteria for studying ice fracture mechanism*

An evaluation of the parameters that define the mechanical
properties of ice, along with the fracture criteria that govern
ice failure, is essential for this study.

The critical values used in this study are given in table 1
where σ<sub>c</sub> is the critical value of the tensile stress σ, τ<sub>c</sub> the
critical value of the shear stress τ dependent on the char-
acteristics of the ice/substrate interface, and G<sub>c</sub> the critical
value of the energy release rate G, also called toughness. The
tensile strength of ice ranges between 0.7 MPa and 3.1 MPa
[31–35], and is strongly influenced by temperature, strain rate,
sample volume, and ice grain size. The study adopted a value
of 3 MPa due to its good correspondence with the conducted
experiments. The critical energy release rate at the ice inter-
face also depends on interface characteristics, such as material
and roughness, as well as the ice properties themselves [36].
Andrews and Lockington [31] estimated an average value of
G<sub>c</sub> = 1.0 J m<sup>-2</sup>, which is commonly used in the literature and
adopted here for the cohesive fracture. Additionally, a value of
0.5 MPa for the critical shear strength is chosen for the compu-
tations, as it is widely reported for freezer ice interacting with
standard metal surfaces [37].
<figure><img src="image_6.png" /><figcaption>**Figure 2.** First flexural mode.</figcaption></figure>

<table><caption>**Table 1.** Critical values for freezer ice.</caption><thead><tr><th></th><th>σ<sub>c</sub> [MPa]</th><th>G<sub>c,coh</sub> [J m<sup>-2</sup>]</th><th>τ<sub>c</sub> [MPa]</th><th>G<sub>c,adh</sub> [J m<sup>-2</sup>]</th></tr></thead><tbody><tr><td>Literature</td><td>[1 – 6]</td><td>[1 – 3]</td><td>[0.5 – 4]</td><td>[0.1 – 2]</td></tr><tr><td>Used here</td><td>3</td><td>1</td><td>0.5</td><td>0.5</td></tr></tbody></table>

The CC of Leguillon [26, 38] states that a fracture occurs when both a stress criterion and an energy criterion are satisfied simultaneously. To this end, this article considers two stages in the fracture analysis.

The first stage is the initiation of fracture on a surface S where the CC is expressed as:

* $\sigma$ or $\tau$ greater than their critical value at any point on the expected crack surface: $\sigma(M) \ge \sigma_c$ or $\tau(M) \ge \tau_c \forall M \in S$

* a value of G computed on the expected crack surface greater than the critical value $G_c$: $G = -\frac{\Delta E_p}{S} \ge G_c$ with $E_p$ the elastic energy stored in the structure.

These two criteria generally vary in opposite ways as the crack surface expands. This leads to the formation of a crack of non-zero initial size when the vibration level becomes high enough.

The second stage is the crack propagation. In this case, the stress concentration at the crack tip enables the stress criterion to be satisfied. Consequently, only the energy release rate criterion needs to be observed. The propagation of the crack can continue if the value of G is greater than a critical value: $G = -\frac{\delta E_p}{\delta S} \ge G_c$.

For the remainder of this article, fracturing surfaces S will be rectangular in shape, of depth b and length a: $S = b \cdot a$.

### 3.4. Fracture mechanisms

Golovin *et al* [39] have shown that the interface length between the ice and the substrate influences the fracture mode and that the transition between these modes occurs at a critical length of this ice block. This was also observed by Gastaldo *et al* [1] who described the fracture mechanisms as a function of ice block size. In this study, the fracture process is modeled for two ice block sizes.

#### Long ice blocks

For this configuration, as illustrated in figure 3, the piece of ice attached to a structure is longer than a critical length, and the fracture mechanism is divided into the four stages described in figure 3:

* (A) Initiation of the cohesive fracture: The cohesive fracture begins at the top of the ice block at the center of the plate. Its

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Fracture mechanism for long ice blocks.</figcaption></figure>

initiation is governed by the tensile stress and energy criteria that must be satisfied simultaneously CC. The initial crack length is denoted $a_{coh,init}$.

* (B) Propagation of the cohesive fracture: After cohesive fracture initiation, the fracture propagates through the thickness towards the ice/structure interface (figure 3) and is governed by the energy criterion.

* (C) Initiation of the adhesive fracture: The adhesive fracture is assumed to be symmetric with respect to the plane of the cohesive fracture (figure 3). The adhesive fracture of initial length $a_{adh,init}$ (figure 3) is initiated from the cohesive fracture and is governed by the CC (shear stress and energy criteria).

* (D) Propagation of the adhesive fracture: The energy criterion governs adhesive fracture propagation.

#### Short ice blocks

In this configuration, as shown in figure 4, the length of ice attached to a structure is considered to be less than a critical length, so that only adhesive fracture occurs (see [1]). This represents the two last states C) and D) presented for the long ice block configuration, with adhesive fracture starting from the edges.

### 3.5. Dynamic model of the structure with piezoelectric actuators

The model of a structure with k piezoelectric actuators/sensors that is reduced to n degrees of freedom can be written as a set of equations (1), which come from [2].

$$ \begin{cases} [M]_{nn} \{\ddot{x}\}_n + [D_s]_{nn} \{\dot{x}\}_n + [K]_{nn} \{x\}_n = [N]_{nk} \{V\}_k, \\ \{Q_k\} = [N]_{kn}^T \{x\}_n + [C_0]_{kk} \{V\}_k \end{cases} \quad (1) $$

$[M]_{nn}$ represents the mass matrix of the system, $[D_s]_{nn}$ the damping matrix of the system, $[K]_{nn}$ the stiffness matrix, $[N]_{nk}$ the electromechanical coupling matrix of the actuators/sensors with structure, $\{V\}_k$ the actuator voltage vector, $\{Q\}_k$ the actuator electric charge vector, $[C_0]_{kk}$ the actuator/sensor capacitance matrix, and $\{x\}_n$ the system displacement vector. The model parameters can be easily calculated analytically or using FE software, except for the damping matrix. For the sake
<figure><img src="image_8.png" /><figcaption>Figure 4. Fracture mechanism for short ice blocks.</figcaption></figure>

of simplicity, this term is calculated using the formula:

$$[D_s]_{nn} = \frac{1}{\{Q_m\}_n} \sqrt{[M]_{nn} [K]_{nn}} \quad (2)$$

where $Q_m$ is the quality factor of the mode that can be measured experimentally.

The structure is analyzed using modal analysis for the different stages of the ice fractures. Since the ice blocks change during the fracture—they may split into several pieces and some may fall off—the terms $[M]_{nn}$, $[K]_{nn}$ and $[N]_{nk}$ vary depending on the ice fracture and must be computed for each fracture state described in the next section, both for the initiation and propagation phases.

In this study, the structure is analyzed for the first flexural mode. It is represented with one degree of freedom (thus $n = 1$) and one out-of-plane modal displacement $x_{\text{mod}}$. In particular, the point at the center of the plate, which is the point of maximum deformation for the first flexural mode, can be selected.

### 3.6. Models for the fracture mechanism

#### Long ice blocks

##### 3.6.1. Cohesive fracture initiation.
The ice layer is discretized into $N$ ice blocks in its thickness (figure 5), resulting in $N$ fracture states, each labeled by the index $i$. $N$ was determined through a convergence study of strain energy and stresses, including shear stress. For the cohesive interface, $N = 31$ points were used, and the same number of points was applied for the adhesive layer to ensure consistent resolution and accuracy. This calculation scheme corresponds to the full FE approach presented in [40]. The results are coupled here to the dynamics of a resonant structure.

For stress evaluation, the contacts between the ice blocks are kept. The critical displacement $x_{\text{coh,init},\sigma,i}$ (for the initiation of cohesive fracture) is obtained when the critical tensile stress $\sigma_c$ is reached in the block $i$. Since the displacements are proportional to the stresses, they can be computed from the modal values with:

$$x_{\text{coh,init},\sigma,i} = \frac{\sigma_c}{\sigma_{\text{mod},i}} x_{\text{mod},i} \quad (3)$$

where $x_{\text{mod},i}$ is the modal displacement at state $i$ and $\sigma_{\text{mod},i}$ the tensile stress in section $i$.

For energy evaluation, the contacts between the ice blocks are removed (without removing the ice blocks). Then the potential energy linked to elastic deformation $E_{p,\text{mod},i}$ is calculated at each state $i$. Only the potential energy stored in the ice

<figure><img src="image_9.png" /><figcaption>Figure 5. Discretization for cohesive and adhesive fracture.</figcaption></figure>

is considered. The critical displacement $x_{\text{coh,init},G,i}$ is obtained when the critical energy release rate $G_c$ is reached in the block $i$. Since the displacements are proportional to the square root of the energy, it can be computed from the modal values:

$$x_{\text{coh,init},G,i} = \sqrt{\frac{G_c}{G_{\text{mod,init},i}}} x_{\text{mod},i} \quad (4)$$

with

$$G_{\text{mod,init},i} = - \frac{(E_p)_{\text{mod},i} - (E_p)_{\text{mod},0}}{b \cdot a_i} \quad (5)$$

where $b$ is the depth of the system, $a_i$ represents the crack length for each fracture state and $(E_p)_{\text{mod},0}$ represents the potential energy without fracture.

The two curves corresponding to these displacements are plotted in figure 6. The critical displacements were evaluated using a FE numerical simulation with post-processing, following the approach described by Doitrand et al [40], referred to as the full FE implementation of the CC. The fracture initiation occurs when the CC is satisfied, i.e. when the critical tensile stress and the critical energy release rate are exceeded, i.e. when the two curves intersect. The abscissa of the intersection point gives the initial crack length $a_{\text{coh,init}}$ and the ordinate the critical displacement $x_{\text{coh,init}}$.

##### 3.6.2. Cohesive fracture propagation.
Only the energy criterion is considered for this phase. The discretization principle in $n$ blocks is the same as for initiation. This time, however, the ice block corresponding to the fracture state is removed at each step.

The critical displacement $x_{\text{coh,prop},G,i}$ obtained when the critical energy release rate $G_c$ is reached in the block $i$ is given by:

$$x_{\text{coh,prop},G,i} = \sqrt{\frac{G_c}{(G_{\text{mod,prop},i})}} x_{\text{mod},i} \quad (6)$$

where

$$G_{\text{mod,prop},i} = - \frac{(E_p)_{\text{mod},i} - (E_p)_{\text{mod},i-1}}{b \cdot (a_i - a_{i-1})} \quad (7)$$

with $(E_p)_{\text{mod},i-1}$ the potential energy at the past fracture state and $a_{i-1}$ the associated past fracture state.

Figure 7 plots the critical displacement $x_{\text{coh,prop},G,i}$ versus the cohesive fracture length during the fracture propagation. It shows that the fracture can fully propagate as the required displacement decreases with the fracture length.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Cohesive fracture initiation—Critical displacements versus cohesive fracture length—Study case with a 3 mm-thick ice layer.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Cohesive fracture propagation—Critical displacements versus cohesive fracture length—Study case with a 3 mm-thick ice layer.</figcaption></figure>

### 3.6.3. Adhesive fracture initiation and propagation.
The modeling approach is similar for the initiation and propagation of adhesive fractures. The ice layer is discretized into *n* blocks along the ice/structure interface. The initiation of the adhesive fracture is modeled by removing the contacts between the ice block and the interface, and the propagation is modeled by removing the ice block while the fracture propagates (figure 8).

To analyze the mechanism involved in the initiation of the adhesive fracture, the critical displacements $x_{\text{adh,init},\sigma,i}$, $x_{\text{adh,init},\tau,i}$ and $x_{\text{adh,init},G,i}$ obtained by exceeding, respectively, the critical tensile stress $\sigma_c$, the critical shear stress $\tau_c$ and the critical energy release rate $G_c$ in the block *i* are calculated similarly as for the cohesive fracture.

Figure 9 shows the three curves corresponding to these critical displacements. The initiation of the adhesive fracture occurs when a CC is satisfied. This is the CC with normal strength and energy that is first met. The abscissa of the intersection point gives the initial crack length $a_{\text{adh,init}}$ and the ordinate the critical displacement $x_{\text{adh,init}}$ (both for adhesive initiation).

To analyze the mechanism involved in the propagation of the adhesive fracture, the critical displacement $x_{\text{adh,prop},G,i}$ (for adhesive propagation) obtained by exceeding the critical

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Example of two steps of the adhesive propagation model.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Adhesive fracture initiation—Critical displacements versus adhesive fracture length—Study case with a 65 mm-long ice layer (Half of the beam is represented).</figcaption></figure>

energy release rate $G_c$ in the block *i*, is computed and plotted in figure 10. This curve shows that, in contrast to cohesive fracture propagation, the critical displacement required for fracture propagation increases with the fracture length (i.e. stable fracture). It also indicates that complete delamination of the ice is not achievable. Rafik et al [30] show how to optimize substrate thickness to promote delamination and mitigate this phenomenon.

### Short ice blocks

In this case, only two fracture states are considered: the initiation and propagation of an adhesive fracture. The method for determining the critical displacements for these two states is similar to that for the long ice block configuration. The initiation of adhesive fractures is governed by the CC with shear stress and energy. The critical displacements calculated for the propagation of adhesive fractures are shown in figure 11. The CC leads to an instantaneous initial fracture over 13% of the interface length for a vibration magnitude of 15 µm. The energy criterion for the propagation of the adhesive fracture shows complete delamination for the same amplitude. Thus, the adhesive fracture starts under the effect of the CC and is completed thanks to the energy stored in the ice. The curve associated with the energy-stored criterion is presented in figure 12.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Adhesive fracture propagation—Critical displacements versus adhesive fracture length—Study case with a 65 mm-long ice layer (Half of the beam is represented).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Adhesive fracture initiation: critical displacements versus adhesive fracture length—Study case with a 30 mm-long ice layer (Half of the beam is represented).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Adhesive fracture propagation: critical displacements versus adhesive fracture length—Study case with a 30 mm-long ice layer (Half of the beam is represented).</figcaption></figure>

# 4. Simulink dynamic model

## 4.1. Model implementation

Figure 13 shows an overall diagram of the dynamic model, with a block representing the structure with the piezoelectric actuators and a block representing the fracture mechanism. The Simulink model consists of three parts (figure 14):

<figure><img src="image_17.png" /><figcaption>**Figure 13.** General representation of the dynamic model of the electromechanical de-icing system.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Simulink model of the de-icing system.</figcaption></figure>

* The fracture mechanism model, which includes two or four stages depending on the ice block size compared to the critical length. This section presents only the case of a large ice block (i.e. four stages). The outputs of each stage are the required displacement to initiate or propagate the fracture and the length of the fracture *a* when propagation is studied. The flowchart presented in figure 14 works as follows. Once the displacement of the plate to initiate the cohesive fracture is reached, its propagation begins and continues in the ice layer until the structure/ice interface is reached (i.e. ice thickness). The end of the cohesive fracture triggers the initiation of the adhesive fracture. Once initiated, the adhesive fracture progresses until the fracture length limit is reached. This limit is the value at which infinite energy is required to propagate the fracture.

* Tables of data generated from modal analyses. These tables contain the parameters $M(a)$, $K(a)$, and $N(a)$ calculated for each fracture stage and each fracture length $a$.

* The dynamic model of the structure with piezoelectric actuators. This model is based on (1), providing outputs such as the central displacement of the plate $x$. In this model, the parameters $M(a)$, $K(a)$, and $N(a)$ come from the tables above.

The blocks for initiating and propagating fractures are implemented similarly for cohesive and adhesive fractures.
The initiation block is mainly based on a comparator block between x and $x_{coh,init,σ,i}$ or $x_{adh,init,G,i}$. Once $x_{c,cohesive}$ or $x_{c,adhesive}$ is reached, the computation continues with the propagation block.

The propagation block is detailed in figure 15. This block is initiated with the index i such that $a_i = a_{coh,init}$ when studying cohesive fractures and with the index i such that $a_i = a_{adh,init}$. It contains an integrator that starts integrating at an average speed of 1000 m s⁻¹ from the initial value $a_{coh,init}$ or $a_{adh,init}$. This model is based on the assumption that the velocity of fracture propagation can be estimated from the speed of wave propagation (Rayleigh). Indeed, Schreurs and Eindhoven [41] estimates the fracture propagation velocity in ice to be about 20% to 40% of the speed of sound $v_s = \sqrt{\frac{G_{icc}}{\rho}}$, yielding $v = 1000 \text{ m s}^{-1}$. It is assumed that the crack tip velocity remains constant. Doitrand [42] has indeed demonstrated that the temporal shape of this velocity has minimal influence. Furthermore, numerical simulations of the proposed model were carried out for propagation speeds ranging from 100 to 3000 m s⁻¹. Little variation in numerical results was observed for values above 300 m s⁻¹. Consequently, a propagation speed of 1000 m s⁻¹ is used for the rest of the study.

## 4.2. Simulation of the fracture initiation and propagation

In the simulation model, a frequency sweep around the resonance mode is used to drive the de-icing system. The voltage amplitude V is computed from the electromechanical coupling relation 8:

$$V = \frac{M\omega_0^2}{NQ_m}x_c, \quad (8)$$

where $x_c$ is the modal displacement, M the modal mass, N the force factor, $\omega_0$ the resonance pulsation, and $Q_m$ the mechanical quality factor. In the simulations, V is set to 100 V for a plate completely covered by ice and 200 V for the case of a plate partially covered by ice (ice block of length 30 mm).

The frequency range results from FE analyses of the iced and non-iced structures. For the plate completely covered by ice, the frequency with ice is around 1100 Hz, and around 550 Hz without ice, so a range of 550 Hz seems adequate to capture all the frequency variations. However, a safety margin is included in case of errors in variation, a greater-than-expected ice thickness (which increases stiffness and thus the frequency), or imperfect ice (which is often the case in reality). The frequency range is therefore set to $\Delta F = 700 \text{ Hz}$. The frequency sweep is carried out by decreasing frequencies to follow the resonant frequency shift of a flexural resonance mode in the presence of fractures. Finally, the frequency range is set to [1200–500] Hz. For the plate partially covered with ice, using a similar approach, the frequency range sweep is set to [800–500] Hz.

Figure 16 shows the simulated evolution of the adhesive and cohesive fractures for a plate entirely covered by ice. The simulation is carried out for frequencies decreasing from 1200 Hz to 500 Hz. The simulation describes the expected mechanism:

<figure><img src="image_19.png" /><figcaption>Figure 15. Fracture propagation model (Stages 2 and 4 of figure 14).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 16. Long ice block configuration: simulated adhesive and cohesive fracture length versus time and versus frequency.</figcaption></figure>

* **First:** At 0.151 s (1092 Hz), the initiation and propagation of the cohesive fracture over the entire thickness of the ice layer occur, alongside the initiation of the fracture.

* **Second:** At 0.67 s (730 Hz), the adhesive fracture propagates over 20.4 mm in both directions, corresponding to 37% of the plate. The fracture model was implemented only for half of the plate. Therefore, the propagation length of 20.4 mm corresponds to a single direction. By symmetry, the total crack propagation is 40.8 mm, which is consistent with the experimental observations. Due to the presence of the fracture and the remaining ice, the resonance frequency drops to around 600 Hz.

Figure 17 shows the same kind of simulated results for a plate partially covered with ice. This time, only the evolution of the adhesive fracture is represented. The simulation is performed for frequencies decreasing from 800 Hz to 500 Hz. The initiation of adhesive fracture starts at 0.55 s and for 719 Hz and propagates over 15 mm and the frequency drops to 557 Hz (at 1.7 s).

## 4.3. Validation of the fracture model

The experimental procedure is detailed below. The PZT ceramics used are PIC 181 from PI Ceramics and are bonded to the structure using Loctite Ablestik 57 °C silver epoxy to ensure reliable mechanical and electrical contact. The
<figure><img src="image_21.png" /><figcaption>**Figure 17.** Short ice block configuration: simulated adhesive fracture length versus time and versus frequency.</figcaption></figure>

electrodes are configured to excite the flexural mode, with polarization and field direction optimized for maximal electromechanical coupling. The electrical impedance and vibratory response were measured to determine the resonance frequency, the mechanical quality factor ($Q_m$), and the electromechanical coupling factor ($N$). The vibrations were measured using a Polytech vibrometer. Ice deposition was performed carefully to minimize thermal stress, keeping the water near 0 °C and adjusting the substrate temperature gradually. Fracture initiation and propagation were recorded using a standard 60 fps camera. PZT excitation was applied around the resonance frequency to induce bending vibrations and study ice detachment.

### Long ice blocks

First, the configuration with the plate completely covered by ice is studied. A 100 V sinusoidal sweep signal ranging from 1200 Hz to 500 Hz is applied for 1 s to observe the ice fracture. Figure 18 shows the result of the de-icing test. The cohesive fracture occurs throughout the entire thickness of the ice. The adhesive failure is partial.

The frequency at which the fracture occurs is well represented in the model (figure 16). The cohesive fracture appears at a resonance frequency of around 1092 Hz for the simulation and around 1105 Hz for the test, and drops to 600 Hz at the end of the 1 s simulation and around 620 Hz at the end of the 1 s test, thus with the same dynamics. Each colored part in the first picture of figure 18 can be compared to the fracture length simulated chart in figure 16.

The final length of the adhesive failure is measured on the two half-plates adjacent to the cohesive failure. It is equal to 25 mm on the first half of the plate and 20 mm on the second half. These values should be compared with the simulated value of 20.4 mm (see figure 18) and show that the model is faithful to the tests.

### Short ice blocks

A 30 mm-long ice block partially covers the plate. In this case, the resonance frequency was measured at 712 Hz (close to the simulated value of 719 Hz). A 200 V sinusoidal signal ran-

<figure><img src="image_22.png" /><figcaption>**Figure 18.** Pictures of the sample and measurement of the adhesive crack length after a de-icing test at 100 V (long ice case).</figcaption></figure>

ging from 800 Hz to 500 Hz is applied for 2 s to observe the de-icing. Figure 19 shows the result of the de-icing test. As expected, a total adhesive fracture is observed, starting from the edges of the ice block and propagating to the center of the block. Regarding the dynamics of the fracture mechanism, the final state of adhesive fracture is reached at 1.56 s in the test and at 1.6 s in the simulation. Each colored part in the first picture of figure 19 can be compared to the fracture length simulated chart in figure 17. Again, the model accurately represents the fracture mechanism.

## 5. Optimization of the frequency sweep signal

### 5.1. Discussion on the sweep time

To get the best performance of the resonant de-icing system, the duration of the sweep signal must be optimized. If the frequency sweep is too fast, there is a risk that the amplitude required for de-icing will not be obtained, as it takes time for a resonant system to reach its steady state. This risk is confirmed by the simulation of figure 20 where the ratio of power to voltage squared (allowing the power to be normalized, whatever the voltage [43]) is represented. Power is the active power absorbed by the transducer during frequency sweeps, obtained by the product of voltage and current, filtered by a low-pass filter of order 2, with a cut-off frequency of 200 Hz. Figure 20 shows that a frequency sweep of 700 Hz generated in 0.5 s does not allow any fracture. This is also confirmed by experimental testing.

Thus, a study of the minimum sweep time required to reach a certain level of acceptable vibration amplitude ($x_r$) compared to the maximum possible vibration amplitude ($x_{max}$) for a given resonance mode is essential. In this section, the objective is to establish a relationship between the normalized response amplitude ($\frac{x_r}{x_{max}}$) of a resonant system and the duration of the frequency sweep excitation signal. Even if second-order systems are well-known, the response to a frequency sweep is not straightforward since the Laplace transform of a frequency sweep cannot be expressed analytically in a simple
<figure><img src="image_23.png" /><figcaption>**Figure 19.** Pictures of the sample and measurement of the adhesive crack length after a de-icing test at 200 V (small ice case).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 20.** Simulation of a de-icing test performed in 0.5 s: no fracture appears.</figcaption></figure>

manner. It is thus not obvious to establish an analytical expression between the response of a resonant system and the sweep time.

## 5.2. Computation of the sweep time

Considering the dynamic equation of the system (1), shown in section 3.5, in one dimension, it can be demonstrated that it corresponds to a harmonic oscillator, which can be expressed as the well-known equation as follows:

$$ \ddot{x} + \frac{\omega_0}{Q_m} \dot{x} + \omega_0^2 x = \frac{N}{M} V \quad (9) $$

where $\omega_0$ is the angular frequency of the mode, $Q_m$ is its quality factor, $N$ is its electromechanical coupling, and $M$ is its modal mass, all detailed in [2].

For forced oscillations with a sinusoidal input signal of constant angular frequency $V = V_0\sin(\omega_0 t)$, the solution to the harmonic equation (9) [44] demonstrates that the vibratory amplitudes of the system increase until they reach a maximum amplitude $x_{max}$ (10):

$$ x_{max} = \frac{NV_0Q_m}{M\omega_0^2} \quad (10) $$

The envelope of the vibratory amplitudes has the form of a first-order step response characterized by the time constant $\tau$:

$$ \tau = \frac{2Q_m}{\omega_0} \quad (11) $$

<table><caption>**Table 2.** Variables of interest for the control of the resonant de-icing system.</caption><thead><tr><th>Variable</th><th>Name</th><th>Unit</th><th>Dimension</th></tr></thead><tbody><tr><td>x<sub>r</sub>/x<sub>max</sub></td><td>Response amplitude</td><td>/</td><td>/</td></tr><tr><td>ΔT/ΔF</td><td>Sweep ramp</td><td>s<sup>2</sup></td><td>T<sup>2</sup></td></tr><tr><td>τ = 2Q<sub>m</sub>/ω<sub>0</sub></td><td>Time constant</td><td>s</td><td>T</td></tr></tbody></table>

In the case of a frequency sweep input signal, the parameter $\tau$ is also crucial for defining the temporal performance of the system controlled by a sweep frequency signal. Other parameters are also required to establish a general relationship between vibration level and sweep time, such as the frequency range $\Delta F$, the sweep time $\Delta T$, and the system response required to de-ice $x_r$, which is estimated thanks to a FE Analysis. These parameters are listed in table 2 with their units and dimensions.

Buckingham's theorem [45] and [46] is used to establish a relation between all these parameters. According to this theorem, if a system is described by a physical equation of the form:

$$ f(x_1, x_2, \dots, x_n) = 0 \quad (12) $$

where $x_i$ are the variables of the problem, then the equation can be rewritten using dimensionless variables denoted $\pi_i$ in the form:

$$ F(\pi_1, \pi_2, \dots, \pi_p) = 0. \quad (13) $$

For an initial equation with $n$ variables and $k$ physical dimensions involved (time, length, mass, ...), then the number of dimensionless variables needed is $p = n - k$.

Using the Buckingham theorem, it is thus possible to compute two non-dimensional numbers $\pi_1$ and $\pi_2$ linking the variables of interest for the control of resonant de-icing systems:

$$ \pi_1 = \frac{x_r}{x_{max}} \quad (14) $$

and

$$ \pi_2 = \frac{1}{T_s^2} \frac{\Delta T}{\Delta F} \quad (15) $$

Now that the number of variables has been reduced to two, $\pi_1$ can be plotted as a function of $\pi_2$ to obtain a relationship between the variables of interest. For this purpose, a design of experiments is computed and presented in table 3. Then simulations are run, $\pi_1$ is plotted as a function of $\pi_2$ (figure 21) and finally, a relationship between $\pi_1$ and $\pi_2$ is computed by regression (16).

$$ \pi_1 = \frac{x_r}{x_{max}} = \frac{1}{1 + \frac{1}{A\pi_2^2 + B\pi_2 + C}} \quad (16) $$

where $A = 0.0319$, $B = 1.2033$ and $C = 0.4955$.

It is also possible to reverse the calculation: determine the sweep time required to reach a given $x_r/x_{max}$ knowing the frequency range and the time constant of the resonant de-icing
<table><caption>**Table 3.** Simulation parameters used in the DOE.</caption><thead><tr><th>Symbol</th><th>Name</th><th>Values</th><th>Unit</th></tr></thead><tbody><tr><td><i>f</i><sub>0</sub></td><td>Resonant freq.</td><td>500–4000</td><td>Hz</td></tr><tr><td><i>Q</i><sub>m</sub></td><td>Quality factor</td><td>10–200</td><td>—</td></tr><tr><td>Δ<i>T</i></td><td>Sweep time</td><td>1, 2</td><td>s</td></tr><tr><td>Δ<i>F</i></td><td>Sweep width</td><td>300, 700</td><td>Hz</td></tr></tbody></table>

<figure><img src="image_25.png" /><figcaption>**Figure 21.** Data and estimated curve for $\pi_1$ function of $\pi_2$.</figcaption></figure>

system. For example, in the case of the long ice, with a $\Delta F = 700$ Hz and a time constant $\tau = 0.021$ s, a sweep time $\Delta T$ of 1 s is needed to reach $\frac{x_r}{x_{max}} = 85\%$ as shown in figure 22. With a sweep time of 0.5 s, the obtained displacement is only 74% of the steady state, which explains why no fracture appears in the simulation of figure 20.

## 6. Real-time application: RCP and HIL

To validate the accuracy and suitability of the proposed dynamic model for the controller development phases, tests are carried out using two approaches: RCP and HIL.

RCP and HIL utilize real-time execution on dedicated hardware platforms, like Speedgoat Baseline and DSpace. In RCP, this hardware interfaces directly with the physical electromechanical system, allowing rapid prototyping and experimental validation. In contrast, HIL replaces the physical plant with a high-fidelity, real-time simulation model of the system, enabling the controller to interact with a virtual plant under real-time constraints. This allows for thorough testing, including fault conditions, to be carried out without risking damage to the hardware. HIL integrates real hardware components with simulated models, to provide robust and repeatable validation of control systems for complex real-time embedded applications [47]). For the resonant de-icing system under study, HIL is particularly useful for validating closed-loop resonance tracking using vector control or adaptive algorithms. This validation could not be performed using transient FE models (e.g. CZM) due to the significant computational cost of these models.

The tests presented here represent frequency sweeps at low and high voltage levels. They enable the model to be evaluated under fractured and non-fractured conditions. A two-phase

<figure><img src="image_26.png" /><figcaption>**Figure 22.** Application of the equation (16) on the long ice case with $\Delta F = 700$ Hz.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 23.** RCP: Power consumption over voltage squared filtered during the de-icing.</figcaption></figure>

control sequence is implemented: detection is first performed at low voltage (20 V) over a wide frequency range (700 Hz) to locate the resonance, followed by de-icing at high voltage (100 V) over a narrower frequency range (200 Hz) around the previously identified resonant frequency. Each phase lasts one second, and the full cycle is repeated twice, resulting in a total duration of four seconds. This control sequence was implemented in real-time RCP and HIL testing. The RCP and HIL test setup includes the following hardware:

* Real-time machine (RTM): Speedgoat Unit real-time target machine with the IO191 I/O module, a 16-bit analog input and output module.

* Workstation (WS): Dell Inspiron 16, 13th Generation Intel® Core i7-13 800 H, 2500 MHz, 14 core.

To illustrate the capacity of the model to reproduce ice fractures in real time, the graphs in figures 23 and 24 show the ratio of active power to voltage squared.

Figure 23 with RCP shows the actual behavior observed in the experimental setup. The peaks visible in the graph correspond to the fractures that occurred during the frequency sweeps. Remarkable is the similarity between the HIL results and the RCP measurements, which demonstrates the accuracy of the dynamic model proposed in this study.
<figure><img src="image_28.png" /><figcaption>**Figure 24.** HIL: Power consumption over voltage squared filtered during the de-icing.</figcaption></figure>

Concerning the real-time performance of the model, for the HIL and RCP tests, the sampling frequency of the RTM was set to 20 kHz, corresponding to a sampling period of 50 µs. Among the various tasks executed by the RTM, the maximum execution time recorded during the tests was 33.84 µs, resulting in a CPU usage of approximately 68%. This leaves a reasonable margin for potential operation at higher frequencies if needed.

## 7. Conclusion

This study presents the development and validation of a dynamic model for a resonant electromechanical de-icing system. The model, derived from modal FE analyses and based on stress and energy release rate calculations, integrates mechanical resonance, fracture dynamics, and electromechanical coupling in a simple but representative model. The model enables the temporal simulation of fracture propagation and the evolution of the resonance frequency due to fracture and ice detachment that occur when the system's first resonance mode is activated. Experiments were carried out to validate this dynamic model, which provides a reliable basis for the development of control laws. The model was used to optimize a frequency sweep signal that controls the resonant electromechanical de-icing system. A law was established to minimize sweep time while maintaining system performance.

Several types of ice exist, such as rime or glaze ice, which can be characterized by their Young's modulus and density [48]. Pure glaze ice was used to simplify the study, as producing rime ice is more challenging. However, this approach can be extended to other types of ice, both in the numerical model (by adjusting the properties of the ice) and in experiments (by using an icing wind tunnel to produce different types of ice). Also, this approach could be used to study de-icing electromechanical systems combined with a coating that reduces the ice adhesion to the substrate [49].

In future work on modeling issues, the dynamic model could be refined to consider the influence of other resonant modes and the effect of losses and temperature on de-icing mechanisms. The model presented has been implemented for

de-icing using flexural mode conditions. It could also be of interest to adapt this work to extensional modes that have demonstrated their de-icing efficiency.

Future work should also address the synthesis of controllers for such resonant de-icing systems, as this is the main application of the dynamic model proposed in this study. Vector control could be a good candidate for tracking resonance frequencies [50]. This kind of controller could be validated on the system studied in this article and could then be used for NACA profiles found on wind turbines, aircraft, or UAVs. The dynamic model of the ice fracture mechanisms in such structures is more complex than that of plates. However, we believe that the control strategy for tracking resonance frequencies, which has been validated on plates, could apply to a wide range of structures.

## Data availability statement

All data that support the findings of this study are included within the article (and any supplementary files).

## Acknowledgments

The authors would like to acknowledge the MeGEP doctoral school, the Institut Supérieur de l'Aéronautique et de l'Espace, the Institut Clément Ader, and the Institut de Recherche et Technologies (IRT) with the project FREEZING 2.

## Author contributions

J Pothin ID 0009-0004-6398-3400

Conceptualization (equal), Data curation (equal), Formal analysis (equal), Investigation (equal), Methodology (equal), Software (lead), Validation (lead), Visualization (equal), Writing – original draft (lead)

A Reysset

Conceptualization (equal), Funding acquisition (equal), Methodology (equal), Resources (equal), Software (equal), Supervision (equal), Validation (equal), Writing – review & editing (equal)

M Budinger

Conceptualization (lead), Formal analysis (equal), Funding acquisition (equal), Investigation (supporting), Methodology (equal), Project administration (equal), Software (supporting), Supervision (lead), Validation (supporting), Visualization (equal), Writing – review & editing (equal)

V Pommier-Budinger ID 0000-0003-3949-3749

Conceptualization (equal), Formal analysis (equal), Funding acquisition (lead), Methodology (supporting), Project administration (lead), Software (supporting), Supervision (lead), Validation (supporting), Writing – review & editing (lead)

G Boscoletto

Investigation (equal), Methodology (equal), Software (equal)
## References

[1] Gastaldo G, Palanque V, Budinger M and Pommier-Budinger V 2022 Stress and energy release rate influence on ice shedding with resonant electro-mechanical de-icing systems *ICAS 2022-33rd Congress of the Int. Council of the Aeronautical Sciences*

[2] Pommier-Budinger V, Budinger M, Rouset P, Dezitter F, Huet F, Wetterwald M and Bonaccurso E 2018 *AIAA J.* **56** 4400–11

[3] Palanque V 2022 Design of low consumption electro-mechanical de-icing systems *PhD Thesis Université de Toulouse ISAE-Supaero*

[4] Palanque V, Pothin J, Pommier-Budinger V and Budinger M 2024 *Chin. J. Aeronaut.* **37** 92–103

[5] Jomaa M, Lévy P E, Vasic D, Costa F and Ali M 2024 *Smart Mater. Struct.* **33** 055043

[6] Wang Y and Wang Y 2025 *Arch. Appl. Mech.* **95** 118

[7] Yang J, Xie J, Chen G and Chen J 2021 *J. Low Freq. Noise, Vib. Act. Control* **40** 1744–63

[8] Andérez González J and Vescovini R 2023 *J. Compos. Sci.* **7** 30

[9] Ghenna S, Giraud F, Giraud-Audine C and Amberg M 2017 *IEEE Trans. Ind. Electron.* **65** 4880–8

[10] Moheimani S O R 2003 *IEEE Trans. Control Syst. Technol.* **11** 482–94

[11] Ling M, He X, Wu M and Cao L 2022 *IEEE/ASME Trans. Mechatronics* **27** 4942–50

[12] Sepehry N, Shamshirsaz M and Bastani A 2010 Electromechanical impedance modeling of cantilever beams with embedded PZT transducers 2010 IEEE/ASME Int. Conf. on Advanced Intelligent Mechatronics (IEEE) pp 569–74

[13] Park G and Inman D J 2007 *Phil. Trans. R. Soc. A* **365** 373–92

[14] Winston H, Sun F and Annigeri B 2001 *J. Eng. Gas Turbines Power* **123** 353–8

[15] Belytschko T and Black T 1999 *Int. J. Numer. Methods Eng.* **45** 601–20

[16] Rybicki E and Kanninen M 1977 *Eng. Fract. Mech.* **9** 931–8

[17] Elices M, Guinea G, Gomez J and Planas J 2002 *Eng. Fract. Mech.* **69** 137–63

[18] Gastaldo G 2025 Optimisation des systèmes de dégivrage hybrides basés sur des revêtements et des actionneurs électromécaniques *Doctoral Dissertation*

[19] Matbou Riahi M 2007 Numerical and experimental studies of the mechanical behaviour at the ice/aluminium interface *PhD Thesis Université du Québec à Chicoutimi*

[20] Xie D and Biggers S B 2006 *Finite Elem. Anal. Des.* **42** 977–84

[21] Palanque V, Villeneuve E, Budinger M, Pommier-Budinger V and Momen G 2022 *Cold Reg. Sci. Technol.* **204** 103679

[22] Chen Y, Wenhao W, Dong W, Li M and Lei L 2014 Numerical study on the adhesion strength between ice and aluminium based on a cohesive zone model *Adv. Appl. Math.* **1A**

[23] Sommerwerk H and Horst P 2017 *Eng. Fract. Mech.* **182** 400–24

[24] Mohanty S, Kumbhar P Y, Swaminathan N and Annabattula R 2020 *Smart Mater. Struct.* **29** 045005

[25] Jailin C, Carpiuc A, Kazymyrenko K, Poncelet M, Leclerc H, Hild F and Roux S 2017 *J. Mech. Phys. Solids* **102** 239–56

[26] Leguillon D 2002 *Eur. J. Mech. A* **21** 61–72

[27] Schulson E M 1999 *JOM* **51** 21–27

[28] Mellor M 1980 Mechanical properties of polycrystalline ice *Physics and Mechanics of Ice: Symposium Copenhagen (Technical University of Denmark)* (Springer) (Accessed 6–10 August 1979) pp 217–45

[29] Veiga C et al 2012 *Rev. Adv. Mater. Sci.* **32** 133–48

[30] Rafik Y, Palanque V, Budinger M, Pommier-Budinger V and Olivier P 2024 *Smart Mater. Struct.* **33** 075008

[31] Andrews E and Lockington N 1983 *J. Mater. Sci.* **18** 1455–65

[32] Yingkui G and Wenyuan M 2015 *J. North China Univ. Water Resour. Electr. Power* **36** 40–43

[33] Petrovic J 2003 *J. Mater. Sci.* **38** 1–6

[34] Druez J, Nguyen D and Lavoie Y 1986 *Cold Reg. Sci. Technol.* **13** 67–74

[35] Reich A, Scavuzzo R and Chu M 1994 Survey of mechanical properties of impact ice 32nd Aerospace Sciences Meeting and Exhibit p 712

[36] Huang X, Tepylo N, Pommier-Budinger V, Budinger M, Bonaccurso E, Villedieu P and Bennani L 2019 *Prog. Aerosp. Sci.* **105** 74–97

[37] Rønneberg S, Laforte C, Volat C, He J and Zhang Z 2019 *AIP Adv.* **9** 055304

[38] Martin E, Vandellos T, Leguillon D and Carrère N 2016 *Int. J. Fract.* **199** 157–68

[39] Golovin K, Dhyani A, Thouless M and Tuteja A 2019 *Science* **364** 371–5

[40] Doitrand A, Martin E and Leguillon D 2020 *Finite Elem. Anal. Des.* **168** 103344

[41] Schreurs P and Eindhoven G 2012 *Department of Mechanical Engineering* (Eindhoven University of Technology)

[42] Doitrand A, Molnar G, Leguillon D, Martin E and Carrère N 2022 *Eur. J. Mech. A* **93** 104483

[43] Budinger M et al 2021 Electromechanical resonant ice protection systems: energetic and power considerations *AIAA J.* **59** 2590–602

[44] Goldstein H, Poole C, Safko J and Addison S R 2002 Classical mechanics

[45] Buckingham E 1914 *Phys. Rev.* **4** 345

[46] Zohuri B 2017 *Dimensional Analysis Beyond the Pi Theorem* vol 1 (Springer) pp 85–128

[47] Karnopp D C, Margolis D L and Rosenberg R C 2012 *System Dynamics: Modeling, Simulation and Control of Mechatronic Systems* (Wiley) (https://doi.org/10.1002/ 9781118152812)

[48] Palanque V, Villeneuve E, Budinger M, Pommier-Budinger V and Momen G 2023 *Cold Reg. Sci. Technol.* **212** 103890

[49] Gastaldo G, Hernández-Rodríguez G, Pommier-Budinger V, Budinger M, Pothin J and Coclite A 2025 *Res. Eng.* **27** 106636

[50] Ghenna S, Giraud F, Giraud-Audine C, Amberg M and Lemaire-Semail B 2015 Vector control applied to a langevin transducer 2015 17th European Conf. on Power Electronics and Applications (EPE’15 ECCE-Europe) (IEEE) pp 1–9