Numerical Simulation of Ice Crystal Trajectories and Fragmentation Dynamics Around the XRF1 Fuselage Nose

Icing on airplane surfaces, probes, and in aircraft engine compressors poses a significant safety concern for the aeronautic industry [8]. It can lead to dramatic events [3] and often causes a considerable reduction in aerodynamic performance [2, 5], significant regions of unsteady flow, inconsistent probe measurements, sudden loss of engine thrust, engine flameout, and even irreversible damage. In fact, several aircraft incidents and dramatic accidents that pertain to ice accretion have been recorded since the onset of aviation [1]. Water droplets that freeze on solid airplane surfaces typically cause icing, but the presence of ice crystals (which eventually mix with supercooled water droplets) in clouds with high particle concentration that melt and accrete on warm surfaces of airplanes is responsible for multiple cases of icing since the 1980s — mainly in aircraft engines and on heated probes. This ice crystal icing (ICI) has garnered significant research interest in recent years, and several international projects have attempted to provide the aeronautical industry with multiphysics numerical simulation tools that meet certification requirements and assess ICI on new products and future generations of aircraft engines. Given the challenge of assessing ICI in ground testing while also meeting flight conditions, the development of these tools is of significant importance — especially since no current numerical tool can successfully address the issue of icing.

Simulation Tool

We employ CEDRE, ONERA’s multiphysics computational fluid dynamics (CFD) code, in our study [4]. Our present computations utilize two solvers: (i) The CHARME solver to simulate the aerodynamic flow field and (ii) the SPARTE solver to predict the ice particle trajectories within a Lagrangian approach. In the framework of the High Altitude Ice Crystals (HAIC) and MUSIC-haic projects, two models were previously implemented in CEDRE that describe the behavior of ice crystal particles in airflow and their complex interaction with walls (fragmentation and accretion) [6, 7]. This study transpires in a sequential manner, first with CHARME and then SPARTE. All of this occurs without feedback impact from the solver in question to the previous one; we do not compute the impact of particle trajectories on the gas flow field.

Gas Solver

We utilize the CHARME solver to compute the flow field by solving the Reynolds-averaged Navier-Stokes (RANS) equations for a gas mixture of two species: air and water vapor. Specifically, we use Florian Menter’s 𝑘−𝜔𝑘−𝜔 turbulence model with shear stress transport (SST) correction and perform the time integration with an implicit first-order Euler scheme. We then carry out the flux calculation with a second-order monotonic upwind scheme for conservation law method and a Harten-Lax-van Leer contact flux scheme. A symmetry plane (XZ) reduces the computational domain to half of the geometry, after which we perform steady computations. The time step was 0.01 seconds.

Lagrangian Particle Solver

We use the SPARTE Lagrangian particle solver to compute ice crystal trajectories and exchange phenomena with the gas phase. The Lagrangian solver tracks each numerical particle individually along its path from the entrance of the domain to the exit, either by melting, impacting walls, or exiting the computational domain.

Drag, gravity, added mass, and Basset history forces all influence particle motion. Gravity and Basset forces are negligible at high particle-to-air density ratios (of the order of 1,000) and small particle sizes (less than a few hundred microns). The particle motion is therefore only governed by the drag force. Our present investigation utilizes the Ganser drag coefficient CD 𝐶𝐷𝐶𝐷, which accounts for non-spherical particles and is valid for high particle Reynolds numbers.

Three distinct phase-change processes can occur along particle trajectories: sublimation, fusion (melting), and evaporation. The particle-air interface’s diffusive mass and heat transfer phenomena (heat conduction and the enthalpy variation due to phase change) drive these three changes. A (potentially partially melted) ice crystal’s interaction with a wall is mainly influenced by particle size, impact velocity, and degree of melting. Three impact regimes can exist depending on the ratio between the particle’s kinetic energy upon impact (driving fragmentation) and its surface energy: elastic particle bouncing, inelastic particle, and particle fragmentation. In the two first regimes, the particle diameter remains unchanged. While internal cracks may form in the second regime, we assume that the particle remains intact; nevertheless, its kinetic energy is much more reduced than in the first regime. In the third regime, we assume that the particle shatters.

Simulation Setup: Flight Conditions

Our investigation considers the XRF1 (an industrial standard research test case) fuselage nose geometry, which Airbus provided in the MUSIC-haic project; it comprises a fuselage and a nose. We seek to assess the ice crystals’ impact models that were developed in the framework of MUSIC-haic on a representative configuration for anemometry purposes. Doing so allows us to assess particle concentration in the vicinity of the fuselage. In fact, engineers must know these values in order to design proper ice protection systems to prevent and remove ice accretion via anti-icing and de-icing methods.

Table 1 summarizes the flight information. We assume that the aircraft is flying at Mach number 0.78 with an angle of attack of 2.05 degrees, which corresponds to an upstream velocity of Ux=Uy =Uy = 0, and Uz=Uz=Uz= 8.54 m/s. It is moving through convective weather that involves ice crystal particles (see Table 2 for the particle characteristics). The described flight point is located within the ICI envelope (defined by regulation materials), which is where ice crystal accretion is most likely to appear (see Figure 1).

Figure 2 illustrates the computational domain, which extends to 10 times the diameter of the fuselage upstream of the nose point in the lateral and vertical directions.

Simulation Setup: Boundary Conditions

The static temperature, static pressure, and velocity are prescribed. At the outlet, the far field pressure is 25,000 pascals. We apply the no-slip condition on walls (nose) except for the fuselage, where we instead use a slip condition to speed up the computation because we are only concentrating on the fuselage’s nose. We also apply an adiabatic temperature condition on the walls (fuselage and nose), consider a symmetry plane to reduce the computational domain, and conduct the computation in a fully turbulent boundary layer. Next, we employ RANS modeling to compute the aerodynamic flow field, with SST 𝑘−𝜔𝑘−𝜔 as a turbulent model. Application of the turbulent LP1 wall function correctly estimates wall fluxes (thermal and momentum) without meshing the boundary layer’s viscous sublayer.

Discussion of Results

Computation of the airflow with the CHARME solver demonstrates that the computational domain’s maximum temperature is 266.7 kelvin (K), which is far below the melting temperature of 273.15 K. The present airflow is thus not hot enough to promote particle melting. Consequently, the remainder of this text focuses on particle trajectories and the fragmentation process.

Figure 3 depicts particle trajectories around the fuselage; the trajectory plot is based on the Lagrangian simulation. Particles that arrive from the upstream position impact the nose’s wall and windshield, and some deviate around the geometry. A particle’s ability to deviate around the fuselage wall depends on the particle response time.

Figure 4 plots the contour of the particles’ mass flux on the fuselage, as well as the contour of particle volume fraction in the symmetry plane and in a plane that is normal to the flow direction. The first contour illustrates the distribution of particles on the fuselage, and the second one highlights the particles’ trajectories around the fuselage. The latter clearly indicates that particles are deflected in this area, yielding a depletion area where almost no particles go through. Figure 4b likewise portrays particle trajectories around the fuselage. It also highlights particle fragmentation and the remission of small-sized particles that result from fragmentation dynamics, and reveals the way in which particles move far from the wall based on their diameter.

Conclusions

Our study considers the numerical simulation of ice crystal particle trajectories and the fragmentation process upon wall surface impact. We aimed to assess MUSIC-haic—a newly coded fragmentation model—in the ONERA multiphysics CFD software, and we conducted the investigation in real flight conditions around a fuselage nose geometry. The results exhibit the particles’ deposition on the fuselage and windshield, which could lead to icing if ice crystals mix with supercooled water droplets or melted particles (on heated surfaces). The newly developed model effectively captures fragmentation, in that it accounts for the splitting process and reinjection of fragments into the flow. It can therefore evaluate the particles’ concentration in the vicinity of fuselage walls.

Moussa Diop delivered a minisymposium presentation on this research at the 2023 SIAM Conference on Computational Science and Engineering, which took place in Amsterdam, the Netherlands, earlier this year.

Acknowledgments: The authors gratefully acknowledge funding from European Union’s Horizon 2020 research and innovation program under agreement No 767560.

References

[1] Cao, Y., Tan, W., & Wu, Z. (2018). Aircraft icing: An ongoing threat to aviation safety. Aerosp. Sci. Tech., 75, 353-385.
[2] Hann, R., Wenz, A., Gryte, K., & Johansen, T. A. (2017). Impact of atmospheric icing on UAV aerodynamic performance. In 2017 workshop on research, education and development of unmanned aerial systems (RED-UAS) (pp. 66-71). Linköping, Sweden: Institute of Electrical and Electronics Engineers.
[3] National Transportation Safety Board. (1996). In-flight icing encounter and loss of control Simmons Airlines, d.b.a. American Eagle Flight 4184 Avions de Transport Regional (ATR) model 72-212, N401AM, Roselawn, Indiana, October 31, 1994. (Aircraft accident report). National Transportation Safety Board. Retrieved from https://www.ntsb.gov/investigations/AccidentReports/Reports/AAR9602.pdf.
[4] Refloch, A., Courbet, B., Murrone, A., Villedieu, P., Laurent, C., Gilbank, P., … Vuillot, F. (2011). CEDRE software. Aerosp. Lab, 2, 1-10.
[5] Siquig, R.A. (1990). Impact of icing on unmanned aerial vehicle (UAV) operations. (AD-A231 191). Naval Environmental Prediction Research Facility. Retrieved from https://apps.dtic.mil/sti/pdfs/ADA231191.pdf.
[6] Soubrié, T., Senoner, J.-M., Villedieu, P., & Neuteboom, M.O. (2022). Ice crystal trajectory simulations in the ICE-MACR compressor rig: Fragmentation and melting dynamics. In 2022 AIAA aviation and aeronautics forum and exposition. Chicago, IL: American Institute of Aeronautics and Astronautics.
[7] Villedieu, P., Trontin, P., & Chauvin, R. (2014). Glaciated and mixed phase ice accretion modeling using ONERA 2D icing suite. In 6th AIAA atmospheric and space environments conference. Atlanta, GA: American Institute of Aeronautics and Astronautics.
[8] Yamazaki, M., Jemcov, A., & Sakaue, H. (2021). A review on the current status of icing physics and mitigation in aviation. Aerospace, 8(7), 188.