Temperature Field and Thermal Stress Analysis of a Composite Wing Electric Heating System with Delamination Damage

1. Introduction

When an aircraft flies in icing weather conditions, ice forms on windward surfaces because of the freezing of supercooled water droplets. Ice accumulation on the wings of an aircraft can negatively impact its lift and manoeuvrability, increasing the risk of accidents and impacting flight performance. In severe cases, this can create safety hazards [1]; therefore, airworthiness standards require must be configured with a corresponding anti-icing system to the wings. Aerothermal anti-icing technology transfers hot engine air to the wing surface through the compressor's pilot line. Compared with traditional aerothermal anti-icing technology, electric anti-icing systems have the advantages of a fast response time, simple system design, flexible heating mode, and good compatibility with composite materials [2]. With the gradual increase in the use of composite materials in the main load-bearing structures of civil airliners [3], electrothermal anti-/de-icing technology has become one of the most promising anti-/de-icing technologies.

Owing to the anisotropic arrangement of the reinforcement material, the mechanical and thermomechanical properties of composite structures are characterized by significant anisotropy. Composite structures with electrically heated anti-icing systems are typically composed of multiple layers, including an insulation layer, a heating layer, and a bonding layer. This multilayer configuration further increases the complexity of the mechanical and thermal properties of the composite structure. Composite structures are susceptible to external impacts and other effects during service, which can cause delamination damage within laminated structures [4,5,6,7]. On the one hand, the multilayer structure formed after the addition of the electrothermal anti-icing system greatly increases the probability of delamination damage and is likely to lead to the emergence of interlayer delamination cracks, and the high temperature of electrothermal heating is prone to causing delamination because of the concentration of temperature and thermal stress [8]. On the other hand, delamination damage can lead to an uneven temperature distribution during electric heating, which affects the anti-icing efficiency. Therefore, the ability of composite electrothermal systems containing delamination damage in terms of the temperature field and mechanics needs to be studied urgently.

In most cases, delamination damage in composite materials is difficult to detect in time via conventional inspection methods, and composite components on aircrafts often continue to operate in a damaged state. This unrepaired damage can significantly affect the mechanical and thermal properties of the material, which in turn poses a potential threat to flight safety. Polla, A [9] et al. investigated the effect of heat on the mechanical properties of Carbon Fiber Reinforced Plastic(CFRPs) and reported that composites containing cracks are affected by temperature and that the absolute strength threshold of the matrix decreases as the temperature approaches the glass transition temperature, which affects the properties of the composite. Tay et al. [10] studied the mechanical properties of composites containing preset delamination and reported that residual thermal stresses have a greater effect on delamination near the surface of the composite. Furthermore, another study reported that composite laminates with preset damage exhibit approximately 14% lower bending strength than undamaged laminates, and experience higher surface temperatures under electric heating conditions[11]. There is an interaction between delamination damage and temperature distribution; thus, numerical calculations of the thermal multiphysics fields are necessary to study the electrothermal anti-icing structure of a wing subjected to delamination damage.

Numerical calculation is an important method for studying changes in material and structural properties in multiple physical fields and is widely used by scholars in calculations for wing anti-/de-icing systems[12,13,14]. Shen et al. [15] established a thermal conductivity equation for anisotropic composites, conducted numerical simulations on wings to investigate the impact of anisotropy on electrothermal anti-icing systems, and obtained the temperature of the anti-icing surface. Lou et al. [16] used a uniform heat flux boundary to heat the leading edge of the inner skin of a wing and calculated the thermal load on the wing skin on the basis of external flow fields and droplet impact; by performing coupled heat transfer calculations between the inner and outer skins, they obtained the temperature field on the de-icing surface. However, none of these methods can accurately simulate multiple delamination events and matrix crack damage in composite structures at the same time. The thermally extended layer-by-layer/3D solid element method used in this work not only simulates the composite material layer-by-layer but also accurately analyses the displacement and stress fields at the delamination leading edges and crack tips, and the two methods can be combined to improve the convergence speed by optimizing the algorithm and increasing the accuracy of the calculation.

In summary, delamination damage in composite materials prevents the structure from achieving its designed stiffness and strength during service. When combined with a multilayered anti-icing system, this will significantly increase the complexity and severity of the issue and put forward higher requirements for theoretical calculations and experimental implementation, which may be the main reason for the lack of relevant research at present. Therefore, this work integrates the study of delamination damage in composite materials with electrothermal anti-icing systems. It utilizes the thermally extended layer-by-layer/3D solid element method to explore the complex interactions between two features. More importantly, considering that delamination is more likely to occur in multilayered anti-icing and de-icing systems, this project conducts thermomechanical multiphysics numerical simulations to study the electrothermal anti-icing structure of composite wing leading edges with delamination damage. This investigation explores the impact of delamination damage at different locations on the temperature field, thermal stress, and displacement due to thermal deformation of the electrothermal anti-icing structure. These findings provide theoretical support for the application of electrothermal anti-icing systems in composite materials.

2. Composite Wing Leading Edge Modelling and Anti-Icing Environment Settings

2.1. Composite Wing Leading Edge Modelling

In this study, the model of the electric heating de-icing system referenced is based on experiments conducted by NASA [17]. This de-icing model employs a multilayered material structure, which consists of, from the top to the bottom, a wear-resistant layer (carbon fibre composite) with a thickness of 0.2 mm, followed by an adhesive layer (H77 adhesive) that is 0.28 mm thick, an electric heating layer (nickel‒chromium alloy) measuring 0.0127 mm, another adhesive layer measuring 0.28 mm, a support layer (glass fibre composite) with a thickness of 0.89 mm, and a thermal insulation layer (silicone foam) with a thickness of 3.43 mm. The physical parameters of each layer of material are shown in Table 1. Considering that the support layer material is anisotropic, it has different thermal conductivities, moduli of elasticity, coefficients of thermal expansion, and Poisson's ratios in the three directions, whereas the wear-resistant layer, adhesive layer, electric heating layer, and thermal insulation layer are isotropic and therefore have the same electrical conductivity in the three directions. Figure 1

Heating units are generally distributed in a regular pattern on wings. A multilayer thermal conduction structure model of a flat plate is set up with reference to the wing configuration, incorporating seven heating units labelled A, B, C, D, E, F, and G, as shown in Figure 2. The length of the flat plate is set to 1.8288 m and is designated the Y direction, whereas its width is 0.3937 m and is designated the X direction. The direction perpendicular to the plate is designated the Z direction. Heating units are placed in the central area of the model, with a length of 0.9144 m in the Y direction and a width of 0.19685 m in the X direction. The widths and heating powers for different zones are detailed in Table 2.

This paper uses the expansion layer-by-layer/3D solid element method [18,19,20,21,22,23,24,25] to investigate the electrothermal de-icing process of composite laminated panels. The upper part of the thermal insulation layer, which is prone to delamination damage, is investigated by the thermal expansion layer-by-layer method, and the thermal insulation layer, which is not prone to delamination damage, is investigated by the eight-node 3D solid element method for simulation. The equilibrium and coordination conditions between the thermal expansion layer-by-layer method and the eight-node 3D solid element method are based on compatibility conditions for temperature, displacement, and internal force balance. Finally, on the basis of the thermally extended layer-by-layer/3D solid element, the heat source and thermal convection boundary conditions are loaded into the model in the form of thermal loads, resulting in a thermodynamic model capable of simulating the electrothermal anti-icing process when applied to damaged composite laminates. The specific implementation process is shown in Figure 4.

The mechanical properties of the materials in each layer are shown in Table 3 and Table 4. The mechanical properties of the OCV-LTX1240 glass fibre unidirectional tape with known anisotropic parameters were chosen to represent the support layer, i.e., the glass fibre resin layer.

Figure 5 illustrates the coordinate system, geometric dimensions, and boundary conditions. To simulate the electric heating structure fixed to the wing by rivets,the displacement boundary condition surrounding the structure is established as fully fixed constraints during the analysis process.

2.2. Composite Material Anti-Icing Environment and Delamination Damage Setting

The meteorological conditions for aircraft icing used for the calculations in this paper are those of Case 67A [26], as shown in Table 5. The droplet collection coefficients and convective environment coefficients are taken from the experimental values.

To compare the effects of the location of delamination damage and the size of the delamination damage region on the anti-icing performance of the electrically heated structure, calculations were carried out in this work for four cases, including two delamination damage locations and two delamination damage areas, as shown in Figure 6.

2.3. Validation of Calculation Methods

Al-Khalil, Morency [17,26], and others have investigated the electrothermal anti-icing process for wings when the meteorological conditions are those of Case 67A using experimental and computational methods and have provided experimental and simulation results for water collection coefficients, overflow water distributions, and temperature distributions of the anti-icing surfaces. Therefore, this calculation condition is also chosen for use in this paper, as shown in Figure 7. The surface temperature distribution obtained by numerical simulation is very close to the change trend reported in the literature, and the simulation results and experiments are similar, which verifies the feasibility of the model and calculation method.

3. Finite Element Equations for Thermally Extended Layer-by-Layer/3D Solid Elements

The thermally extended layer-by-layer/3D solid element method was used to simulate the electrically heated anti-icing process for composite laminates [18]. A brief description of the method is given in this section.

3.1. Heat Mixed Variational Theorem

The heat mixing variational theorem is used as a starting point for the theoretical derivation [27]. The heat mixing variational principle is given by the following equation:

$$\begin{array}{c}{\displaystyle \underset{V}{\int }{\sigma }_{\alpha \beta }}\delta {\epsilon }_{\alpha \beta }\mathrm{d}V+{\displaystyle \underset{V}{\int }\left(\rho \frac{{\partial }^2{u}_{\alpha }}{\partial t^2}\right)}\mathrm{d}V-{\displaystyle \underset{V}{\int }\delta }{u}_{\alpha }{f}_{\alpha }\mathrm{d}V-{\displaystyle \underset{{\Gamma }_F}{\int }\delta }{u}_{\alpha }{F}_{\alpha }\mathrm{d}\Gamma \hfill \\ -{\displaystyle \underset{V}{\int }{h}_{\alpha }}\delta e_{\alpha }\mathrm{d}V-{\displaystyle \underset{V}{\int }\delta }\theta ({\Theta }^0\dot{\eta }-s)\mathrm{d}V-{\displaystyle \underset{\Gamma }{\int }\hslash }\delta \theta \mathrm{d}\Gamma =0\hfill \end{array}$$

(3-1)

where σαβ, εαβ, uα, and eα represent the stress tensor, strain tensor, displacement vector, and temperature gradient vector, respectively; $\eta$, ${\Theta }^0$, and $\rho$ represent the entropy density, constant positive reference temperature, and density, respectively; $\theta$ represents the temperature change relative to the reference temperature ${\Theta }^0$; fα and Fα represent the body force and surface force on the boundary ${\Gamma }_F$; s represents the internal heat source; hα represents the heat flux vector; and $\hslash$ represents the heat flux on the boundary.

The force and thermal boundary conditions can be written as follows:

$$\begin{array}{l}{\sigma }_{ij}{n}_j=\begin{array}{ccc}{F}_i& on& {\Gamma }_F\end{array}\\ u_i=\begin{array}{ccc}\overline{u_i}& on& {\Gamma }_H\end{array}\\ h_i{n}_i=\begin{array}{ccc}\hslash & on& {\Gamma }_H\end{array}\\ \theta =\begin{array}{ccc}\overline{\theta }& on& {\Gamma }_{\theta }\end{array}\end{array}$$

(3-2)

where $n_i$ is the unit normal vector of the boundary.

The gradient relationship between the strain field and the thermal field is given by the following equation:

$$\begin{array}{l}{\epsilon }_{ij}={\displaystyle \frac{1}{2}}(u_{i,j}+u_{j,i})\\ e_i=-{\theta }_i\end{array}$$

(3-3)

The constitutive equation is given by the following equation:

$$\begin{array}{c}{\sigma }_{ij}=C_{ijkl}{\epsilon }_{ij}-{\lambda }_{ij}\theta \hfill \\ \eta ={\lambda }_{ij}{\epsilon }_{ij}+\varrho \theta \hfill \\ h_i={\kappa }_{ij}{e}_j\hfill \end{array}$$

(3-4)

where$C_{ijkl}$ is the matrix of elastic stiffness constants,${\lambda }_{ij}$ is the stress temperature matrix, $\varrho$ is the coefficient of thermal expansion, and ${\kappa }_{ij}$ is the coefficient of thermal conductivity.

3.2. Finite Element Equation for Thermally Extended Layer-by-Layer/3D Solid Elements

The part above the thermal insulation layer is simulated using the thermal expansion layer-by-layer method. For a detailed description of the thermal expansion layer-by-layer method, please refer to [23]. The finite element equation utilized in this method is as follows:

$$\begin{array}{l}\left[\begin{array}{cc}{\mathit{M}}_{\alpha \beta \zeta \eta ke\kappa t}^{uu}& 0\\ 0& 0\end{array}\right]\left\{\begin{array}{c}{\ddot{\tilde{\mathit{U}}}}_{\beta \eta et}\\ 0\end{array}\right\}+\left[\begin{array}{cc}0& 0\\ {\mathit{C}}_{\beta \zeta \eta ke\kappa t}^{\theta u}& {\mathit{C}}_{\zeta \eta ke\kappa t}^{\theta \theta }\end{array}\right]\left\{\begin{array}{c}{\dot{\tilde{\mathit{U}}}}_{\beta \eta e\iota }\\ {\stackrel{.}{\tilde{\mathit{\Theta }}}}_{\eta e\iota }\end{array}\right\}+\\ \left[\begin{array}{cc}{\mathit{K}}_{\alpha \beta \zeta \eta ke\kappa \iota }^{uu}& {\mathit{K}}_{\alpha \zeta \eta ke\kappa \iota }^{u\theta }\\ 0& {\mathit{K}}_{\zeta \eta ke\kappa \iota }^{\theta \theta }\end{array}\right]\left\{\begin{array}{c}{\tilde{\mathit{U}}}_{\beta \eta en}\\ {\tilde{\mathit{\Theta }}}_{\eta en}\end{array}\right\}=\left\{\begin{array}{c}{\mathit{F}}_{\alpha \zeta k\kappa }^u\\ {\mathit{F}}_{\zeta k\kappa }^{\theta }\end{array}\right\}\end{array}$$

(3-5)

In the equation, $\kappa =m,s,h$; $\iota =n,g,f$; $m,n=1,2,\cdots N_E$;$s,g=1,\cdots N_E^P$; and$h,f=1,\cdots N_E^Q$. M, K, and C are the mass, stiffness, and damping matrices, respectively; see the Appendix A for specific expressions. $F^u$ and $F^{\theta }$ are mechanical and temperature load vectors, respectively, which can be expressed as follows:

$$\begin{array}{c}{\mathit{F}}^u={\displaystyle \underset{\Omega }{\int }\left({\mathit{q}}_b\mathit{\delta }{u}_3^{N+1}+{\mathit{q}}_{\mathit{t}}\mathit{\delta }{u}_3^0\right)}\mathrm{d}x\mathrm{d}y\hfill \\ {\mathit{F}}^{\theta }={\displaystyle \underset{\Omega }{\int }\left({\mathit{F}}_{\zeta k}{\mathit{\psi }}_m+h_t\mathit{\delta }{\mathit{\theta }}^0\right)}\mathrm{d}x\mathrm{d}y\hfill \end{array}$$

(3-6)

where ht is the heat flux from the upper surface of the composite material.

Eight-node 3D solid elements are used to model the thermal insulation layer with the following finite element formula for the 3D thermally coupled problem:

$$\left[\begin{array}{cc}{\mathit{M}}^{uu}& 0\\ 0& 0\end{array}\right]\left\{\begin{array}{c}\ddot{\mathit{u}}\\ 0\end{array}\right\}+\left[\begin{array}{cc}0& 0\\ {\mathit{C}}^{\theta u}& {\mathit{C}}^{\theta \theta }\end{array}\right]\left\{\begin{array}{c}\dot{\mathit{u}}\\ \dot{\mathit{\theta }}\end{array}\right\}+\left[\begin{array}{cc}{\mathit{K}}^{uu}& {\mathit{K}}^{u\theta }\\ 0& {\mathit{K}}^{\theta \theta }\end{array}\right]\left\{\begin{array}{c}\mathit{u}\\ \mathit{\theta }\end{array}\right\}=\left\{\begin{array}{c}{\mathit{F}}^u\\ {\mathit{F}}^{\theta }\end{array}\right\}$$

(3-7)

The total finite element control equations are assembled via the dynamic thermal method proposed by Li [28] et al. All the degrees of freedom (displacements and temperatures) in the control equations for the portion above the thermal insulation layer are divided into two groups: (1) the degrees of freedom that contact the thermal insulation layer (external degrees of freedom) and (2) the other degrees of freedom that do not contact the thermal insulation layer (internal degrees of freedom). Therefore, equations (3−5) can be transformed as follows:

$$\begin{array}{l}\left[\begin{array}{cccc}{\mathit{M}}_{11}^p& {\mathit{M}}_{12}^p& 0& 0\\ {\mathit{M}}_{21}^p& {\mathit{M}}_{22}^p& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{\ddot{\mathit{u}}}_1^p\\ {\ddot{\mathit{u}}}_2^p\\ {\ddot{\mathit{\theta }}}_1^p\\ {\ddot{\mathit{\theta }}}_2^p\end{array}\right\}+\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ {\mathit{C}}_{\theta u,11}^p& {\mathit{C}}_{\theta u,12}^p& {\mathit{C}}_{\theta \theta ,11}^p& {\mathit{C}}_{\theta \theta ,12}^p\\ {\mathit{C}}_{\theta u,21}^p& {\mathit{C}}_{\theta u,22}^{\mathit{p}}& {\mathit{C}}_{\theta \theta ,21}^p& {\mathit{C}}_{\theta \theta ,12}^p\end{array}\right]\left\{\begin{array}{c}{\dot{\mathit{u}}}_1^p\\ {\dot{\mathit{u}}}_2^p\\ {\dot{\mathit{\theta }}}_1^p\\ {\dot{\mathit{\theta }}}_2^p\end{array}\right\}\\ +\left[\begin{array}{cccc}{\mathit{K}}_{uu,11}^p& {\mathit{K}}_{uu,12}^p& {\mathit{K}}_{u\theta ,11}^p& {\mathit{K}}_{u\theta ,12}^p\\ {\mathit{K}}_{uu,21}^p& {\mathit{K}}_{uu,22}^p& {\mathit{K}}_{u\theta ,21}^p& {\mathit{K}}_{u\theta ,22}^p\\ {\mathit{K}}_{\theta u,11}^p& {\mathit{K}}_{\theta u,12}^p& {\mathit{K}}_{\theta \theta ,11}^p& {\mathit{K}}_{\theta \theta ,12}^p\\ {\mathit{K}}_{\theta u,21}^p& {\mathit{K}}_{\theta u,22}^p& {\mathit{K}}_{\theta \theta ,21}^p& {\mathit{K}}_{\theta \theta ,22}^p\end{array}\right]\left\{\begin{array}{c}{\mathit{u}}_1^p\\ {\mathit{u}}_2^p\\ {\mathit{\theta }}_1^p\\ {\mathit{\theta }}_2^p\end{array}\right\}=\left\{\begin{array}{c}0\\ {\mathit{F}}_2^p\\ 0\\ 0\end{array}\right\}+\left\{\begin{array}{c}{\mathit{T}}_1^p\\ 0\\ 0\\ 0\end{array}\right\}+\left\{\begin{array}{c}0\\ 0\\ {\mathit{Q}}_1^p\\ {\mathit{Q}}_2^p\end{array}\right\}\end{array}$$

(3-9)

where the superscript p refers to the extended layer-by-layer method. The subscript 1 indicates the degrees of freedom and loads (external degrees of freedom) that are at the interface between the extended layer-by-layer method and the 3D solid element method. Subscript 2 indicates degrees of freedom and loads that are not at the interface between the extended layer-by-layer method and the 3D solid element method (internal degrees of freedom). M, K, and C are the mass, stiffness, and damping of the thermomechanical extended layerwise method (TELM), all of which are given in specific form in the Appendix. The control equation (3−8) for the 3D solid element method can be similarly rewritten as follows:

$$\begin{array}{c}\left[\begin{array}{cccc}{\mathit{M}}_{11}^{\mathrm{s}}& {\mathit{M}}_{12}^{\mathrm{s}}& 0& 0\\ {\mathit{M}}_{21}^{\mathrm{s}}& {\mathit{M}}_{22}^{\mathrm{s}}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{\ddot{\mathit{u}}}_1^{\mathrm{s}}\\ {\ddot{\mathit{u}}}_2^{\mathrm{s}}\\ {\ddot{\mathit{\theta }}}_1^{\mathrm{s}}\\ {\ddot{\mathit{\theta }}}_2^{\mathrm{s}}\end{array}\right\}+\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ {\mathit{C}}_{\theta u,11}^{\mathrm{s}}& {\mathit{C}}_{\theta u,12}^{\mathrm{s}}& {\mathit{C}}_{\theta \theta ,11}^{\mathrm{s}}& {\mathit{C}}_{\theta \theta ,12}^{\mathrm{s}}\\ {\mathit{C}}_{\mathit{\theta }u,21}^{\mathrm{s}}& {\mathit{C}}_{\theta u,22}^{\mathrm{s}}& {\mathit{C}}_{\theta \theta ,21}^{\mathrm{s}}& {\mathit{C}}_{\theta \theta ,22}^{\mathrm{s}}\end{array}\right]\left\{\begin{array}{c}{\dot{\mathit{u}}}_1^{\mathrm{s}}\\ {\dot{\mathit{u}}}_2^{\mathrm{s}}\\ {\dot{\mathit{\theta }}}_1^{\mathrm{s}}\\ {\dot{\mathit{\theta }}}_2^{\mathrm{s}}\end{array}\right\}\hfill \\ +\left[\begin{array}{cccc}{\mathit{K}}_{uu,11}^{\mathrm{s}}& {\mathit{K}}_{uu,12}^{\mathrm{s}}& {\mathit{K}}_{u\theta ,11}^{\mathrm{s}}& {\mathit{K}}_{u\theta ,12}^{\mathrm{s}}\\ {\mathit{K}}_{uu,21}^{\mathrm{s}}& {\mathit{K}}_{uu,22}^{\mathrm{s}}& {\mathit{K}}_{u\theta ,21}^{\mathrm{s}}& {\mathit{K}}_{u\theta ,22}^{\mathrm{s}}\\ {\mathit{K}}_{\theta u,11}^{\mathrm{s}}& {\mathit{K}}_{\theta u,12}^{\mathrm{s}}& {\mathit{K}}_{\theta \theta ,11}^{\mathrm{s}}& {\mathit{K}}_{\theta \theta ,12}^{\mathrm{s}}\\ {\mathit{K}}_{\theta u,21}^{\mathrm{s}}& {\mathit{K}}_{\theta u,22}^{\mathrm{s}}& {\mathit{K}}_{\theta \theta ,21}^{\mathrm{s}}& {\mathit{K}}_{\mathit{\theta }\mathit{\theta },22}^{\mathrm{s}}\end{array}\right]\left\{\begin{array}{c}{\mathit{u}}_1^{\mathrm{s}}\\ {\mathit{u}}_2^{\mathrm{s}}\\ {\mathit{\theta }}_1^{\mathrm{s}}\\ {\mathit{\theta }}_2^{\mathrm{s}}\end{array}\right\}=\left\{\begin{array}{c}0\\ {\mathit{F}}_2^{\mathrm{s}}\\ 0\\ 0\end{array}\right\}+\left\{\begin{array}{c}{\mathit{T}}_1^{\mathrm{s}}\\ 0\\ 0\\ 0\end{array}\right\}+\left\{\begin{array}{c}0\\ 0\\ {\mathit{Q}}_1^{\mathrm{s}}\\ {\mathit{Q}}_2^{\mathrm{s}}\end{array}\right\}\hfill \end{array}$$

(3-10)

where the superscript s refers to a 3D solid element. M, K, and C are the mass, stiffness, and damping of the 3D solid element, the specific forms of which are also given in the Appendix A.

The displacement conditions, temperature conditions, and interface force balance conditions between the extended layer-by-layer and 3D solid element interface are simplified; the insulation layer and the support layer are considered rigidly connected, and their specific expressions are shown in equations (3-11).

$$\begin{array}{c}{\mathit{u}}_1^{\mathrm{p}}={\mathit{u}}_1^{\mathrm{s}},{\mathit{T}}_1^{\mathrm{s}}+{\mathit{T}}_1^{\mathrm{p}}=0,\\ {\mathit{\theta }}_1^{\mathrm{p}}={\mathit{\theta }}_1^{\mathrm{s}},{\mathit{Q}}_1^{\mathrm{p}}+{\mathit{Q}}_1^{\mathrm{s}}=0\end{array}$$

(3-11)

A schematic of the method is shown in Figure 8:

We can assemble the final control equation by using the compatibility conditional equation (3-11), which is given by the following equation:

$$\left[\begin{array}{cc}{\mathit{M}}^{11}& 0\\ 0& 0\end{array}\right]\left\{\begin{array}{c}\ddot{\mathit{U}}\\ \stackrel{\mathrm{..}}{\mathit{\theta }}\end{array}\right\}+\left[\begin{array}{cc}0& 0\\ {\mathit{C}}^{21}& {\mathit{C}}^{22}\end{array}\right]\left\{\begin{array}{c}\dot{\mathit{U}}\\ \stackrel{.}{\mathit{\theta }}\end{array}\right\}+\left[\begin{array}{cc}{\mathit{K}}^{11}& {\mathit{K}}^{12}\\ {\mathit{K}}^{21}& {\mathit{K}}^{22}\end{array}\right]\left\{\begin{array}{c}\mathit{U}\\ \mathit{\theta }\end{array}\right\}=\left\{\begin{array}{c}{\mathit{F}}^1\\ {\mathit{F}}^2\end{array}\right\}$$

(3-12)

4. Temperature Distribution of Composites with Multiple Delamination Damage Events

Figure 9 presents the temperature distributions of the composite laminate electric heating structures under different delamination conditions. The analysis indicates that when delamination occurs between the support layer and the adhesive layer (i.e., Case 1 and Case 4), temperature transfer cannot occur along the thickness direction, and the temperature distribution in the thickness direction is uneven, ultimately resulting in a nonuniform thermal stress distribution, which adversely affects structural integrity and safety.

Similarly, when delamination occurs between the wear-resistant layer and the adhesive layer (i.e., Case 2 and Case 3), the maximum temperature reaches approximately 500 °C at t = 50 s. This is because the influence of delamination damage causes the structure in the delaminated area to behave as a thermal conduction model with a high-power heat source but without convective heat transfer, leading to a sharp increase in temperature, which is unacceptable in engineering applications. The operational temperature range for the polyimide resin, the matrix material most commonly used in aircraft composite structures, is between 288 °C and 316 °C [29], at which point the member is already in a state of failure; therefore, no further calculations are carried out for the later time steps.

Furthermore, a comparison of Case 2 and Case 3 reveals that with the expansion of the delamination area, the maximum temperature increases from 500 °C (Case 3) to 550 °C (Case 2), and the low-temperature region on the outer surface further expands. This phenomenon occurs because the expansion of the delaminated area obstructs temperature transfer along the thickness, further impeding the diffusion of heat from the heating zone. Consequently, the larger the delamination area is, the higher the maximum temperature and the lower the average temperature of the anti-icing surface, which ultimately leads to premature failure of the component and low anti-icing efficiency.

Figure 10 shows the temperature distribution on the anti-icing surface under different delamination conditions when t = 300 s and Y = 0.9144 m. A comparison of Case 1 and Case 3 indicates that, with the same delamination area, different locations of delamination lead to different anti-icing efficiencies for the anti-icing surface of the electrically heated composite material structure. Moreover, a comparison of Case 2 and Case 3 reveals that when the delamination region is between the electric heating layer and the anti-icing surface, the impact of delamination significantly affects the anti-icing efficiency. Notably, a marked increase in temperature is observed in front of and behind the low-temperature area of the anti-icing surface, which is due to the impact of delamination damage, leading to a sharp increase in the temperature of the delamination area below the surface and an increase in the temperature of the anti-icing surface around the delamination area.

A comparison of the results with those of the undamaged case reveals that for Case 1 and Case 4, when delamination damage occurs below the electric heating layer, the delamination area does not affect the temperature distribution on the anti-icing surface, and the temperature distribution curves almost overlap. This phenomenon occurs because delamination does not affect the upwards transfer of heat in this configuration; therefore, the heating effect of the electrically heated area on the surface is not affected.

To react to realistic situations with multiple delamination damage events, the multiple delamination damage events shown in Figure 11 are selected for this study, and the anti-icing performance and the thermal performance are compared with those of the undamaged case. Figure 12 presents cloud diagrams of the temperature distribution in the composite laminate electric heating structure with multiple instances of delamination damage.

Because delamination damage occurs both above and below the electrically heated layer at the same time, the upwards and downwards transfer of temperature is obstructed. Consequently, the structure exhibits characteristics of both Case 2 and Case 4, with low-temperature regions on the anti-icing surface, a steep increase in temperature in the local area, and the lower part of the delamination area is almost not heated. As shown in Figure 12, the maximum temperature of the structure reaches 700 °C at t = 25 s. Compared with that shown in Figure 10, the maximum temperature in the case of multiple damage events is much greater than that in the case of a single damage event. This phenomenon arises because the damage interrupts the thermal transfer pathways in both directions. Moreover, the regions between the two delamination damage zones consist of an electric heating layer and an adhesive layer, where the electric heating layer is relatively thin and the thermal conductivity of the adhesive layer is very low. This configuration results in the blockage of heat transfer and causes severe local heating.

5. Surface Temperature and Thermal Stress Distribution of Airfoils with Preset Damage

In this paper, on the basis of Li [28] et al., the existing thermal expansion layer-by-layer/3D solid element method is applied to study the electrically heated anti-icing process of a wing with an NACA0012 airfoil. The effects of delamination damage on the anti-icing performance and thermal deformation of a wing's electroheating structure are analysed.

Figure 13 presents the geometric dimensions and boundary conditions of the wing leading-edge model utilized in this study. To simulate a wing fixed by rivets, the solid support boundary condition is set on the right side of the wing, constraining displacements in all three directions to zero. The corresponding constant heat source is applied at the location of the heating layer shown in Figure 13. To simulate the atmospheric icing environment in which the electrically heated structure of the composite laminate is located, the environmental conditions, as well as the mechanical and thermophysical parameters of each layer, remain consistent with those described earlier. Both the environmental temperature and the initial temperature are set to 251.6 K.

To compare the effects of the location of the delamination damage and the size of the delamination damage on the anti-icing performance of the electrically heated wing structure, two delamination damage locations and two delamination damage areas were considered, as shown in Figure 14.

An analysis of Figure 15 shows that when delamination damage occurs between the wear-resistant layer and the adhesive layer, as in Case 1 and Case 4, a phenomenon similar to that observed in earlier flat composite materials is evident. In this case, the temperature in the delaminated region increases sharply. Furthermore, compared with the maximum temperature of 550 °C recorded in the previous laminated anti-icing system, the maximum temperature of the wing's electrothermal anti-icing system, which has the same delamination damage location and area, decreases to approximately 450 °C. This finding indicates that when delamination damage occurs above the electric heating layer, its impact on the wing structure is less significant than its effect on flat structures. At the same time, in this case, a high-temperature region with temperatures ranging from approximately 200 °C to 450 °C emerges around the damage. If the temperature in this region can be detected in time, further deterioration of the situation can be prevented by manually switching off the electric heating anti-icing system. A comparison of Cases 2 and 3 reveals that when delamination damage occurs between the support layer and the adhesive layer, the temperature distribution differs significantly from that observed in flat structures. Notably, the support layer and thermal insulation layer surrounding the delaminated region of the wing do not exhibit the low-temperature zones observed in the flat configuration. Compared with that under undamaged conditions, the temperature only decreases by approximately 10 °C. This is due to the influence of the curvature factor, which causes delamination damage to have a relatively small effect on curved components.

As shown in Figure 16, a comparison of the temperature distributions for five different cases at t = 300 s and X = 0.9144 m provides significant insights. The temperature curves for the undamaged condition, Case 2, and Case 3 nearly coincide, indicating that when delamination damage occurs beneath the electric heating layer, delamination damage has less influence on the surface temperature of the wing. In contrast, comparisons of Cases 1 and 4 under undamaged conditions demonstrate that when delamination damage occurs above the electric heating layer, the surface temperature of the wing significantly varies, with the region near the centre point remaining almost unheated, whereas a marked increase in temperature is observed in the area surrounding the delaminated region on the wing surface. Moreover, a comparison with the composite flat structure in Figure 9 shows that when delamination occurs between the corrosion-resistant wear layer and the adhesive layer and the delamination area remains constant, the peak temperature near the low-temperature region of the wing's anti-icing surface is greater than that observed in the flat structure.

Figure 17 shows the displacement distribution in the direction perpendicular to the wing surface (u3 direction) at the position where the abscissa X is 0.9144 m after heating for 300 s for five different cases. The analysis reveals that the displacement distributions for the nondamaged case, Case 2, and Case 3 are relatively similar. However, Case 1 and Case 4 exhibit significant displacements near the centre due to thermal deformation, which is caused by the increased temperature difference resulting from delamination damage.

6. Conclusions

This study considers the heat conduction issues associated with the actual heating process, which incorporates both a constant heat source and heat flow on the anti-icing surface. This provides a supplement to the thermal expansion layer-by-layer/3D solid finite element method by deriving a new theoretical framework that accommodates the corresponding boundary conditions. The computational analysis focuses on the temperature field, displacement distribution in the axial direction, and thermal stress distribution during the electrothermal anti-icing process of preventing ice from forming on composite laminated plates with preset damage electrothermally. These findings indicate that when delamination damage occurs above the electric heating layer, the anti-icing structure in the delaminated region can be modelled as a thermal conduction case with a high-power heat source and negligible convective heat transfer, which exceeds the maximum usage temperature of the material and significantly impacts the thermodynamic performance of the structure. Conversely, when delamination damage occurs below the electric heating layer, the delamination damage has almost no effect on the temperature distribution on the anti-icing surface. Furthermore, multiple damage events have a much more pronounced effect on heat transfer obstruction, ultimately resulting in a maximum temperature that is considerably higher than that observed when considering only a single delamination event.

For aircraft wing structures with preset damage, when delamination occurs above the electric heating layer, the maximum temperature of the structure reaches 450 °C at 50 s, and the structure temperature is far greater than the material tolerance temperature. Consequently, the temperature of the anti-icing surface in the delaminated region approaches the ambient temperature, which severely reduces the anti-icing efficiency. Conversely, when delamination occurs below the electric heating layer, its impact on the temperature distribution of the anti-icing surface is negligible. Moreover, the hindrance of heat transfer due to multiple types of damage becomes more pronounced, resulting in greater detrimental effects on the structure's thermal performance.