# 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 signiﬁcant 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], signiﬁcant regions of unsteady ﬂow, inconsistent probe measurements, sudden loss of engine thrust, engine ﬂameout, 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 signiﬁcant 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 ﬂight 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 ﬂow ﬁeld 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, ﬁrst 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 ﬂow ﬁeld.

**Table 1.**Flight conditions.

### Gas Solver

We utilize the CHARME solver to compute the ﬂow ﬁeld 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 \(k-\omega\) turbulence model with shear stress transport (SST) correction and perform the time integration with an implicit ﬁrst-order Euler scheme. We then carry out the ﬂux calculation with a second-order monotonic upwind scheme for conservation law method and a Harten-Lax-van Leer contact ﬂux 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.

**Table 2.**Ice crystal characteristics.

### 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.

**Figure 1.**Location of flight conditions inside the CS-25 App. P envelope: A Federal Aviation Administration document that provides aircraft manufacturers with the necessary requirements for aircraft certification. Figure courtesy of ANDHEO.

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 inﬂuenced 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 ﬁrst 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 ﬁrst regime. In the third regime, we assume that the particle shatters.

### Simulation Setup: Flight Conditions

**Figure 2.**The computational domain extends to 10 times the diameter of the fuselage upstream of the nose point in the lateral and vertical directions. Figure courtesy of ANDHEO.

Table 1 summarizes the ﬂight information. We assume that the aircraft is ﬂying at Mach number 0.78 with an angle of attack of 2.05 degrees, which corresponds to an upstream velocity of \(\textrm{Ux} =\) 238.64 meters/second (m/s), \(\textrm{Uy =}\) 0, and \(\textrm{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 ﬂight point is located within the ICI envelope (deﬁned 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 ﬁeld 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 ﬂow ﬁeld, with SST \(k-\omega\) as a turbulent model. Application of the turbulent LP1 wall function correctly estimates wall ﬂuxes (thermal and momentum) without meshing the boundary layer’s viscous sublayer.

**Figure 3.**Particle trajectory around the fuselage. Figure courtesy of ANDHEO.

### Discussion of Results

Computation of the airﬂow 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 airﬂow 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 ﬂux 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 ﬂow direction. The ﬁrst 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 deﬂected 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.

**Figure 4.**Particle trajectories around the fuselage.

**4a.**The particles’ mass ﬂux distribution on the fuselage.

**4b.**The resulting particle fragmentation dynamics and re-emission. Figure courtesy of ANDHEO.

### 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 ﬂight 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 ﬂow. 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] Reﬂoch, 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.