
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)11715-X
10.1016/j.heliyon.2024.e35684
e35684
Research Article
Environmental hazard assessment of forest fire sites to firefighting aircraft—part I: Canyon wind and temperature distribution
Luo Yaojing luoyj022@avic.com
lyj@st.gxu.edu.cn
ab⁎
Huang Lingcai huanglc003@avic.com
a
Shi Lei shil078@avic.com
a
Bao Guihao baogh001@avic.com
a
Dai Fei daphige@buaa.edu.cn
b
a AVIC General Huanan Aircraft Industry Co., Ltd., Guangdong, Zhuhai, 519090, China
b School of Electronic and Information Engineering, Beihang University, Beijing, 100000, China
⁎ Corresponding author. AVIC General Huanan Aircraft Industry Co., Ltd., Guangdong, Zhuhai, 519090, China. luoyj022@avic.comlyj@st.gxu.edu.cn
10 8 2024
30 8 2024
10 8 2024
10 16 e3568412 4 2024
29 7 2024
1 8 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
Wildfires have caused immense damage to the environment, property and human safety in recent years. Fortunately, the deployment of firefighting aircraft, particularly water bombers, has emerged as one of the most effective strategies for combating wildfires. However, the intricate environment of forest fire sites, characterized by thermal turbulence and canyon winds, presents a formidable flight risk for firefighting aircraft. To address this issue, this study conducted a comprehensive risk assessment of firefighting aircraft operating in a wildfire environment, analyzing the impact of thermal turbulence and canyon winds using a finite element simulation method. This approach not only bridges the research gap in the field of flight safety but also evaluates the environmental risks encountered by firefighting aircraft when entering forest fire sites. Our findings underscore thermal turbulence and canyon winds as potential hazards for aircraft flying over mountainous terrain, elucidating the effect of thermal turbulence on aircraft engine intakes and quantifying the total lift loss incurred during encounters with canyon winds featuring non-uniform airflow velocities. Specifically, thermal turbulence can induce instability and vibration in aircraft engines, whereas canyon winds can generate updrafts and downdrafts that may compromise aircraft structures or lead to lift loss. Furthermore, we cite several references to emphasize the multifaceted risks associated with the forest fire site environment, encompassing temperature gradients, thunderstorms, and air pollution. Such comprehensive wildfire data can be invaluable in assessing the flight safety of firefighting aircraft during water-dropping missions in forest fire scenarios.

Keywords

Wildfires
Firefighting aircraft
Risk assessment
Thermal turbulence
Canyon wind
==== Body
pmc1 Introduction

Nowadays, as global climate change intensifies, forest fires are becoming more frequent and severe around the world [[1], [2], [3], [4], [5]]. Unpredictable wildfires, which can be caused by lightning strikes, human-ignited fires or even hot weather, are a serious problem, threatening air quality, animal habitats and human health. However, most wildfire accidents are human-caused [[6], [7], [8], [9], [10]]. To minimize the human and environmental damage caused by forest fires, one effective way to combat forest fires is through the use of aerial firefighting resources, including firefighting aircraft [[11], [12], [13]]. However, the use of firefighting aircraft in forest fires presents significant challenges due to the high-risk environment in which they operate [[14], [15], [16]]. This paper focuses on the assessment of high-risk forest fire environments for firefighting aircraft, including thermal turbulence and canyon winds, and the effects of these environments on aerial firefighting, such as lift loss and aircraft engine intake temperature.

As global climate change becomes increasingly unstable, it annually increases the likelihood of forest fires [17,18]. It is important to examine the effects of forest fires and improve the reliability of firefighting aircraft [19]. For instance, a study analyzes the 2017 Knysna fires [20], South Africa's largest wildland-urban interface fire disaster, which utilized aerial reconnaissance, geo-located databases, satellite imagery, and site visits to examine fire spread mechanisms, firebreak efficiency, and the influence of construction materials, vegetation, and weather conditions. Furthermore, utilizing remote sensing and geographic information systems (GIS) technology, another study analyzed seven variables, including climate, topography, land use, environmental protection type, anthropogenic factors, fire defense, and fire data (severity and area), through Partial Least Squares-Path Models. This approach suggests that by contextualizing each study area, tailored measures to prevent and mitigate damage from forest fires can be devised [21]. Additionally, some researchers have explored the potential of unmanned aerial vehicles (UAVs) in extinguishing forest fires. A recent study also proposed the use of a swarm of UAVs equipped with water-dropping mechanisms, cameras, and sensors for forest firefighting. These UAVs could guide the deployment of firefighting resources effectively [22].

In particular, we analyze the complex environment of forest fires and their impact on firefighting aircraft. We focus on the distribution of flow velocities in canyon winds, the propagation properties of wildfire temperatures, and the effect of thermal turbulence on firefighting aircraft. Our findings suggest that canyon winds, in particular, can cause dangerous turbulence and wind shear that significantly affect the stability of firefighting aircraft. In addition, the simulation results show that thermal turbulence caused by the heat of the wildfire can create unpredictable and hazardous operating conditions for the engines of firefighting aircraft, which in turn increases the risk of aerial accidents during firefighting operations [[23], [24], [25], [26]].

We also highlight the risks and challenges that firefighting aircraft face when operating in high-risk environments during forest fires. Our analysis, utilizing a finite element simulation approach, provides an intuitive visualization of the impact of environmental factors, such as canyon winds and thermal turbulence, on firefighting aircraft. By analyzing simulation data, we can better prepare and equip firefighting aircraft to operate in high-risk environments, ultimately reducing risk and enhancing the effectiveness of aerial firefighting operations. In other words, the simulation results offer valuable insights into the dynamics of canyon wind and temperature propagation at forest fire sites, and have the benefit of improving forest fire management and prevention efforts, which are both significant and necessary.

The remainder of the paper is organized as follows. The simulation setup is presented in Section 2. Section 3 presents the simulation results, focusing on the velocity and temperature profiles of the flow, and compares the results in terms of various factors. Section 4 discusses the limitations and novelties of our study. Conclusions and future research directions are discussed in Section 5.

2 Forest fire site simulation setup methods

This section details the simulation and accuracy setup methods of the forest fire site, including the simulation geometric model, boundary conditions, computational fluid dynamics (CFD) model and initial values of fire characteristics and environmental factors. We also present here the simulation procedure.

2.1 Simulation conditions and parameters

In our study, two kinds of simulation geometric models were applied: (a) a two-dimensional (2-D) mountain model and (b) a three-dimensional (3-D) valley model. Mesh grid statistics are shown in Table 1. The simulation models were constructed with random mountain topography as a reference. Based on the 2-D mountain model and the 3-D valley model, the boundary conditions, initial values, and other mesh parameters of the geometric models were set according to the actual environmental conditions at the forest fire site [27]. Additionally, we adopted two physical field interfaces: the turbulence realizable k-ɛ interface and the heat transfer in fluids interface.Table 1 Mesh grid statistics for the 2-D mountain and 3-D valley models.

Table 1Message	2-D mountain model	3-D valley model	
Size	7.886 km2	0.0303 km3	
Number of elements	12 529	54 277	
Number of degree of freedom for solving	37 615	165 230	
Mesh vertices	7523	11 307	
Triangles	10 519	8846	
Quads	2010	/	
Tetrahedra	/	54 277	
Edge elements	505	417	
Vertex elements	7	10	

The forest fire simulation was carried out on a small server equipped with an Intel(R) Xeon(R) W-2225 CPU with a reference speed of 4.1 GHz and 64 GB of random access memory (RAM). The computer is also equipped with an NVIDIA RTX A4000 graphics card featuring 32 GB of video RAM. The operating system used was Windows 10.

2.1.1 Parameters of the 2-D mountain model

The 2-D mountain model was used to simulate canyon winds propagating over sloping terrain. The space surrounding the mountain profile constitutes the target analysis area. Consequently, a 2-D geometric simulation model featuring an irregular mountain profile in the x-y plane was established and meshed, as depicted in Fig. 1(a) and (b), respectively. It is important to note that we focus solely on the flow features occluded by the mountain profile, rather than the mountain profile itself. As is typical, the velocity of the flow increases as the mountain profile elevates the direction of the flow.Fig. 1 (a) 2-D mountain model for simulating the propagation of canyon winds; (b) Meshing effect for finite element analysis of canyon wind flow; (c) 3-D valley model for analyzing the propagation of wildfire temperature; (d) Meshing effect of the 3-D valley model.

Fig. 1

In particular, we briefly introduce the 2-D mountain model. The adapted 2-D mountain model in Fig. 1(a) has the advantage of reducing modeling time and improving the efficiency of data acquisition, compared to the 3-D valley model [see Fig. 1(c)]. In general, a denser mesh grid produces a more accurate simulation, although the calculation time is longer. Conversely, a less dense mesh grid results in a less accurate simulation but with a shorter calculation time. In our 2-D mountain model, the size of the effective calculation area shown in Fig. 1(b) is 7.886 km2. The complete mesh grid comprises 12 529 domain elements and 37 615 degrees of freedom for solving (plus 1 internal degree of freedom), and the mesh density at the mountain profile boundary of the model is high. The critical thermal properties of the material “Air” used in both the 2-D mountain model and the 3-D valley model are given in Fig. S1 of the Supplementary Material.

2.1.2 Parameters of the 3-D valley model

The 3-D valley model was used to simulate the temperature effect of forest fires in valley terrain [see Fig. 1(c)], featuring a surface area of 0.325 km2 and a volume of 0.0303 km3. We focus here on the temperature distribution properties of forest fires. Before the simulation begins, the 3-D valley model necessitates setting boundaries and initial conditions within the turbulence realizable k-ɛ and heat transfer in fluids interfaces, including airflow inlet/outlet velocities, material parameters, and simulation model mesh domain scales. Since the heat of a forest fire is known to be concentrated at the boundary of the burning zone, we designate a closed projected circular curve of length 1557.2 m as the source of heat along the boundary. Furthermore, the material parameters, which encompass the thermal properties of “Air” and “Soil” in the simulation models, were calibrated to accurately replicate a real-world forest fire. The mesh domain scale of the simulated model was designed to optimally represent the actual dimensions of the forest fire and its vicinity, taking into account our computer hardware capabilities (detailed in Section 2.1).

Similarly, the total cubic size of the 3-D valley model in Fig. 1(c) and (d) is about 0.0303 km3, which has x × y × z dimensions of 0.55 × 0.55 × 0.1 km3. The full mesh grid contains 54 277 domain elements and 165 230 degrees of freedom for solving, with an additional 18 131 internal degrees of freedom. The mesh is dense at the boundary of the linear heat source model. In the Supplementary Material, the critical properties of the material “Soil” used in the 3-D valley model are given in Fig. S2, and the “Air” properties are the same as in Fig. S1.

2.2 Governing equations

The governing equations for the turbulence realizable k-ɛ interface and the heat transfer in fluids interface are analytical functions. These governing equations are used to calculate the velocity and temperature profiles of the flow in our simulations.

2.2.1 Setup method for the turbulence realizable k-ɛ interface

The turbulence realizable k-ɛ interface in our simulation is a combination of three transport equations that govern the behavior of turbulent flow, which belongs to the Reynolds-Averaged Navier-Stokes (RANS) framework [28]. The realizable k-ε turbulence model is an extension of the standard k-ε model, incorporating built-in realizability conditions [29]. According to the general theory of turbulent flows, the force equilibrium equations that turbulent flows satisfy are as follows.(1) ρa(u⋅∇)u=F+∇⋅(−pn+KF)+gρa,

(2) ρa∇⋅u=0,

(3) KF=(μ+μT)[∇u+(∇u)T],

In equations (1), (2), (3), where F is the bulk force vector of the fluid; ρa is the air density; n is the unit vector; u is the flow velocity vector; p stands for the pressure; μT is the turbulent viscosity coefficient; KF represents the viscous stress tensor, and g is the acceleration vector of gravity.

In fact, the realizable k-ɛ is a turbulence flow model that is commonly used in engineering application analysis at present. It adds two additional independent variables on the basis of the general single-phase flow model, namely, turbulent kinetic energy k and turbulent dissipation coefficient ε, that is:(4) k(u⋅∇)ρa=Pk−ρaε+∇⋅[(μ+μTσk)∇k],

(5) ε(u⋅∇)ρa=C1ρaεSij−Cε2ρaε2k+νε+∇⋅[(μ+μTσε)∇ε],

Specially, equation (4) for the turbulent kinetic energy is the same as for the standard k-ε turbulence model, while equation (5) is solved for the turbulent kinetic energy dissipation rate ε,

where,(6) μT=Cμρak2ε,

(7) Pk=μT{∇u:[∇u+(∇u)T]},

with,(8) Cμ=1A0+kε6cos{13arccos[48SijSjkSki)(SijSij)3]}SijSij+ΩijΩij,

and,(9) C1=max{0.43,(η5+η)},

(10) η=k2Sij:Sijε,

(11) Sij=12[∇u+(∇u)T],

(12) Ωij=12[∇u−(∇u)T].

In equations (6), (7), (8), (9), (10), (11), (12), Pk represents the generating term for the turbulent kinetic energy k, while Cμ, C1, and η are coefficients related to the generation of the turbulent kinetic energy dissipation rate ε. Additionally, μ denotes the dynamic viscosity coefficient, ν represents the kinematic viscosity, Sij, Sjk, and Ski are the mean strain-rate tensors, and Ωij is the rotation-rate tensor. A0, σk, σε, and Cε2 are constants used in the turbulence flow realizable k-ε model, as detailed in Table S1 of the Supplementary Material.

It is evident that equations (4), (5) employ two turbulence length scales: k, which characterizes the energy of the turbulence, and ε, which quantifies the rate of dissipation of that energy. These equations solve for the local values of k and ε, as well as the scales associated with fluid velocity u and variations in turbulent kinetic energy Pk. Based on the relevant theories of the realizable turbulent k-ε interface in simulations, it is well-established that the flow characteristics of canyon wind can be described using CFD equations.

In particular, the reason why we adopt the realizable k-ε model instead of others, such as the shear stress transport (SST) and k-ω turbulence models, is that: (a) not all CFD models exhibit good convergence during the simulation process; (b) no CFD model can simulate the real fluid turbulence phenomenon with absolute accuracy (even for SST and k-ω), and the meshing accuracy of the simulation geometric model will significantly affect the simulation results. Therefore, the differences in fluid simulation effects among various turbulence models are not the key focus of our simulations; (c) in our 2-D mountain model, the airflow region we focus on is above hills, not valley winds or below hill height, so the differences in airflow simulation over hills among various turbulence models, such as negative pressure gradient, gravity, and backflow effects, are negligible.

The difference between the SST model and the k-ε model is that the SST model combines the advantages of both the k-ε and k-ω models, featuring a mixing term [(1–Km) part] for switching the specialized cross-diffusion term [30,31]. However, the SST model can switch to the k-ε model when Km = 1. Based on these reasons, we choose the realizable k-ε model over the SST or k-ω models.

Furthermore, the advantage and disadvantage of the realizable k-ε model compared to the k-ω model are that the k-ω model can indeed improve flow solution stability near strong inverse pressure gradients [32,33]. However, when the turbulent kinetic energy k in the free flow experiences slight changes, the turbulent viscosity coefficient μ changes dramatically, leading to obviously unreasonable flow calculations (the simulation process may be interrupted due to non-convergence). In contrast, the realizable k-ε model is not overly sensitive to variations in the turbulent kinetic energy of the free incoming flow and thus converges better during the simulation.

In our simulation examples, the inlet average airflow velocity is assumed to be constant and represents a fully developed flow (assuming an incompressible fluid). Meanwhile, the outlet airflow velocity, with a static pressure set at 0 Pa, is determined based on the CFD simulation settings, considering both wind speed and direction. Notably, the inflow/outflow fluid speed constitutes a crucial aspect of boundary conditions, particularly in the turbulence realizable k-ɛ interface.

For instance, in the 2-D mountain model [see Fig. 1(a)], we establish various air-inflow velocity conditions for comparison, spanning from 3 m/s to 45 m/s at 3 m/s intervals. The top edge of this model is designated as an open boundary. In contrast, the 3-D valley model [see Fig. 1(c)] features an air-inflow velocity set at 6 m/s. On the side faces, we designate open boundaries, excluding those directly related to the inflow/outflow and valley area boundaries (comprising a total of 6 surfaces).

2.2.2 Setup method for the heat transfer in fluids interface

At this interface, a comprehensive set of components for governing heat transfer equations is provided, encompassing the definition of fluid properties, inflow/outflow boundaries, heat flux, thermal insulation, and the inclusion of a linear heat source in the 3-D valley model. The temperature distribution equations defined within the fluid domains adhere to the heat convection-diffusion equation, which may incorporate additional contributions such as those from a linear heat source and boundary heat insulation [34]. These equations are as follows.(13) ∇⋅(Qc+Qr)+ρaCp(∂T∂t+u⋅∇T)=Tβp(∂p∂t+u⋅∇p)+Qv+Qo,

(14) βp=−1ρa⋅∂ρa∂T,

(15) Qv=τ:∇u,

In equations (13), (14), (15), where Qc and Qr are the heat flux through conduction and radiation respectively; Cp is the specific heat capacity at constant pressure, T is the temperature, t is the time variable, βp is the coefficient of thermal expansion, Qv and Qo are the viscous dissipation in the fluid and heat sources other than viscous dissipation, respectively, and τ is the viscous stress tensor [35].

In particular, the setting of inflow/outflow boundary conditions at the interface of heat transfer in the fluid is also required, which is treated in our simulations as the coupling of heat transfer to the fluid. At the inlet boundary of a fluid domain, the inflow boundary condition defines a heat flux that accounts for the energy normally carried by the fluid flow if the channel upstream to the inlet is defined as follows.(16) k∇T⋅n=ρaΔHfu⋅n,

(17) ΔHf=∫TusTinCpdT.

In equations (16), (17), where ΔHf is the enthalpy variation, Tus is the upstream temperature, and Tin is the inlet temperature. For the 2-D mountain model [see Fig. 1(a)], we set the airflow inlet and outlet on the right and left side of the mountain profile, respectively; the static pressure at the outlet boundary is defined as 0 Pa, and the normal inflow velocity ranges from 3 m/s to 45 m/s, with an interval of 3 m/s. For the 3-D valley model [see Fig. 1(c)], the airflow inlet and outlet are arranged on the opposite side walls of the air domain, while a linear heat source is mapped to the valley surface. Here, we set the power density of the linear heat source at 910 000 W/m and the inflow velocity at 6 m/s.

3 Simulation results

Unlike some controllable experiments, natural forest fire sites are complex situations that contain multiple harsh physical conditions, including some negative factors such as the turbulent flow of canyon winds and the high temperatures of the flames. Here, we analyze the results of full-dimensional simulations, including 2-D and 3-D visualizations of cloud maps.

3.1 Canyon wind effects for the 2-D mountain model

The stability of firefighting aircraft flying during water drops on nearby wildfires is known to be affected by canyon winds. In Fig. 2, we compared the simulation results of inflow velocities (ranging from 3 m/s to 45 m/s with an interval of 3 m/s) and the maximum airflow velocity at the area of the high-speed jet stream above the mountain profile, which is usually the most dangerous area when a firefighting aircraft is buzzing. As an example, Fig. 3 shows the visualization of airflow velocity distributions (with an inflow velocity of 18 m/s) represented by different colors with flow lines. Note that the visualization method for velocity distributions remains the same for each inflow velocity condition. Thus, despite the different inflow velocities, there is no discernible difference in the color gradients of the visualizations.Fig. 2 Dependence of the inflow on the maximum flow velocity in the 2-D mountain model.

Fig. 2

Fig. 3 Cloud map of the velocity distribution of the canyon wind flow.

Fig. 3

Fig. 4 illustrates the distribution of airflow velocity across multiple sampled lines, with an inflow velocity of 6 m/s. Note that all vertical/horizontal transversal lines are marked with capital letters. Fig. 5(a) and (b) show the curves of airflow velocity as a function of the length of the vertical and horizontal transversal lines, respectively, where we analyze the sampling data from the vertical/horizontal transverse lines in the flow velocity cloud map. In Fig. 5(a), we can see that decays commonly appear in the curves of “U→O,” “T→N,” and “S→M” from 0 km to 0.4 km of the length of the vertical transverse line, because the turbulent viscosity coefficient exists in the “Air” material, which can generate a leeward slope vortex, resulting in a sharp decrease in airflow velocity at the left bottom of the 2-D mountain model. This flow characteristic has been commonly used to improve the maneuverability of the aircraft's rudder surface and increase the stalling angle of attack [36]. Then, the curve of “V→P” increases rapidly from 0 km to 0.1 km of the length of the vertical transverse line, which indicates that the sampling data is acquired outside the mountain profile (since the mountain profile area is empty, part of the sampling line “V→P” is invalid). This suggests that the flow over the mountain is relatively turbulent and not as smooth as it might be on the right side of the mountain profile. However, it is important to note that the velocity of the flow is not the only factor to consider when evaluating the effect of canyon winds on a mountainous area while a firefighting aircraft is skimming over it.Fig. 4 Airflow velocity profiles obtained with multiple sampled lines at an inflow velocity of 6 m/s.

Fig. 4

Fig. 5 Curves depicting the variation of airflow velocity with respect to the length of (a) vertical transversal lines and (b) horizontal transversal lines.

Fig. 5

What's more, the curves of “X→R” and “W→Q” are relatively flat since they are not located in the turbulent zone of sharp changes in airflow velocity. In Fig. 5(b), it is clear that all the curves begin at 6 m/s, which is the initial airflow velocity set in our simulation. As for the sharp change of the “L→F” curve in the range of 1.5 km to 3 km, it is formed because the airflow is blocked by the mountain profile, causing a low velocity zone on the left side of the mountain profile. Note that the peak value points along the horizontal transverse lines of the curves “K→E,” “J→D,” “I→C,” “H→B,” and “G→A” gradually shift to the right. This is caused by the airflow turbulence area across the mountain profile, indicating that a banded shape high-speed airflow velocity distribution zone in the transverse direction is approximately consistent with a Gaussian distribution [37]. According to Bernoulli's principle [38], the pressure is lower where the airflow rate is higher, which is further influenced by the boundary stress from the mountain profile.

In addition, we also collected vertical and horizontal lateral sampling data from Fig. 4 at 3 m/s intervals ranging from 3 m/s to 45 m/s, as shown in Fig. 6, Fig. 7, respectively. In Fig. 6(a) and (b), it can be seen that on the right bottom boundary of the 2-D mountain model, a low airflow velocity area exists due to air viscosity and the blocking effect of the mountain profile, which causes the airflow velocity curves to be consistently lower than the preset inflow velocity from 0 km to 1.1 km (for “X→R”) and from 0 km to 0.75 km (for “W→Q”). In Fig. 6(c), initial jumps appear in the airflow velocity curves. There are two reasons for this: (a) an empty area of the mountain profile has invalidated part of the “V→P” line, so the airflow velocity values start from 0 m/s; (b) air viscosity effects are concentrated at the boundary of the mountain profile, leading to a rapid increase in airflow velocity. For the rest of the “V→P” line, from about 0.2 km to 1.6 km, we can see that the airflow velocity curves maintain a relatively stable downward trend, indicating a positive correlation between inflow velocity and the rate of airflow velocity decline. In Fig. 6(d), (e), and 6(f), it is obvious that the initial valley values of each curve are caused by the low velocity airflow zone in the left bottom corner of the mountain profile. However, the airflow velocity increases rapidly from approximately 0.4 km–0.9 km (for “U→O”), 0.4 km–1.1 km (for “T→N”), and 0.4 km–1.3 km (for “S→M”), indicating that a high-speed airflow zone (with airflow speeds exceeding the preset inflow velocity) is generated above the mountain profile.Fig. 6 Curves of the airflow velocity as a function of the vertical transverse line length.

Fig. 6

Fig. 7 Curves of the airflow velocity as a function of the horizontal transverse line length.

Fig. 7

As seen in Fig. 7(a), due to Bernoulli's principle, the 2-D mountain profile facing the wind generates a peak point of airflow velocity along the horizontal transverse line “L→F,” which emerges approximately 2.06 km from the horizontal axis of each curve, irrespective of the initial airflow velocity settings. Concurrently, Fig. 7(b) through 7(f) illustrate that as the longitudinal height of each horizontal transverse line (“K→E,” “J→D,” “I→C,” “H→B,” and “G→A”) increases, the resistance of the 2-D mountain profile to the airflow progressively diminishes, and correspondingly, the peak points of airflow velocity along the horizontal axis of each curve also gradually shift towards higher distances. Specifically, the peak points are observed at 2.23 km for “K→E” in Fig. 7(b), 2.50 km for “J→D” in Fig. 7(c), 2.79 km for “I→C” in Fig. 7(d), 2.93 km for “H→B” in Fig. 7(e), and 3.13 km for “G→A” in Fig. 7(f).

From the above analysis, there is no doubt that an instantaneous loss of lift will be experienced by the firefighting aircraft as it traverses the 2-D mountain profile (with the closer the vertical distance to the mountain peak, the greater the lift loss), as evidenced in Fig. S3 of the Supplementary Material. In accordance with the aerodynamic theory presented in Refs. [39,40], a sudden alteration in the aircraft's lift is inevitable, and the variation in the aircraft's lift can be mathematically described by the following equations:(18) ΔFR=S⋅Ca2ρa(ΔVL2−ΔVR2),

(19) {ΔVR=Vap−VaiΔVL=Vap−Vam.

In equations (18), (19), ΔFR (N) represents the lift variation of the aircraft, S (m2) denotes the frontal area, and ΔVR and ΔVL (m/s) are the relative velocities of the aircraft with respect to the air on the right side (under initial inflow velocity conditions) and the left side (under maximum airflow velocity conditions), respectively, of the 2-D mountain profile (see Fig. S3 in the Supplementary Material). Vap is the ground velocity of the aircraft prior to entering the 2-D mountain model area, Vam is the maximum velocity of the airflow relative to the 2-D mountain profile along the flight path of the aircraft, and Ca is the coefficient of air resistance. For this analysis, we set Ca = 1, S = 100 m2, Vap = 100 m/s, ρa = 1.29 kg/m3, and select the horizontal transverse line “K→E” for analyzing the values of Vam obtained from Fig. 7(b). It is noteworthy that the flow velocity is assumed to be 0 m/s outside the 2-D mountain model. Based on these parameters, we calculated ΔFR under various inflow velocity conditions, as presented in Table 2.Table 2 Variation of aircraft lift under different inflow velocities along the horizontal transverse line “K→E.”

Table 2ΔFR (N)	ΔVR (m/s)	ΔVL (m/s)	Vam (m/s)	Vai (m/s)	Vap (m/s)	
−9721	97	96.22	3.78	3	100	
−18760	94	92.44	7.56	6	
−27002	91	88.67	11.33	9	
−34681	88	84.89	15.11	12	
−41678	85	81.11	18.89	15	
−47993	82	77.33	22.67	18	
−53530	79	73.56	26.44	21	
−57585	76	69.88	30.22	24	
−62759	73	66.00	34.00	27	
−66349	70	62.22	37.78	30	
−69183	67	58.45	41.55	33	
−71414	64	54.67	45.33	36	
−72962	61	50.89	49.11	39	
−73769	58	47.12	52.88	42	
−73959	55	43.34	56.66	45	

In Table 2, all parameters are based on reasonable assumptions and encompass the results of calculations under extreme inflow velocity conditions. It is evident that as Vai increases linearly, the magnitude of |ΔFR| increases nonlinearly due to the term [(ΔVL)2 – (ΔVR)2] in equation (18), which represents a nonlinear curve. We can postulate that the minimum flight speed of a firefighting aircraft during a water-dropping firefighting mission is approximately 50 m/s. When Vai = 39 m/s, our simulation results indicate that ΔVL = 50.89 m/s, which is close to 50 m/s, suggesting that the aircraft may encounter insufficient lift (losing approximately 72 962 N) and unpredictable wind forces that severely affect flight stability. Consequently, when the speed of canyon winds exceeds 39 m/s, a firefighting aircraft should avoid approaching the forest fire site (particularly above the peaks) to prevent lift loss or airflow disturbance caused by the microburst phenomenon.

Specifically, based on the wind speed data [[41], [42], [43]] provided by the NOAA National Severe Storms Laboratory [44], we can understand that a microburst is a small, concentrated downburst that generates an outward burst of strong winds at or near the surface, with maximum wind speeds sometimes exceeding 100 mph (approximately equivalent to 161 km/h or 45 m/s). This is the rationale behind setting 45 m/s as the upper limit of the inflow velocity in the 2-D mountain model. However, owing to the utilization of a steady-state solver, the intricate details of the canyon wind velocity distribution in the 2-D mountain model and the temperature diffusion in the 3-D valley model under multi-time step characteristics cannot be simulated. Instead, only different fixed values can be set, such as the inflow velocity and the power density of the linear heat source.

3.2 Temperature propagation for the 3-D valley model

The analysis of temperature propagation for the 3-D valley model offers a comprehensive assessment of thermal behavior in complex terrain settings. By exploring the characteristics of temperature distribution within the 3-D landscape of the valley [see Fig. 1(c)], we aim to gain insights into how various factors, such as topography, mesh grid density, and canyon wind conditions, influence the thermal patterns observed in the 3-D valley model. These insights have practical implications for wildfire environmental assessment in such regions.

Here, we set a linear heat source to simulate the wildfire temperature distribution. The power density of the linear heat source is adjusted to 910 000 W/m, with a heat radius of 3 m. Meanwhile, the inlet airflow velocity is set at 6 m/s, resulting in a maximum temperature of approximately 1200.2 K, which is very close to the wildfire values reported in the literature [45]. However, due to the convective heat dissipation effect, the temperature magnitude in the windward side of some hills is lower than the presupposed initial value (293.15 K), reaching a minimum of 284.63 K.

Fig. 8(a) shows the temperature magnitudes nephogram in different colors, where the darker color indicates a closer temperature to 1200.2 K, and the lighter color represents a lower temperature, approaching 284.63 K. In practice, a forest fire burns in a ring-like form [see Figs. S4(a) and S4(b) in the Supplementary Material] [46,47], which justifies the use of a linear heat source in the 3-D valley model simulation. This simulation is calculated using the “stationary study” solver, simulating the harsh environment faced by a firefighting aircraft flying above the valley [see Figs. S4(c) and S4(d) in the Supplementary Material] [48]. As can be seen in Fig. 8(a) (depicting the same cloud map but from a different visual angle), the isosurface of the temperature distribution is influenced by the inflow propulsion and appears to extend in the positive y-axis direction, implying that most of the wildfire thermal energy is concentrated on the lee side of the valley ridge.Fig. 8 (a) Steady-state nephogram illustrating temperature magnitude in different colors (with an initial temperature set at 293.15 K); (b) Steady-state thermal expansion of air along a linear heat source, resulting from the heat of combustion generated by the wildfire. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 8

Conversely, according to equation (14), the setting of the power density on the linear heat source affects air density (ρa), resulting in a negative correlation between temperature and air density and a positive correlation between temperature and the thermal expansion coefficient (βp). In other words, the heat generated by the wildfire causes the air to expand, thereby subtly influencing the airflow velocity. Consequently, an isosurface of airflow velocity is observed along the linear heat source, reaching a maximum of 4.17 m/s, as illustrated in Fig. 8(b).

What's more, when a firefighting aircraft is flying near a forest fire, the aircraft engines (such as turboprop and turbofan engines) will inhale the heated air and even wood particles present in the smoke. Malfunctions may occur as soon as the air temperature at the engine inlet exceeds the acceptable operating range of the aircraft's engines. For instance, the typical acceptable intake temperatures for turboprop engines range from approximately 213 K–323 K, yet forest fires can heat the air, causing it to rise significantly above 323 K. Consequently, the engine's performance can be severely compromised. One of the primary issues that arises is the potential for engine knocking or detonation. Excessive heat in the inlet air can lead to uncontrolled ignition of the fuel-air mixture before the engine has completed its proper combustion cycle, resulting in destructive knocking sounds within the engine. Not only does heated air reduce engine efficiency, but it can also inflict long-term damage. Additionally, the elevated temperature of the inlet air can disrupt the efficiency of the engine compressor. Compressors are designed to operate within a specific temperature range, and their performance deteriorates when faced with excessively hot air, which can give rise to a phenomenon known as compressor stall [[49], [50], [51]], where the airflow through the engine becomes unstable. In short, the compressor stall phenomenon not only reduces engine thrust but also imposes additional stress on engine components.

4 Discussions

In addition to mechanical problems, the high temperature of the inlet air can adversely affect the electronic systems of aircraft engines. Many modern aircraft engines rely on sophisticated sensors and electronic control laws to regulate various parameters. Excessive air heated by wildfires can cause sensitive electronic components to overheat, leading to erratic turbine behavior, erroneous readings, and even system failures. To mitigate the challenges posed by wildfire-heated air, some aircraft are equipped with specialized inlet protection systems [52,53], which incorporate filters or barriers that prevent large debris and even ice particles from entering the inlet of aircraft engines. However, these aircraft engine inlet protection systems do not block heated air and smoke. In practice, firefighters operating aircraft intentionally minimize the time spent in the hottest and most adverse areas during flight missions, and they may choose flight paths that allow them to quickly enter and exit areas affected by hot air, thereby reducing prolonged exposure to excessive heat conditions. As always, even though pilots are trained to closely monitor engine parameters and react swiftly to any indications of degraded operational performance or malfunctions, the heat and smoke from fires can still accelerate the aging of aircraft engines and their components, ultimately leading to performance degradation or failure.

As a matter of fact, our current geometric models (including the 2-D mountain model and the 3-D valley model) and their design and simulation results are subject to five main limitations. First, real forest fire site data were obtained from a number of other references, rather than being collected by us, and thus a comparative analysis of experimental study data and simulation data under identical working conditions is lacking. Second, the simulation parameters, such as the inflow velocity and power density of the linear heat source, are determined based on the current forest fire site data [20,21] that can not represent the complex and multiply conditions of forest fire around the world. Third, the constructed geometrical model lacks relevant reference criteria, leading to randomness in the geometrical features of the mountain profile within the model. Fourth, the mesh density directly influences the accuracy and stringency of the simulations, thereby enlarging the error margin in our simulation results. More importantly, our simulations use a steady-state solver rather than a transient solver, and thus neglect dynamical changes in the forest fire environment, such as the varying intensity and direction of the fire, which can significantly affect the temperature distribution.

Therefore, it is crucial to conduct experimental studies and provide specific examples to complement the simulation data for ensuring their accuracy. The temperature and canyon winds in a forest fire environment can be better characterized by collecting real-time data from actual firefighting operations. Regarding the novelty of our research, which resides in the analytical method, we have the following: (a) In the 2-D mountain model, we first focus on the canyon wind velocity distribution curve across multiple transversal lines (see Fig. 4), and then estimate the lift loss of an aircraft skimming the mountain under different inflow velocity conditions; (b) In the 3-D valley model, the turbulence realizable k-ε interface and the heat transfer in fluids interface are coupled to acquire wildfire temperature propagation simulation data under the influence of canyon winds. Additionally, the thermal expansion effect of air along a linear heat source caused by the wildfire heat of combustion is also taken into account. In contrast to previous studies [11,17,23,24,43], there is currently less research focusing on the environmental risk analysis of firefighting aircraft entering forest fire sites using the same analytical method.

5 Conclusions

In our study, we have simulated the challenges faced by firefighting aircraft operating in wildfire situations, focusing in particular on the hazards posed by extreme heat and turbulence in the air currents. Through the application of finite element simulation methods, our study emphasizes the necessity for enhanced aircraft protection against the risks that arise in forest fire environments.

Our simulation examples focus on high-risk conditions encountered in natural forest fire sites, specifically investigating the effects of canyon winds and temperature profiles on firefighting aircraft during water-dropping missions. These simulations revealed two main aspects:

First, simulations utilizing a 2-D mountain model examined the effect of canyon winds on the lift force of the aircraft. By evaluating the airflow velocity profile at various inflow velocities ranging from 3 m/s to 45 m/s, with a 3 m/s interval, we observed that the turbulent dynamic viscosity coefficient μ indicates perturbations in the flow on the leeward side of the 2-D mountain model. Bernoulli's principle accounts for pressure variations, suggesting a sudden change in lift as the aircraft traverses a mountain peak under the influence of canyon winds. Notably, these lift variations, calculated across the horizontal transverse line “K→E,” reach a maximum airflow velocity of 56.66 m/s, potentially leading to a lift loss of approximately 73 959 N for the aircraft.

Second, our simulations of the temperature distribution in a 3-D valley model indicate that the temperature spread behavior in a random hill topography depends on the inflow velocity and direction settings. Thermal concentrations mainly occur on the lee side of valley ridges, suggesting that firefighting aircraft must avoid approaching such high-temperature areas.

As such, our further research is aimed at exploring the effect of lightning on firefighting aircraft in extreme fire conditions. In the near future, the development of simulation techniques and predictive models may provide crucial information to aid in the risk assessment and mission planning of firefighting aircraft.

Data availability statement

Data will be made available on request.

CRediT authorship contribution statement

Yaojing Luo: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis. Lingcai Huang: Writing – review & editing, Supervision, Project administration, Methodology, Conceptualization. Lei Shi: Writing – review & editing, Supervision, Methodology, Funding acquisition, Conceptualization. Guihao Bao: Writing – review & editing, Supervision, Formal analysis, Data curation. Fei Dai: Writing – review & editing, Validation, Supervision.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Supplementary data

The following is/are the supplementary data to this article:Multimedia component 1

Multimedia component 1

Multimedia component 2

Multimedia component 2

Multimedia component 3

Multimedia component 3

Multimedia component 4

Multimedia component 4

Multimedia component 5

Multimedia component 5

Acknowledgement

This work was supported in part by the Research and Verification of Aircraft Lightning and HIRF Protection Design Technology under Grant MJZ5-2N22 and in part by the Research on Key Technology of High Seaworthiness of Surface Aircraft under Grant MJZ5-3N21 .

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e35684.
==== Refs
References

1 Sabri Y. El Kamoun N. Forest fire detection and localization with wireless sensor networks Networked Systems. NETYS. Lecture Notes in Computer Science 2013 331 334 Berlin, Heidelberg
2 Wu F. Chu A. Che M. Zhou Bi-objective scheduling of fire engines for fighting forest fires: new optimization approaches IEEE Trans. Intell. Transport. Syst. 19 4 2013 1140 1151
3 Alexander F. Udo N. Assessment of social vulnerability to forest fire and hazardous facilities in Germany Int. J. Disaster Risk Reduc. 87 2023 03562
4 Marc G. Rupert S. Cornelius S. Increasing aridity causes larger and more severe forest fires across Europe Global Change Biol. 29 6 2022 1648 1659
5 Usha M. Dimri A.P. Sandhya F. Forest fires and climate attributes interact in central Himalayas: an overview and assessment Fire Ecology 19 1 2023 1 18
6 Jolly W.M. Cochrane M.A. Freeborn P.H. Holden Z.A. Brown T.J. Williamson G.J. Bowman D.M.J.S. Climate-induced variations in global wildfire danger from 1979 to 2013 Nat. Commun. 6 1 2013 1 11
7 Arndt N. Vacik H. Koch V. Arpaci A. Gossow H. Modeling human-caused forest fire ignition for assessing forest fire danger in Austria iForest 6 2013 315 325
8 Balch J.K. Bradley B.A. Abatzoglou J.T. Nagy R.C. Fusco E.J. Mahood A.L. Human-started wildfires expand the fire niche across the United States Proc. Natl. Acad. Sci. USA 114 11 2013 2946 2951
9 Stanton Florea Military aircraft equipped with Modular Airborne Firefighting Systems (MAFFS) mobilized to assist with wildfire suppression efforts National Interagency Fire Center 2023 [Online]. Available: https://www.nifc.gov/sites/default/files/media-news-release/nr_nifc_maffs_activation_08.03.23.pdf
10 Berčák R. Holuša J. Trombik J. Resnerová K. Hlásny T. A combination of human activity and climate drives forest fire occurrence in central europe: the case of the Czech republic Fire 7 4 2024 1 18
11 Stonesifer C. Cal B. Bayham J. Calkin D. Belval E. Firefighting aircraft: understanding current practices to shape future response in a changing world Environmental Sciences Proceedings 1 2022
12 Wildfires kill 5 in Luhansk region, firefighting aircraft engaged - Ukrainian Interior Ministry Inter: Russia & CIS Military Newswire 2020
13 Ardil C. Aerial firefighting aircraft selection with standard fuzzy sets using multiple criteria group decision making analysis International Journal of Transport and Vehicle Engineering 17 4 2023 136 145
14 Fung M.K. Schaller M.F. Hoff C.M. Katz M.E. Wright J.D. Widespread and intense wildfires at the Paleocene-Eocene boundary Geochemical Perspective Letters 10 2013
15 Miandashti N. Kamravaei S. Machichi K.I. Khadour F. Henderson L. Naseem M. Lacy P. Melenka L. Extremely high-risk air quality for recruited RCMP officers in 2016 fort McMurray wildfires: is there a health impact from short periods of intense wildfire smoke exposure Am. J. Respir. Crit. Care Med. 199 2019
16 Mehmet C. Özge I.P. Mehtap O.K. Ilker A. Suhrabuddin N. Masoud D. Saye N.C. GIS-based forest fire risk determination for Milas district, Turkey Nat. Hazards 119 2023 2299 2320
17 Liu Z. Eden J. Dieppois B. Blackett M. A global view of observed changes in fire weather extremes: uncertainties and attribution to climate change Climatic Change 173 14 2022 1 33 35811834
18 De Rigo D. Libertà G. Durrant T.H. Vivancos T.A. San-Miguel-Ayanz J. Forest Fire Danger Extremes in Europe under Climate Change: Variability and Uncertainty (Doctoral Dissertation) 2017 Publications Office of the European Union 1 72
19 Woodman S. Bearman C. Hayes P. Aviation safety and accident survivability: where is the need for aviation rescue fire fighting services greatest? Saf. Sci. 173 2024 106465
20 Natalia F.Q. Lesley G. Willem S.C. Patrick R. Ryan H. Ashton M. Armandt V.S. Richard W. Analysis of the 2017 Knysna fires disaster with emphasis on fire spread, home losses and the influence of vegetation and weather conditions: a South African case study Int. J. Disaster Risk Reduc. 88 2023
21 Rodriguez-Jimenez F. Lorenzo H. Acuna-Alonso C. Alvarez X. PLS-PM analysis of forest fires using remote sensing tools. The case of Xurés in the Transboundary Biosphere Reserve Ecol. Inf. 75 2023
22 Flores P. Pablo R. Ahmed R.L. Marco A.I. Mohammad S.A.C. Pascual Wild hopper: a heavy-duty UAV for day and night firefighting operations Heliyon 8 6 2022 e09588
23 Usui K. Iwasaki T. Yamazaki T. Ito J. Numerical simulations and trajectory analyses of local ‘karakkaze’ wind: a case that could have contributed to an aircraft accident at narita airport on 23 march 2009 SOLA 0 2022
24 Regmi G. Shrestha S. Maharjan S. Khadka A.K. Regmi R.P. Kaphle G.C. The weather hazards associated with the US-bangla aircraft accident at the tribhuvan international airport, Nepal Weather Forecast. 5 2020
25 Kohyu S. Iwao M. Kunio K. Yang K.T. A numerical study of water dump in aerial fire fighting Fire Saf. Sci. 8 2005 777 787
26 Struminska A. Filippone A. Flight performance analysis of aerial fire fighting Aeronaut. J. 1–29 2024
27 Natalia E. Viacheslav P. Viktor R. Roman F. Gennadiy R. Andrei T. Assessment of health risk of the baikal region population associated with the wildfire air pollution: approaches, modelling, digital environment Emerging Contam. 9 1 2023
28 Pawel M. Dariusz K. Experimental investigation on the performance of the prototype of aircraft Opposed-Piston engine with various values of intake pressure Energy Convers. Manag. 269 2022
29 Shih T.H. Liou W.W. Shabir A. Yang Z. Zhu J. A new k-ϵ eddy viscosity model for high Reynolds number turbulent flows Comput. Fluid 24 3 1995 227 238
30 Menter F.R. Two-equation eddy-viscosity turbulence models for engineering applications AIAA J. 32 8 2012 1598 1605
31 Syaiful Hasna N. M. S. K Tony S.U. Agus S. Maria F.S. Numerical simulation of heat transfer enhancement from tubes surface to airflow using concave delta winglet vortex generators Results in Engineering 16 2022
32 Dendy A. Rizwanul F.I.M. Nura M.M. COMPARISON OF STANDARD k-ε AND SST k-ω TURBULENCE MODEL FOR BREASTSHOT WATERWHEEL SIMULATION Journal of Mechanical Science and Engineering 7 2 2020 1 6
33 Chethan R.P. Devendra D. Bhaskar P. Nitin B. Numerical investigation of heat transfer in aircraft engine blade using k-ε and SST k-ω model IOP Conf. Ser. Mater. Sci. Eng. 1013 2021 012027
34 Sarjati S. Jnana R.K. Kamalini D. B Sree S.P. Kishanjit K.K. Turbulence modelling for depth-averaged velocity and boundary shear stress of a dense rigid grass bed open channel AQUA-Water Infrastructure, Ecosystems and Society 72 9 2023 1748 1769
35 Papanastasiou T. Georgiou G. Alexandrou A.N. Viscous Fluid Flow 2021 CRC press
36 Anthony D.G. Anya R.J. Karen M. Jonathan W.N. Marilyn J.S. Review of rotating wing dynamic stall: experiments and flow control Prog. Aero. Sci. 137 2023
37 Chen J.F. Yuan Y. Gao H. Zhou T.Y. Gaussian distribution-based modeling of cutting depth predictions of kerf profiles for ductile materials machined by abrasive waterjet Mater. Des. 227 2023
38 Kai S. Alexander T. Jan S. Leonard G. Franz D. Klaus D. A novel gripper for battery electrodes based on the Bernoulli-principle with integrated exhaust air compensation Procedia CIRP 23 2014
39 Dai J.Q. Zhao Y. Li Z.M. Sentiment-topic dynamic collaborative analysis-based public opinion mapping in aviation disaster management: a case study of the MU5735 air crash Int. J. Disaster Risk Reduc. 6 2013 315 325
40 Xu Z.F. Zhao K. Chen M. Du C. Liu M. Research on the influence of air resistance in spacecraft moment of inertia measurement J. Phys. Conf. 1634 1 2020
41 Salah S. Alsamamra H.R. Shoqeir J.H. Exploring wind speed for energy considerations in eastern Jerusalem-Palestine using machine-learning algorithms Energies 15 7 2022 2602
42 Siyavash F. Soheil R. Roozbeh P. Erfan A. Mehdi N. Exploring wind energy potential as a driver of sustainable development in the southern Coasts of Iran: the importance of wind speed statistical distribution model Sustainability 13 2021 7702
43 Wenzheng Y. Gao Y. Zhengyu Y. Xin Y. Mingxuan Z. Hanxiaoya Z. Poisson-gumbel model for wind speed threshold estimation of maximum wind speed Comput. Mater. Continua (CMC) 73 1 2022 563 576
44 The NOAA National Severe Storms Laboratory (NOAA) Severe weather 101: types of damaging winds [Online]. Available: https://www.nssl.noaa.gov/education/svrwx101/wind/types/
45 Amici S. Spiller D. Ansalone L. Miller L. Wildfires temperature estimation by complementary use of hyperspectral PRISMA and thermal (ECOSTRESS & L8) J. Geophys. Res.: Biogeosciences 127 12 2022
46 Heartbroken! You will never see these beautiful places destroyed by humans again [Online]. Available https://world.huanqiu.com/article/9CaKrnKn3Al
47 Brazilian Amazon has reached 36 771 fire spots this month, an increase of 175% over the same period of last month [Online]. Available https://www.sohu.com/a/335952575_292978
48 “Media Focus | 3.33 million+ views, “Kunlong” sweeps the screen!” (in Chinese) Accessed: July. 29, 2024. [Online]. Available: https://mp.weixin.qq.com/s/FAuz-vFai_f8omVv5OyOaA.
49 Guan D. Liu Y. Zhao D. Du J. Dong X.S.D. Experimental mode decomposition investigation on 3-stage axial flow compressor stall phenomena using aeroacoustics measurements Aero. Sci. Technol. 139 2023
50 Hutchings J. Hall C.A. In-stall compressor performance and the effects of Reynolds number J. Turbomach. 144 8 2022 081005
51 Akhlaghi M. Azizi Y. Nouri N.M. Estimations of compressor stall and surge using passage stall behaviors Machines 10 8 2022 706
52 Przybyła B.S. Przysowa R. Zapałowicz Z. Implementation of a new inlet protection system into HEMS fleet Aircraft Eng. Aero. Technol. 92 1 2020 67 79
53 Richard M. Ian R. Numerical optimisation of a helicopter engine inlet electrothermal ice protection system SAE International Journal of Advances and Current Practices in Mobility 2 1 2019 265 271
