Development of a Computational Fluid Dynamics Model for Ice Formation: Validation and Parameter Analysis

1. Introduction

Ice formation is a serious hazard for civil aviation aircraft. Ice can represent one of the main obstacles for flight safety as reported by several investigation [1,2,3,4]. The main outcomes of ice formation on the external parts are associated with performance reduction, risk of flow separation and stall, loss of authority on the control surfaces and blocking of the probe inlets.

Ice formation on aircraft occurs when the airplane or the helicopter is flying through visible moisture (rain, clouds and so on) [5,6]. By varying size of droplets and Air Temperature is possible to obtain different ice types. The three main types are rime, clear and mixed, here described:

• Rime ice: it forms with air temperature ranging from -40°C to -15°C. water droplets suspended in air immediately freeze after impact; due to the rapidity air its trapped inside ice resulting a white ice crystalline and brittle. It is easy to remove with de-icing system such as inflatable boots. Rime ice can affect seriously the aerodynamic performance due to the irregular horn shaped protrusion that affect the adhesion of the boundary layer of the airstream [7,8,9,10].
• Clear Ice: occurs at higher temperatures, ranging from -10°c to 0°C with larger supercooled water droplets. Water remaining liquid runs back as a thin film and freezes progressively. Ice formed have no air cavities and so the final ice appears translucent. Clear ice is considered the most tenacious to remove and the most critical for balancing due to the high density [11,12].
• Mixed Ice: between -15°C and -10°C mixed ice forms. It is basically a blend of the previous two with the worst characteristics of both; glaze ice is surrounded by thin feather-shaped rime ice formations [13].

Figure 1. Ice accretion on wing leading edge.

Figure 1. Ice accretion on wing leading edge.

Ice accumulation on an aircraft can led to several adverse effects, including [5,14,15]:

• Increased weight: Ice accumulation on an aircsraft's wings, tail, and other surfaces can increase its weight, which can affect its performance and fuel consumption;
• Reduced lift: Ice can change the shape of an aircraft's wings and reduce its ability to generate lift, which reduce the performance and make it harder to take off, climb, and maintain altitude.
• Increased drag: Ice can increase the drag on an aircraft's wings and other surfaces, which can reduce its speed and fuel efficiency.
• Risk of stall: the protrusions created by the ice formation can led to detachment of the boundary layer and so can induce a stall phenomenon. Furthermore stall warning system are designed to operate with clean airfoil. The profile change due to ice accumulation anticipates the stall effect without giving to pilot prodromal advices;
• Loss of control: Ice can affect an aircraft's control surfaces, such as the ailerons, elevators, and rudder, making it more difficult to maneuver and potentially leading to a loss of control.
• Engine problems: Ice can accumulate on an aircraft's engines, disrupting their airflow and potentially causing them to malfunction or stall.
• Reduced visibility/ Air Data corruption: Ice accumulation on an aircraft's windshield, windows, and air data sensors can reduce visibility and corrupt air data (pressure altitude, air speed, vertical speed) reducing the safety of the flight.

In order to permit safe flight into icing conditions anti-ice and de-ice system have been developed [16]. Deicing equipment are designed to remove formed ice using different methods; among the others the more commonly installed are pneumatically inflation of boots, thermal heating or by induced vibration in the leading edge. Deicing systems operates with a duty cycle automatically controlled. The low power required makes this equipment ideal for aircraft with low excess of power available like General Aviation or Commuters. For example surface deformation de-icing are the widely installed. These systems can use pneumatic boots [17] glued on the leading edge or Electro Magnetic Impulses (EIDI) [18].

Figure 2. De Icing System (a) Pneumatic Anti Icing system[16]; (b) EIDI [19].

Figure 2. De Icing System (a) Pneumatic Anti Icing system[16]; (b) EIDI [19].

Anti-icing comprises a different design philosophy: prevent icing formation instead of removing it periodically. Two approach are adopted among this family: running wet is the lowest power consuming, preheating the contact surface enough to prevent icing and full evaporative that vaporizes by contact the water droplet as soon as it gets in touch.

Different power source are used to pre heat the leading edge of the wing: electro thermal is one of the more adopted in more electric aircraft exploiting a conductive mesh embedded into the composite leading edge panel [20]. Copper meshes are widely used in helicopter blades [21] and for small surfaces such as Outside Air Temperature (OAT) probes or control surfaces. For larger surfaces to be anti-iced the widely system used is still the thermal pneumatic anti-ice. In this system, turbine compressor or heat exchangers with engine exhaust supplies hot air and then this air is delivered through pipes to the leading edges of the aircraft [22]. A schematic section of this system is reported in Figure 3. Different method are present at the same time on liner aircraft for different icing positions and devices. Table 1 reports briefly the most common technique applied in civil liners for different icing locations.

The novel solution, proposed by the authors [23,24], consist of a sandwich panel with a lattice core manufactured with laser bed powder fusion [25]. An artistic view is reported in Figure 4.

This novel solution permits to obtain several advantages:

• Thermal efficiency: air passage inside of the lattice core create a high turbulence that enhances the heat exchange. Moreover the punctual control over the lattice permits to customize the heat diffused where needed;
• System efficiency: integrating the system function inside of the structure permits an important weight saving;
• Reduction in construction and maintenance cost: the panel is recyclable and is a single piece plug and play without need of welding, joining or special gluing;
• Easy to rescale and manufacture: permits to install this system even in small scale aircraft or UAV with an enhancement in the flight safety and operability

This novel system permits to customize locally the heat flux in order to prevent icing formation or to melt ice already formed. For this reason is important to develop a tool that simulates accurately icing deposition and accretion and is able to integrate the local amount of heat required to have a running wet situation or a full evaporative solution.

2. Materials and Methods

Analytical or semi empirical methods are unlikely to capture the phenomenon of ice growth in detail. For this reason, numerical simulations have been developed capable of providing highly reliable predictions. Currently there are several codes developed by different countries:

• LEWICE (USA) [26];
• FENSAP-ICE (Canada) [27];
• CAPTA (France) [28];
• MULTIICE (Italy) [29];
• TRAJICE (UK) [30];

In this paper, a novel ice formation model will be implemented using STAR-CCM + software. The novelty of this tool states in the capability to extrapolate the amount of heat needed to prevent ice forming. The choice of the software lies in the will, in the future, to make this part of the simulation dialogue autonomously with the internal flow tool in order to have a closed loop optimization.

Ice formation phenomenon is not stationary and depends mainly on the motion field, the collision of the aqueous particles, the boundary layer and the resolution of thermodynamic equations that govern solidification.

The amount of mass that solidifies can be expressed as the integral in time of the particles dispersed in a volume that impact on a surface A at a velocity v. The quantity of particles per unit of volume are reported as w, the volumetric concentration. This integral becomes:

$$\mathrm{M}={{\displaystyle \int }}^{\text{}}\mathrm{wVA} \mathrm{dt}$$

(1)

However, this equation would require that all impacting particles solidify. This approximation, too coarse, can be corrected by inserting three terms: η_1,2,3.

Collision efficiency is represented by η₁; it has been observed that due to the different inertia, some particles in the collision trajectory do not impact, being transported by the aerodynamic current. This effect is governed by the relationship between the size of the average particle and the characteristic size of the body. Particles of small size do not interfere with large bodies but tend instead to accumulate on tapered surfaces, for example on the tip of the propeller blades.

The second efficiency η₂ is characterized by the collection efficiency. This parameter permits to take into account the particles that collide with the aerodynamic profile but do not adhere to it. Super-cooled large droplets for example collide and flooding adhere perfectly. On the other hand, some particularly dry and icy snow formations tend to rebound after collision without increasing the ice layer on the profile.

The third and last useful term is the accretion efficiency. This value reports the post-impact of the particles. Due to the different latent heat, different icing effectiveness are observed: some drops of large dimensions can give rise to fluid slides that come out of the profile thus reducing the mass of ice in formation.

The formula 1 can thus be rewritten as [31]:

$$\mathrm{M}={{\displaystyle \int }}^{\text{}}{\mathsf{\eta }}_1{\mathsf{\eta }}_2{\mathsf{\eta }}_3\mathrm{wVA} \mathrm{dt}$$

(2)

Ice growing, rime and glaze is reported in Figure 5

The mathematical exposition that follows is related to Glaze Ice due to its higher risk.

Solidification of liquid film due to balance of energy between air and drop of water can be expressed as:

$${\mathrm{Q}}_{\mathrm{lat}}{+\mathrm{Q}}_{\mathrm{kin}}{+\mathrm{Q}}_{\mathrm{aero}}{=\mathrm{Q}}_{\mathrm{evap}}{+\mathrm{Q}}_{\mathrm{con}}{+\mathrm{Q}}_{\mathrm{sens}}$$

(3)

Q_lat reports the latent heat of the freezing particles of water, and can be expressed by:

$${\mathrm{Q}}_{\mathrm{lat}}=(1-\mathsf{\lambda }){\dot{\mathrm{m}}}_{\mathrm{w}}{\mathrm{L}}_{\mathrm{f}}$$

(4)

The parameter λ represent the liquid fraction, L_f the specific latent heat and ${\dot{m}}_w$ is the mass flow of the liquid portion of the air stream [32].

The second term reported in Equation 3 represent the kinetic energy lost during the impact. It is expressed as follow:

$${\mathrm{Q}}_{\mathrm{kin}}=\frac{{\dot{\mathrm{m}}}_{\mathrm{w}}{\mathrm{V}}^2}{2}$$

(4)

Third term is evaluated as Equation 4.

$${\mathrm{Q}}_{\mathrm{aero}}=\frac{{\mathrm{hRv}}^2}{{2\mathrm{c}}_{\mathrm{p}-\mathrm{a}}}$$

(4)

Last term evidences the heat lostdue to evaporation and it is given by:

$${\mathrm{Q}}_{\mathrm{evap}}=\frac{{\mathrm{hkL}}_{\mathrm{vap}}\left({\mathrm{e}}_{\mathrm{s}}{-\mathrm{e}}_0\right)}{{\mathrm{c}}_{\mathrm{p}-\mathrm{a}}\mathrm{p}}$$

(4)

Q_con represents the heat lost due to convection with air. It is function of ambient temperature T₀ and surface temperature T_s; h is the convective heat transfer.

$${\mathrm{Q}}_{\mathrm{con}}{=\mathrm{h}(\mathrm{T}}_{\mathrm{s}}{-\mathrm{T}}_0)$$

(5)

Last equation define Q_sens and so the heat transferred because of the difference of temperature between airfoil surface and droplet. T_d, in Equation 6 is the droplet water temperature.

$${\mathrm{Q}}_{\mathrm{sens}}={\dot{\mathrm{m}}}_{\mathrm{w}}{\mathrm{c}}_{\mathrm{p}-\mathrm{w}}{(\mathrm{T}}_{\mathrm{s}}{-\mathrm{T}}_{\mathrm{d}})$$

(6)

The mass flow rate of the liquid phase. m_w is connected with the Liquid Water Content (LWC) a, macroscopic parameter determined in meteorology. It is usually expressed in g/m³ and the usual value ranges from 0 to 0.6 [6]. Median Volumetric Diameter (MVD) states the average diameter of the droplet. Most critical situations occur with MVD ranges from 1.5 to 50 μm. Larger diameter are associated with glaze ice [33]; there is a relation between low LWC to evidence large MVD, reported in Figure 6. Larger droplets increases the collection efficiency since they are less affected by the aerodynamic flow and tend to contact the airfoil.

Temperature is also a macroscopic parameter that influences the type of ice formed on airfoils. If it is around 0°C it maintains longer the liquid state of the particles and so creates a glaze ice. Lower temperatures instead freezes immediately and product soft rime, as reported by Figure 7.

The CDF model followed three step: ice accretion model, anti-ice model and de-ice model.

2.1. Ice accretion model

Ice accretion model is time dependent, unsteady and is structured into 5 phases:

• Flow field resolution: external flow is solved through Navier Stokes Equations.
• Dispersed Phase: calculation of the trajectories of liquid phase is entered using Dispersed Multiphase Model (DMP) present in STAR CCM+.
• Fluid Film Resolution: evaluates the impacting particles and determines freezing or liquefaction by solving Conjugate heat Transfer (CHT) thermodynamic balance.
• Ice thickness calculation; identify the local thickness formed in step 3.
• Morphing: modify the flow grid subtracting the solidified ice from the air domain.

Dispersed Multiphase Model uses equations reported previously for both phase. The approach is identified as “one way”: the external flow can influence the liquid phase through heat and aerodynamic resistance but the latter cannot influence the gas phase.

The resolution of the fluid film is based on enthalpy. Enthalpy for the liquid solid transition film H*_f is calculated as follows:

$${\mathrm{H}}_{\mathrm{f}}^{*}{=\mathrm{H}}_{\mathrm{f}}+\left({1-\mathsf{\alpha }}_{\mathrm{s}}^{*}\right){ \mathrm{H}}_{\mathrm{lat}}$$

(7)

Where H_f is the sensible Enthalpy and $1-{\alpha }_s^{*}$ represent the portion of transition film occupied by the liquid phase. ${\alpha }_s^{*}$Depends on the normalized temperature:

$${\mathsf{\alpha }}_{\mathrm{s}}^{*}=\{\begin{array}{c}{1 \mathrm{if} \mathrm{T}}^{*} \text{<} 0\\ \mathrm{f}\left({\mathrm{T}}^{*}\right){ \mathrm{if} 0 \text{<} \mathrm{T}}^{*} \text{<} 1\\ {0 \mathrm{if} 1 \text{<} \mathrm{T}}^{*}\end{array}$$

(8)

T* identify the normalize temperature, evaluated as follows:

$${\mathrm{T}}^{*}=\frac{{\mathrm{T}-\mathrm{T}}_{\mathrm{sol}}}{{\mathrm{T}}_{\mathrm{liq}}{-\mathrm{T}}_{\mathrm{sol}}}$$

(9)

Ice thickness growth is given by:

$$∆{\mathrm{h}}_{\mathrm{s}}^{\text{'}}=\{\begin{array}{cc}-\frac{∆\mathrm{H}\left(\mathrm{T}\right)-∆{\mathrm{H}}_{\mathrm{liq}}}{{\mathrm{H}}_{\mathrm{f}}} (\mathrm{h}+∆{\mathrm{h}}_{\mathrm{s}})& \mathrm{if} \mathrm{T} \ge { \mathrm{T}}_{\mathrm{liq}}\\ { \mathrm{h}}_{\mathrm{f}}{\mathsf{\alpha }}_{\mathrm{s}}^{*}& \mathrm{if} \mathrm{T} \le { \mathrm{T}}_{\mathrm{liq}}\end{array}$$

(10)

This model does not permit the reduction of the icing formed. Thought is impossible to simulate the deicing function but only the anti-icing as is. The

The magnitude of the ice formed reflect then in the morphing displacement for each phase:

$${\mathrm{d}}_{\mathrm{face}}=\mathrm{K}∆{\mathrm{h}}_{\mathrm{s}}{\mathrm{n}}_{\mathrm{face}}$$

(11)

K is called Solid Time Step Factor. Allows to widen if necessary the displacement magnitude obtaining an acceleration on the simulations.

2.1. Anti ice Model Setup

Anti-ice model has been implemented starting from ice accumulation model inserting some modifications. A metal solid part has been added, as reported in Figure 8 in purple.

Adding new metal part requires implementing a new heat exchange interface between metal profile and the shell with embedded the thin film model. The new schematization is reported in Figure 9.

Two local interface between 3D solid mesh and 2D shell mesh requires a perfect and stable correspondence between faces. To permit this correspondence to remain stable and persistent also during unsteady calculation the removal of the frozen film thickness has been disabled permanently.

2.2. De ice Model Setup

This model has the aim to evaluate the amount of time necessary to liquefy the ice deposited on the leading edge of the airfoil. Rather than being a deicing for in-flight operation this tool enable to evaluate the possibility to remove ice formed on the airfoil during pre-flight warm-up. The interest for this feature relies into the possibility to install the novel anti ice system on Unmanned Aerial Vehicles (UAV) operated from airfield with no supporting crew.

To simulate de-icing on ground a Volume of Fluid (VOF) has been adopted. An eulerian approach has been selected with a thermodynamic solidification and liquefaction model implemented. Material properties of liquid and solid phase are reported in Table 2.

The geometric part is comprises by a metallic airfoil with a superimposed layer of ice with pre-settled thickens. The airfoil has been divided in two parts to simulate the heated leading edge and non-heated back part.

Figure 10. De-ice model setup: (a) Heated Part; (b) Non Heated Part.

Figure 10. De-ice model setup: (a) Heated Part; (b) Non Heated Part.

The advantage of adopting VOF approach lies in the ease to the access at the concentration of the two phases. It is simple to observe the liquefaction in the front part and the enlargement to the rear half.

3. Results

This section will report the results collected from the model descripted in Section 2 Materials and Methods. In the first section the ice accretion model has been validated with NASA TP 2000-210031 Ice Accretion and Icing Effects for Modern Airfoils [10]; in the second section the result obtained for anti-icing will be discussed and in the third the result for the de-icing function.

3.1. Ice accretion model Validation

To validate the ice accretion model NASA Glenn Research Centre results were used [10]. FAR 25 appendix C conditions were considered. A Business Jet airfoil was chosen to be representative due to the similarity with the final application. The section of the profile is reported in Figure 11.

Five simulation have been used to compare the experimental evidences with the ice accretion simulations obtained by the model evaluation. This simulation differs for external airspeed, temperature, MVD Droplet, LWC and angle of attack. In this manner a uniform information on the performance of the model can be evaluated representative of different flight condition and different icing ambient. The property of the five tests are reported in Table 3.

The results collected are reported, graphically, in the graphs are showing the experimental point evaluated by NASA in red and the numerical simulation evaluated by Politecnico di Torino in Blue.

The results computed with the Ice Accretion Model are consistent with the data obtained by NASA. The lower camber line reports always a good correspondence with the experimental points in all test. The upper camber icing behavior replicates satisfactorily the experimental evidences with lower accretion time and lower angle of attack at high speed (i.e. Figure 12 a). Increasing the angle of attach it and reducing the airspeed it is possible to evidence some underestimation in Figure 12c; the estimation seems to reduce proportionally with the AOA as reported also in Figure 12d. overall, however, the model prediction are consistent with the experimental outcomes in shapes and magnitude in all the tests performed.

3.2. Anti icing Results

The model, validated in the ice formation and upgraded with the anti-icing feature has been tested with the conditions of Table 4 to evaluate if the heat flux imposed at the leading edge was sufficient for a running wet hypothesis.

The choice of temperature and speed has been imposed to match with experimental setup campaign conducted during winter 2022 on Aircraft.

Results obtained from the run are reported in Figure 13.

Figure 13 a represent the impingement rate, the particle of water that directly contact the surface. Due to the AOA of 6° the lower camber is the most affected with 0.033 kg/m²*s maximum. The same trend is confirmed by the second image, which reports the thickness of the fluid film, reaching a peak of 0.01mm in the lower camber. Figure 13c confirms that with the heat flux imposed there is no ice thickness accumulation on the leading edge, providing that the heat flux imposed is sufficient for a running wet anti-icing condition. After half of the chord the ice, tend to accumulate on both upper and lower profiles in minor entity. The last image reported evidences a high temperature on the leading edge (more than 9°C in le leading edge with a peak of 21.51 °C compared with the external flow temperature of -14.75°C. the upper profile, spanned by higher flow velocity presents lower temperature compared to the lower profile.

Figure 14 presents the upper and lower temperatures of the airfoil external profile varying the external air temperature. The curves present the same trend passing from -14°C (blue) to -6°C (green) and final to -2°C. it is interesting to evidence that the temperature curve remain unchanged in shape but only gets shifted upwards with towards higher temperatures.

3.3. De Icing Results

De-ice simulation is aimed to verify the possibility for the anti-icing system to remove formed ice. This task is intended to verify this possibility for the operation of UAV in airport without human personnel assistance. The same heat flow for anti-icing condition has been imposed with 5mm ice formed and a simulation time of 7 minutes.

Table 5. De-icing conditions applied.

Group	Parameters	Test Case
Aluminium	Static Temperature[°C]	-4
	Heat Flux [W/m2]	1*104
Ice /Water	Water Density [kg/m3]	997
	Ice Density [kg/m3]	917
	Water thermal conductivity [W/m*K]	0.62
	Ice thermal conductivity [W/m*K]	2.33
	Initial ice layer thickness[mm]	5

Table 5. De-icing conditions applied.

Group	Parameters	Test Case
Aluminium	Static Temperature[°C]	-4
	Heat Flux [W/m2]	1*104
Ice /Water	Water Density [kg/m3]	997
	Ice Density [kg/m3]	917
	Water thermal conductivity [W/m*K]	0.62
	Ice thermal conductivity [W/m*K]	2.33
	Initial ice layer thickness[mm]	5

Figure 15. De Icing Results after 7 minutes: (a) Temperature Profile; (b) Solid Volume Fraction.

Figure 15. De Icing Results after 7 minutes: (a) Temperature Profile; (b) Solid Volume Fraction.

After 420 seconds, the temperature at the leading edge is around 58°C and the ice is liquefying on the correspondent area. Gradually, also the rear part is starting to melt, as evidenced by Figure 15b. These preliminary results reported is cautelative because it does not consider the fall down of ice part once the layer attaching to the airfoil is liquefied due to vibration and gravity.

5. Conclusions

The purpose of the current study was to present a novel, validated, icing model that is suitable also for anti-icing evaluations. The investigation has identified that the model is suitable to predict quantitatively and qualitatively icing phenomena on airfoil and is expandable to the anti-icing requests. The most important limitation lies in the fact that the anti-icing requires cross validation with experimental benches that are under construction at the moment. Further research should be undertaken to validate this aspects and to examine more closely the de-icing phenomena.

6. Patents

This work is part of the development of the patent of a novel anti-icing system for aircraft use [24]