Computational Icing Analysis on NASA’s SIDRM Geometry to Investigate Collection Efficiency

Computational Icing Analysis on NASA’s SIDRM Geometry to Investigate

Collection Efficiency

Eric A. Stewart

Naval Air Warfare Center Aircraft Division , Patuxent River, MD, 20670, USA

Tadas P. Bartkus

Ohio Aerospace Institute, Cleveland, OH, 44142, USA the same wall will stick to the water film layer and begin to lower the

Abstract

wall temperature through convection and evaporative cooling. Over time, the wall surface temperature can decrease enough for ice to C omputational icing analysis results were compared to experimental adhere between the wall and water film surface. Ice will continue to icing tunnel data including aerothermal ( e.g., dry air) and supercooled accumulate until aerodynamic forces cause the ice to break from the water droplet rime - ice conditions from tests conducted in early 2022 surface and shed off the su rface [ 1 ]. Previous studies in 2018 were at the NASA Icing Research Tunnel (IRT). The Simulated Inter - conducted with a NACA 0012 airfoil test article by NASA Glenn compressor Duct Research Model (SIDRM) test article was used in Research Center in their 2nd Fundamental Ice Crystal Icing Physics this study, and its geometry represents the inter - compressor duct Test at NASA Glenn’s Propulsion Systems L aboratory to better region of a turbofan engine. The test article’s purpose is to study the understand ice crystal ic ing (ICI) [2 - 4 ].

physics of supercooled water icing and ice crystal icing. This study compared three different icing codes : FENSAP - ICE (Eulerian The Simulated Inter - compressor Duct Research Model (SIDRM) test approach), LEWICE3D (Lagrangian approach ), and GlennICE article represents the inter compressor duct region of a turbofan engine.

(Lagrangian approach ) . All three icing codes were conducted on Fig ure 1 shows the outer shroud wall’s complex geometry inside the SIDRM’s complex body flow - field and compared to different inter compressor duct with inlet guide vane “struts”. The SIDRM test experimental supercooled water rime runs. The test article article represents this axisymmetric wall by use of a 2D extruded wall instrumentation (pressure taps, thermocouples, etc.) and 3D laser scans polynomial shape with a single strut junc tion at the centerline of the of final ice shapes were used to compare against the different icing test article’s 1.8 2 m ( ~ 6 ft) span. The struts on SIDRM are NACA0012 code simulations. The o verall objectives are to understand how the at 0 ° AoA with 101.6 mm ( 4 in ) chord s . SIDRM is a symmetric icing codes handle capturing collection efficiency on the complex test geometry about its centerline . The test article has heated zones on the article’s unheated surfaces. In the a erothermal cases, pressure tap main section that can be controlled independently, and the test article readings matched the CFD results, but dry air CFD underpredicted is instrumented with heat flux gauges, thermocouples, and pressure thermocouple readings. Collection effic iency results from all three taps. The ability for the test article to heat its outer - mold - line ( OML ) icing codes matched well together on the main body leading edge, surfaces is necessary when trying to simulate the warm inner surfa ces main body slope, and the st rut leading edges. Rime experimental of the turbofan inter compressor duct region. SIDRM went through a collection efficiency was calculated from the raw 3D laser scan of the series of icing physics tests in the Icing Research Tunnel (IRT) targeted ice shape using a traditional equation. Results showed different at producing icing code development data [ 5 - 6 ]. W ith well - matches to the icing codes at the main body nose, ramp, and strut characterized icing conditions these past tests will help in developing leading edges. Al l three icing codes underpredicted the final ice shape future computational icing models to study ice crystal and supercooled using a single - shot constant ice density approach, with more difficult y water droplet accretion on a complex sloped wall surface and at a strut coming from the strut leading edge ice shape due to the swept wing junction with that wall.

like flow field. NASA ’s overall goal for this effort is to dev elop computational icing tools to assist in the design and certification of engines for flight in icing conditions.

Introduction

Aircraft engine ice crystal ingestion is a growing concern that can be attributed to many engine power - loss events with possi ble failure modes of engine stall, rollback, flameout, and physical damage to the compressor blades [ 1 ] . Ice crystals ingested into the turbofan engine will deflect and break - apart coming off the fan stage. The portion of ice crystals that enter the lower stages of the compressor will face warmer than freezing temperature gradients that will melt the ice - crystal and generate a mixed - phase condition. The ice - crystals will eventually impact compressor shroud wall surfaces and form a thin Figure 1 . How SIDRM's test article geometry relates to a turbofan's inter - water film layer at the wall surface. Additional ice - crystals impacting compressor duct outer shroud wall [ 7 ].

Page 1 of 17 In addition to targeting ice crystal tunnel runs, super cooled liquid water droplet runs were conducted . Supercooled liquid icing tests provide a baseline for icing code simulation validation data on the new 3D geometry . One of the more difficult validation data to collect is collection efficiency . C ollection efficiency is the fraction of water mass flux ac cumulated on a point of the test article surface compared to the upstream freestream droplet water mass flux [ 8 ]. The impingement or collection efficiency, β, will be higher for super cooled liquid water droplets compared to ice crystals due to ice crystal s having the ability to bounce off the surface if they can not stick to a water film surface layer . More information on the SIDRM test article , the icing tests, and the icing cloud characterization tests can be found in Bartkus et al. [ 5 ].

Th is paper will cover an icing study analyzing different icing codes involving FENSAP - ICE (Eulerian approach) [ 9 ] , LEWICE3D (Lagrangian approa ch) [1 0 ] , and GlennICE (Lagrangian approach) [1 1 ] on SIDRM’s complex body flow - field and be compared to experimental icing tunnel rime runs. The test article i nstrumentation and 3D laser scans of the final ice shapes in the run will be used to validate and compare the icing code simulations. This paper is in conjunction with a paper by Bartkus et al. [ 12 ] that provides supercooled liquid icing results from the SIDRM icing tests in 2022 .

Figure 3 . Pressure tap (Pxxx) and thermocouple (Txxx and STx) locations on SIDRM’s instrumented side.

Experimental Data Analysis

A depiction of the SIDRM test article mounted on the turn table in side the IRT tunnel test section can be seen in Fig. 2. For more information on the experiment setup and measurement instrumentation see Bartkus et al. [ 5 ].

Figure 2 . SIDRM test article setup inside the NASA IRT’s test section .

Experimental Setup For aerodynamic characterization , SIDRM has a total of 64 pressure taps that are split between an upper and lower row on the test article.

Fig ure 3 shows a 2 D side view of the nominal location s of the upper and lower pressure tap rows (labeled as Pxxx) on what is referred to as the “instrumented side” . Pressu re taps are on both the upper surface Figure 4 . View of pressure tap and thermocouple locations on SIDRM’s and lower surface side s of the test article , and some are placed staggard leading edge.

at the leading - edge (LE) nose shown in Fig . 4. Note the vertical line During the a ero thermal runs , the test article was swept through angle seen near the coordinate system is not at the exact leading edge of the ° of attack changes (0, 1, 2, 3, and 4 ) and tunnel airspeed changes of 50, main body and is only the intersection of the surface panels. Also note 100, 150, 200, and 230 knots . Total temperature was targeted to remain this global coordinate system is different than the coordinate system constant during the airspeed sweep, but this ended up varying due to presented in Bartkus et al. [ 12 ] due to the needs of the icing codes.

various limitations and testing time.

Thermocouple’s T201, T3 01, and T401 are located at the leading edge.

Fig ure 3 shows the thermocouples (labeled as Txxx and STx) on the Simulated Dry Air Aerothermal Validation Cases instrumented side. There are also five additional thermocouples on the other side or “non - instrumented side, ” and they are referred to as the The computational icing code process for FENSAP - ICE, LEWI CE3D, “b ack side” thermocouples in this paper .

and GlennICE requires calculating the “ dry air ” ( no icing cloud and no humidity ) Navier - Stokes flow solution and inputting that into each Page 2 of 17 icing code. The dry air c omputational fluid dynamics (CFD) to calculate the true adiabatic wall surface temperature inside the icing simulations were conducted, and they were compared to the pressure code. The tunnel walls were simulated with an inviscid boundary tap and thermocouple measurements in SIDRM’s aerothermal runs. condition due to the assum ption that the boundary layer effect s on the An aerothermal run is different than a n icing run in the IRT, with the test article centerline icing cloud and pressure tap row measurements main difference being that no icing cloud is produced by the spray bars are negligible . Static pressures and total temperatures in the during the run. Aerothermal runs are preferred for validation because simulations matched the IRT tunnel measurements during the runs.

during an icing run all pressure taps are provided positive pressure to Aerothermal runs were recorded for 30 s ec after t he tunnel reaches prevent icing within the pressure taps. Table 1 l ists the aerothermal equilibrium and an average of that 30 s ec was used as the boundary runs that were chosen for validation. condition input values. SIDRM’s separation near the trailing edge required running a few thousand iterations until the coefficient of lift Table 1 . Aerothermal Cases Used in Validation Study and drag plat eau ed and reached a steady - state convergence.

Tunnel Air Tunnel Cobra Fig ure 6 details the slice locations on the main body analyzed for the Speed Test ID T P AoA Test Date total static CFD results for the aerodynamic pressure comparisons . Note the O centerline and x - axis of SIDRM goes through the NAC A0012 strut s . (#) ( C) (Knots) (Pa) (Deg) The maximum thickness of the strut protrudes out to ± 6.096 mm (± 2/16/2022 3 5.23 150.3 95626.7 0 0.24 in) in the z - axis “main body spanwise” direction . Slices taken at 2/16/2022 10 5.43 200.1 92854.9 0 ± 19.05 mm (± 0.75 in ) were intended to investigate the pressure field 2/16/2022 12 5.26 230.2 90825.7 0 change from the strut at the junction between the strut root and main body. Slices taken at 152.4 mm ( 6 in ) and 304.8 mm ( 12 in ) are for 4/11/2022 12 3.73 150.5 94654.8 4 investigating areas not affected by the strut and to see how the convergence of the results look spanwise. Fig s . 7 - 9 shows good The SIDRM test article geometry, including the NACA0012 101.6 mm agreement between dry air CFD predict ions ( line curves ) and all 64 (4 inch) chord struts, was meshed inside the IRT’s rectangular tunnel experimental pressure taps (symbols) for tunnel airspeed s of 150, 200, test section . SIDRM’s high flow blockage causes the preference of a and 230 knots with the test article at 0 ° angle of attack (AoA). In rectangular prism tunnel grid over a spherical far - field grid due to addition, Fig. 1 0 for the 150 knots and the test article at 4° showed the constrained air flow between the test article and tunnel walls . The inlet simulation ma tched the pressure tap readings on the pressure and and outlet of the IRT tunnel was extended about 12 body len gths from suction sides. For the plots, pressure coefficients are plotted against the the test article to ensure convergence stability and to not rely as much normalized chord length ( x represents the chordwise direction and c is on the pressure outlet value behind SIDRM. Angle of attack changes the chord length). The CFD slice data from the upper and lower were created by rotating the test article and wake planes and re - surface s of the main body are l aying right on top of each other because generating a new volume mesh. The mesh wa s composed of using a the simulation is at 0°. The upper and lower pressure tap rows are hybrid mesh of prism layers (with y+ ~1 at the fastest tunnel speed ) at identified by different symbols. For Figs . 7 - 9 at 0° AoA, t hey should the surface and an unstructured tetrahedral for the rest of the volume overlap with each other, but s ome locations are off in the non - separated grid. Grid r efinement was placed around the test article and flow regions. The simulation predict s flow separation , where the Cp downstream to capture the wake . The mesh for the 0 ° AoA can be seen flattens to the trailing edge, at x /c = 0.8 2, and these variance s in slice in Fig. 5 and shows the tunnel geometry and refined mesh around the locations show the simulation having difficulty converging the test article. Total mesh is around 1 3.8 million cells.

separat ion location in these steady - state simulations.

Figure 5 . Tunnel grid mesh and refined mesh around the SIDRM test article.

Dry air CFD simulations were r u n in ANSYS Fluent [ 13 ] using the Spalart - Allmaras Turbulence model with viscous heating and curvature correction enabled. The viscous heating option was activated for the viscous dissipation terms and to better simulate the thermal energy generated by the fluid viscous shear at the wall [ 13 ]. This viscous heating option was selected because turbulent heat fluxes are needed to capture the ice accretion. A velocity inlet and pressure outlet boundary condition were applied to the tunnel grid, and a no - slip surface and fixed surface temperature on the test article. FENSAP - ICE require s a wall thermal boundary condition of around Total Temp+10 ° K wall surface temperature , while GlennICE requires two Figure 6 . Dry air CFD pressure contour in Pa units with slice locations shown flow solutions with different thermal wall temperatures at the surface for 2 0 0 knots and 0° AoA case.

Page 3 of 17 4in Strut SIDRM 150knts and 3.7degC at 4deg AoA with Cp vs. x/c 4in Strut SIDRM 150knts and 5.2 degC with Cp vs. x/c -3 -3.5 0.75 inch Slice 0.75 inch Slice -0.75 inch Slice -0.75 inch Slice -3 -2.5 6 inch Slice 6 inch Slice 12 inch Slice 12 inch Slice -2.5 exp_lower taps -2 exp_lower taps exp_upper taps exp_upper taps -2 Strut LE and TE Strut LE and TE -1.5 -1.5 Suction -1 -1 -0.5 -0.5 0.5 0.5 Pressure Coefficient, Cp Pressure Coefficient, Cp 1 Pressure 1.5 1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x/c x/c Figure 7 . Dry air CFD to experimental Cp vs. x/c comparison at 150 knots and Figure 10 . Dry air CFD to experimental Cp vs. x/c comparison at 150 knots and 0 ° AoA. 4 ° AoA.

To get the dry air CFD surface temperature, another CFD simulation 4in Strut SIDRM 200knts and 5.4 degC with Cp vs. x/c is needed where the test article wall boundary condition is changed -3.5 from a fixed wall temperature to a zero - heat flux condition. ANSYS 0.75 inch Slice -0.75 inch Slice Fluent calculates what the adiabatic wall temperature should be, and -3 6 inch Slice that result is used to compare to the experimental therm o couple 12 inch Slice -2.5 exp_lower taps measurements . Fig ures 1 1 - 1 3 shows the trend between the prediction exp_upper taps and measurements to be in agreement between CFD predictions (line -2 Strut LE and TE curves) and all 20 main body experimental thermocouples (symbols) -1.5 for tunnel airspeeds of 150, 200, and 230 knots with the test article at 0° angle of attack (AoA). Temp erature d ifferences can be from the test -1 article not reaching equilibrium in the aerothermal run and from the -0.5 tolerance of the thermocouples themselves. In a ddition, the thermocouples were imbed ded flush with the surface, and they are likely measuring warmer temperatures due to being closer to the 0.5 Pressure Coefficient, Cp heaters on the inner - mold - line side . The 4 ° AoA case shown in Fig. 1 4 shows a ± 0.5°K agreement with the instrumen ted side and back side thermocouples to the CFD .

1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x/c Figure 8 . Dry air CFD to experimental Cp vs. x/c comparison at 200 knots and 4in Strut SIDRM 150knts and 5degC with T vs. x/c 0 ° AoA.

4in Strut SIDRM 230knts and 5.3 degC with Cp vs. x/c 278.5 -3.5 0.75 inch Slice -0.75 inch Slice -3 6 inch Slice 12 inch Slice -2.5 exp_lower taps exp_upper taps -2 Strut LE and TE 277.5 0.75 inch Slice -1.5 -0.75 inch Slice 6 inch Slice -1 12 inch Slice exp Instr. Side TCs -0.5 Surface Temperature, degK exp Back Side TCs Strut LE and TE 276.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.5 x/c Pressure Coefficient, Cp Figure 11 . Dry air CFD to experimental surface temp comparison at 150 knots and 0 ° AoA.

1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x/c Figure 9 . Dry air CFD to experimental Cp vs. x/c comparison at 230 knots and 0 ° AoA.

Page 4 of 17 4in Strut SIDRM 200knts and 5.4 degC with T vs. x/c Icing Code Methodology The icing codes use d the same dry air flow solutions from the ANSYS 278.5 Fluent run as inputs. Collection efficiency and ice shapes for supercooled liquid water droplets were calculated for SIDRM’s geometry using FENSAP - ICE , LEWICE3D , and GlennICE.

277.5 The ANSYS FENSAP - ICE v2021 R1 code uses a Eulerian approach for droplets and for solving the heat and mass balance for the control 277 volum e. FENSAP - ICE’s DROP3D and ICE3D in sequence were the 0.75 inch Slice 276.5 main modules used in this analysis with this icing code . DROP3D is -0.75 inch Slice 6 inch Slice a 3D finite element Eulerian water droplet/ice crystal impingement 12 inch Slice solver that was used to simulate different cloud droplet/ice - crystal exp Instr. Side TCs exp Back Side TCs distributions. ICE3D takes the data from DROP3D and uses a 3D 275.5 Strut LE and TE Surface Temperature, deg K finite volume water runback and ice accretion solver. It can produce water film thickness, 3D ice accretion shapes, and final surface 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 temperatures [ 9 ] . For collection efficiency it is defined a s , x/c ⃗ ⃗ 𝑳𝑾𝑪 𝑽 ∙ 𝒏 ⃗ ⃗ Figure 12 . Dry air CFD to experimental surface temp comparison at 200 knots 𝒍𝒐𝒄𝒂𝒍 𝒅 Eq - 1 𝜷 = − and 0 ° AoA.

( 𝑳𝑾𝑪 ) 𝑽 ∞ ∞ ⃗ where is the 𝑉 droplet velocity vector , 𝑉 is freestream velocity, and 𝑑 ∞ 𝑛 ⃗ is the surface normal vector. Liquid water content (LWC) local and 4in Strut SIDRM 230knts and 5.3 degC with T vs. x/c freestream are designated with 𝐿𝑊𝐶 and 𝐿𝑊𝐶 . To calculate the 𝑙𝑜𝑐𝑎𝑙 ∞ ice accretion and water runback, ICE3D uses a thin - film Shallow - Water Icing Model (SWIM), which is a thermodynamic model based on a system of partial differential equations [ 14 ] . The SWIM uses the convection heat transfer to define the cooling effect and shear stress from the flow solver to calculate the runback from the airflow solution as an input, rather than using an empirical equation. Heat transfer coefficient is derived by t he CFD turbulence model.

FENSAP - ICE utilizes native Fluent exports and does not need any processing to read in the dry air solution. The guideline for default 276 0.75 inch Slice droplet convergence used is 1e - 10 for residual change in total beta -0.75 inch Slice (β) , but some smaller size dr oplet runs oscillated in the 1e - 07 range 6 inch Slice 12 inch Slice for these SIDRM simulations .

exp Instr. Side TCs exp Back Side TCs Surface Temperature, deg K Strut LE and TE The LEWICE icing software was developed by the Icing Branch at the NASA Glenn Research Center for the analysis of the ice accretion 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 process on airframes, and the code uses a Lagrangian approach for x/c tracking individual stream tubes representing a small finite area and Figure 13 . Dry air CFD to experimental surface temp comparison at 230 knots an inviscid potential flow calculation. LEWICE3D v3.6.3 was and 0 ° AoA.

developed for 3D complex flows and can utilize the flow solutions from other CFD solver programs, but the softwa re uses a two - dimensional cut approach using the two - dimensional LEWICE 4in Strut SIDRM 150knts and 3.7degC at 4deg AoA with T vs. x/c software for calculating the ice growth at the interest area. A Monte - 277.5 Carlo approach is used for water droplet stream tubes . Here collection efficiency is calculated as the freestream area of the impinging upstream stream tubes divided by the area of the surface 277 Instr. Side cell impacted . LEWICE3D uses a single time step approach to calculating the ice profile shape. The software calculates the 276.5 thermodynamics of the freezing process of s uper - cooled droplet Pressure Suction impingement to evaluate the effects of various icing conditions [ 10 ] .

The Fluent dry air solution was imported into Tecplot, and all needed variables were converted from cell centered to node centered. Then 0.75 inch Slice Tecplot was outputted as a .s zplt file. A Tecplot to LEWICE3D -0.75 inch Slice converter was used to nondimensionalize the variables and put them 6 inch Slice 275.5 12 inch Slice into a format that LEWICE3D’s T rajgrid could read. A constant ice exp Inst. Side TCs Surface Temperature, deg K density , RHOI variable name, was set for most of the analysis , but a exp Back Side TCs Strut LE and TE swept wing ice densit y model, IDENSITYM=1, was used for one case as an additional comparison. This model calculates the ice 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x/c density for swept wing flows based off correlations to experimental data with curvature correction factors [1 5 ].

Figure 14 . Dry air CFD to experimental surface temp comparison at 150 knots and 4° AoA.

Page 5 of 17 The new er GlennICE (Glenn Icing Computational Environment) v2. 2 icing software also uses a Lagrangian approach but utilize s a newer IRT to 7-Bin 29.9 MVD Droplet Distribution adaptive refinement algorithm tracking method [ 15 - 17 ] designed for 0.9 3D viscous flow solutions. It can utilize the discretized flow field generated from the CFD flow solvers. GlennICE mainly requires 0.8 IRT PSD pressure, density, (u, v, w) velocity components, and (x, y, z) surface 0.7 Computed 7-bin wall shear component variables. For the ANSYS Fluent solutions, 0.6 GlennICE needs an insert variable of Cell Wall Distance that is used in the adaptive refinement algorithm. The heat transfer coefficient or 0.5 the wall heat flux with the surface temperature are needed for ice grow 0.4 calculations. In this analysis, temperature ( only needs surface temperature), heat transfer coefficient (ht c) , and surface heat flux were 0.3 Cumulative Fraction output from Fluent. Ice accretion is solved for the entire 3D surface, 0.2 and heat transfer coefficient can be provided from the CFD turbulence 0.1 model or calculated in GlennHT using two constant surface temperature runs [ 11 ] . GlennICE has a built in Tecplot to GlennICE 0 30 60 90 120 150 180 210 240 converter to read in the .szplt flow solution generated for LEWICE3D .

Droplet Diameter (microns) Droplets used the adaptive refinement algorithm targeting a high Figure 15 . IRT PSD for 29.9 μm (2019 Calibration) droplet distribution vs.

fraction contained tolerance on the boundary surfaces of interest to calculated 7 - bin distribution .

ensure convergence. For the ice growth, the McClain roughness model with turbulent htc augmentation of 5.0 and an ideal rime limit of 0.009 An example of collection efficiency for GlennICE for 200 knots at 29.9 were used for the final ice s hapes.

μm MVD can be seen in Fig. 1 6 . The slice locations plotted in the following figures can be seen in Fig. 1 6 . Note the STxx thermocouples Supercooled Liquid Collection Efficiency Analysis on the instrumented strut are at about y = 135.6 mm ( 5.34 in ) and y = 145.8 mm ( 5.74 in ) .

This analysis focuse s on supercooled liquid droplet rime icing runs and shows the results for three rime ice tunnel runs . The variables that define an icing condition are the cloud median volumetri c diameter (MVD), total air temperature (T ), angle of attack (AoA) , liquid total water content (LWC), relative humidity, and ice accretion time (t). For all rime icing runs, a 𝑇 of around - 17°C was targeted . The IRT is 𝑡𝑜𝑡𝑎𝑙 close to fully saturated as a recirculating tunnel and the relative humidity can be approximated as 100%. The selected icing runs chosen to run in the icing codes can be seen in Table 2.

Table 2 . SIDRM Rime Icing Runs Simulated Accretion Tunnel Air Tunnel Time Speed T AoA P MVD LWC total static Test ID O 3 (#) (min) ( C) (Knots) (Deg) (Pa) ( μ m) (g/m ) UG3548 7 -17.1 200.2 0 91292.5 30 0.45 UG3538 5 -17.0 149.8 0 94141.2 30 0.45 UG3540 5 -17.0 150.5 4 94104.7 30 0.45 For the IRT particle size distribution ( PSD ) data (from 2019 Calibration) for 29.9 μm MVD had 36 - bins [1 8 ] . To save Figure 16 . Collection efficiency contour from GlennICE for 200 knots at computational calculation for the runs , an interpolated 7 - bin 29.9μm MVD.

distribution was calculated from the data. This new droplet distribution can be seen in Fig. 1 5 as the green points vs. the red current IRT Figures 1 7 - 2 2 show surface air temperature comparison before the distribution. The c umulative f raction is interpolated from a Langmuir spray system is turned on for UG3548 at 200 knots , UG3538 at 150 D distri bution's percent c umulative LWC w eight. Only the weight knots , and UG3540 for 150 knots at 4° AoA. The recorded data average needs to sum to 100% . Table 3 contains the droplet before the cloud is around 30 sec. The surface temperature for the z - distribution information used in the icing codes.

axis slices on the main body are shown in Figs. 1 7 , 19 , and 2 1 in comparison to the experiment thermocouples . In addition, the surface Table 3 . Interpolated 7 - bin Droplet Distribution temperature for the y - axis slices on the instrumented strut are shown Drop Dia. ( μm ) Weight (%) in Figs . 1 8 , 2 0 , and 2 2 . The simulation is underpredicting by 1. 3 °K 106.7 5 on the main body and 0.7°K on the strut leading edge. These K - type thermocouples have a tolerance of ±1.1°K for this aerothermal 64.4 10 temperature range .

46.4 20 29.9 30 18.6 20 10.7 10 7.5 5 Page 6 of 17 Instrumented Strut Dry Air 30s 150knts and -17degC with T vs. Z (in) Dry Air 30s 200knts and -17degC with T vs. x/c 256.5 257.5 ST1 Y=5.34 inch ST2 Y=5.74 inch 256.4 exp Strut TCs 256.3 256.5 256.2 256.1 255.5 255.9 254.5 0.75 inch Slice -0.75 inch Slice 255.8 Surface Temperature, deg K 6 inch Slice 12 inch Slice 255.7 exp Instr. Side TCs -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Surface Temperature, deg K 253.5 exp Back Side TCs Z (in) Strut LE and TE Figure 20 . Strut predictions compared to thermocouple measurements for 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 UG3538 at 150 knots and 0 ° AoA.

x/c Figure 17 . Main body CFD predictions compared to surface temperature measurements for UG3548 at 200 knots and 0 ° AoA.

Dry Air 30s 150knts -17degC at 4deg AoA with T vs. z/c 257.5 Instrumented Strut Dry Air 30s 200knts and -17degC with T vs. Z (in) 256.6 ST1 Y=5.34 inch 256.5 ST2 Y=5.74 inch exp Strut TCs 256.4 256.5 256.3 Instr. Side 256.2 256 256.1 Pressure 255.5 Suction 255.9 0.75 inch Slice 255 -0.75 inch Slice 255.8 6 inch Slice 255.7 12 inch Slice exp middle TCs Surface Temperature, deg K 254.5 Surface Temperature, deg K 255.6 exp backside TCs Strut LE and TE 255.5 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Z (in) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x/c Figure 18 . Strut predictions compared to thermocouple measurements for UG3548 at 200 knots and 0 ° AoA.

Figure 21 . Main body CFD predictions compared to surface temperature measurements for UG3540 at 150 knots and 4 ° AoA.

Dry Air 30s 150knts and -17degC with T vs. x/c Instr. Strut Dry Air 30s 150knts 4deg -17degC with T vs. Z (in) 256.8 ST1 Y=5.34 inch 256.7 ST2 Y=5.74 inch exp Strut TCs 256.5 256.6 256.5 256.4 256.3 256.2 256.1 255.5 255.9 0.75 inch Slice -0.75 inch Slice Surface Temperature, deg K 255.8 6 inch Slice 255.7 12 inch Slice -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Surface Temperature, deg K exp Instr. Side TCs Z (in) exp Back Side TCs Strut LE and TE Figure 22 . Strut predictions compared to thermocouple measurements for 254.5 UG35 40 at 150 knots and 4 ° AoA.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x/c Figure s 23 - 2 5 show the general trends of collection efficiency Figure 19 . Main body CFD predictions compared to surface temperature produced by FENSAP - ICE for UG3548 at 200 knots . Composite measurements for UG3538 at 150 knots and 0° AoA.

curves shown refer to combin ing and weight averag ing 7 - bin runs together. Increasing in the positive Y - direction is the same going further outward on the span of the instrumented strut toward the Page 7 of 17 endcap . Collection efficiency increases with increasing y - direction in Fig. 2 3 on the instrumented strut because local velocity increases further away from the mai n body surface.

Figure 25 . FENSAP - ICE individual droplet bin beta on the instrumented strut leading edge at Y = 9 in slice for UG3548 200knts 0° AoA.

Experimental Collection Efficiency Analysis for 0 ° AoA Figure 23 . FENSAP - ICE composite beta on the instrumented strut at UG3548 After the rime icing runs, 3D laser scans are taken of the accreted ice 200knts 0° AoA.

shape profile on SIDRM’s surface. The final ice profile on the main body and struts w ere scanned after each icing test run using a 3D laser To see the trend of changing droplet sizes using FENSAP - ICE , the scanner. The scan ned area covered the center +/ - 5 - inch span , and the UG3548 200knts 0 ° AoA results are plotted in Fig s . 2 4 - 2 5 for each scanner has an accuracy range of +/ - 0.00 3”. From comparing the 3D bin and the final weight averaged “ composite ” . The beta max peak laser scans to the clean geometry, ice thickness normal to the surface refers to collection efficiency on SIDRM’s leading edge and the represented by ∆ , can be calculated. For rime cases where freezing smaller peaks “humps” are higher collection efficiency striking on fraction 𝑛 = 1 then all the impinging water would freeze on contact the steep ramp geometry. Increasing droplet size increases collection with the surface. In the se cases, rime ice thickness can be expressed as efficiency at the slope and that secondary beta peak location moves the following [ 19 ] : further back on the chord as droplet size increases. For the instrumented strut, Fig. 25 shows increasing droplet size increases 𝑳𝑾𝑪 ∙ 𝑽 ∙ 𝜷 ∙ 𝒕 ∆ = Eq - 2 beta striking the strut leading edge.

𝝆 𝒊𝒄𝒆 where liquid water content defined by 𝐿𝑊𝐶 , freestream velocity by 𝑉 , collection efficiency by 𝛽 , and icing duration by 𝑡 . The ice density will be assumed constant as 𝜌 = 917 𝑘𝑔 𝑚 ⁄ unless otherwise stated 𝑖𝑐𝑒 and is a valid assumption used in all three icing codes for strictly rime simulations. With the experimental rime ice thickness, the collect ion efficiency can be calculated by rearranging Eq.2 .

∆ ∙ 𝝆 𝒊𝒄𝒆 ( 𝜷 ) = Eq - 3 𝒆𝒙𝒑 𝑳𝑾𝑪 ∙ 𝑽 ∙ 𝒕 𝐿𝑊𝐶 will be assumed as the freestream tunnel constant value . This method for calculating collection efficiency from a rime ice shape has been performed by Tsao and Lee [20] on a swept NACA 0012 wing, and by Tsao and Porter [21] on a Common Research Model Midspan Wing Section model. However, in this test article the upstream leading edge of SIDRM will a ffect the downstream 𝐿𝑊𝐶 impacting the sloped wall. This 𝐿𝑊𝐶 effect can be seen in Fig. 2 6 from the FENSAP - ICE Figure 24 . FENSAP - ICE individual droplet bin beta on the main article result where freestream 𝐿𝑊𝐶 is green and yellow - red represent larger surface at Z = 6 in for UG3548 200knts 0° AoA.

𝐿𝑊𝐶 concentration th a n free stream . A 𝐿𝑊𝐶 higher concentration is striking the main body ramp and the near the root of the strut.

Page 8 of 17 from root to tip, and the local flow field behaves more like a swept wing as the air flow is influenced by the main body geometry .

Figure 26 . FENSAP - ICE LWC contour at UG3548 200 knots at 0° AoA for the composite to see higher concentration zones.

Even though both 𝑉 and 𝐿𝑊𝐶 can vary on the steep ramp slope, a constant freestream method was explored for Eq. 3 . In the constant Figure 28 . Dry air flow field for UG3538 150 knots at 0 ° AoA.

freestream method, a constant freestream tunnel v elo city (103 m/s for this UG3548 case) and constant freestream tunnel 𝐿𝑊𝐶 (0.45 g/m for Figure s 29 - 30 show using the method used on the experimental ice this case) were used. Fig ure 27 shows the results for back calculating thickness applied to the strut at a lower slice at Y= 6 in (152.4 mm) ⁄ experimental collection efficiency with a constant 𝜌 = 917 𝑘𝑔 𝑚 𝑖𝑐𝑒 and an upper slice at Y = 11 in (279.4 mm) . GlennICE and FENSAP - ⁄ (shown purple) or 𝜌 = 450 𝑘𝑔 𝑚 (shown orange) compared to the 𝑖𝑐𝑒 ICE predicted the same collection efficiency results on the strut, but three icing code composite beta results . All three icing codes are on LEWI CE3D predicted a lower beta max (Z = 0). Using 𝜌 = 𝑖𝑐𝑒 top of each other with only slight differences at the beta peaks. Th e 917 𝑘𝑔 𝑚 ⁄ for the exp beta created a better match to the icing codes.

exp beta curves shows that the main body leadi ng edge collection efficiency is captured well (where beta max is at Y = 0 ) with 𝜌 = 𝑖𝑐𝑒 UG3548 Sim and Experiment Composite Beta vs. Z (in) for 29.9 MVD on Instr Strut at Y = 6in ⁄ 917 𝑘𝑔 𝑚 , but the slope collection efficiency (strut root LE at Y= GlennICE ±4.94 in) is not captured when comparing to computational. Switching LEWICE3D 0.9 FENSAP-ICE Exp Beta (917 kg/m^3 Ice Density) ⁄ to 𝜌 = 450 𝑘𝑔 𝑚 we see that beta max drops significantly, but the 0.8 𝑖𝑐𝑒 Exp Beta (450 kg/m^3 Ice Density) slope β is captured better after the peak on the slope. These two ice 0.7 densities bound what the experimental ice density might be. However, 0.6 this collection efficiency capture difference could just be due to ice on 0.5 the main body nose affecting the slope collection efficiency as icing 0.4 0.3 duration time increases.

0.2 Collection Efficiency, Beta 0.1 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 Z (in) Figure 29 . Simulation and experimental beta at the instrumented strut slice Y=6 in for UG3548 200 knots rime run.

UG3548 Sim and Experiment Composite Beta vs. Z (in) for 29.9 MVD on Instr Strut at Y = 11in GlennICE LEWICE3D 0.9 FENSAP-ICE Exp Beta (917 kg/m^3 Ice Density) 0.8 Exp Beta (450 kg/m^3 Ice Density) 0.7 0.6 0.5 0.4 0.3 0.2 Collection Efficiency, Beta 0.1 Figure 27 . Simulation and experimental beta at the main body slice (Z = 4 in) 0 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 Z (in) for UG3548 200 knots rime run.

Figure 30 . Simulation and experimental beta at the instrumented strut slice The same method was applied to an instrumented strut slice at Y = 11 Y=11 in for UG3548 200 knots rime run.

in (279.4 mm) . Constant velocity and constant 𝐿𝑊𝐶 were kept to the tunnel freestream values to not change the definition of collection Figure 3 1 shows the results on a lower speed rime icing run at 150 efficiency. The local velocity is slower than tunnel frees tream near the knots, and the trends are identical to the 200 knots case. The oscillating strut leading edge due to the pressure field of the main body. Figure 2 8 scatter lines in the Exp Beta curves are due to th e data being extracted shows the velocity gradient increases going up the instrumented strut from the raw 3D laser scan data. Where there is os cillation, there are ice feathers on the ramp slope. The base true average slope ice Page 9 of 17 thickness would be in the “valleys” of the ice feather formations. Heat Transfer Coefficient Distribution at 0° AoA and 150 knots Techniques like maximum combined cross - section (MCCS) could be used to smooth out the feathers in future analysis. Results fo r heat transfer coefficient ( htc ) distribution on the main body can be seen in Fig. 34 for UG3538 150 knots run. GlennICE and FENSAP - ICE have a similar heat transfer coefficient due to both utilizing the data from the CFD solver solution. LEWICE3D uses the Integral Boundary Layer Method for calculating the heat transfer coefficient with modifications for using edge velocity, temperature, and pressure from the volume CFD solution [1 5 ]. In Fig. 34 , LEWICE3D is calculating a large dip in heat transfer coefficient at the stagnation point on the test article’s nose.

HTC vs. Y (in) for 150knts -17degC on Main Article at Z =4 in GlennICE K) LEWICE3D FENSAP-ICE Figure 31 . Simulation and experimental beta at the main body slice (Z = 4 in) 400 for UG3538 150 knots rime run.

For the strut slices, the c alculat ed exp erimental beta in Fig s . 3 2 - 3 3 , ⁄ shows the 𝜌 = 917 𝑘𝑔 𝑚 curve match ing to the icing code 𝑖𝑐𝑒 predictions improves around beta max and the extents for this 150 Heat Transfer Coefficient, W/(m knots case compared to the 200 knots case . This could just be due to -10 -8 -6 -4 -2 0 2 4 6 8 10 the icing duration for UG3548 200 knots case being 7 mins icing Y (in) duration vs. the UG3538 150 knots at 5 mins icing duration time. So, Figure 34 . Heat transfer coefficient vs. Y (in) simulation for UG3538 150 knots a better match due to the smaller ice shape on the test article.

0° AoA run at the main bod y Z = 4 in slice.

Figures 35 - 36 show the heat transfer coefficient (htc) on the instrument UG3538 Sim and Experiment Composite Beta vs. Z (in) for 29.9 MVD on Instr Strut at Y = 6in strut slices for the UG3538 150 knots run. Figure 3 5 (Y= 6 in) shows GlennICE 0.9 LEWICE3D LEWICE3D calculating a lower htc along the chord , while GlennICE FENSAP_ICE 0.8 Exp Beta (917 kg/m^3 Ice Density) and FENSAP - ICE show a sharper peak formation at the leading edge. Exp Beta (450 kg/m^3 Ice Density) 0.7 The GlennICE , LEWICE3D, and FENSAP - ICE models are using a 0.6 ⁄ constant ice density 𝜌 = 917 𝑘𝑔 𝑚 . Figure 36 (Y= 11 in) shows 𝑖𝑐𝑒 0.5 higher max values for heat transfer coefficien t than the lower root (Y 0.4 = 6 in) slice due to the higher velocity at that section cut.

0.3 0.2 Collection Efficiency, Beta HTC vs. Z(in) for 150knts -17degC on Strut at Y = 6 in 0.1 0 GlennICE -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 K) LEWICE3D Z (in) FENSAP-ICE Figure 32 . Simulation and experimental beta at the instrumented strut slice Y = 6 in for UG3538 150 knots rime run.

UG3538 Sim and Experiment Composite Beta vs. Z (in) for 29.9 MVD on Instr Strut at Y = 11in GlennICE LEWICE3D 0.9 FENSAP_ICE Exp Beta (917 kg/m^3 Ice Density) 0.8 Exp Beta (450 kg/m^3 Ice Density) 0.7 0.6 0.5 0.4 Heat Transfer Coefficient, W/(m 0.3 0 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.2 Z (in) Collection Efficiency, Beta 0.1 Figure 35 . Heat transfer coefficient vs. Z (in) simulation for UG3538 150 knots -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0° AoA run at the strut LE lower root Y = 6 in slice.

Z (in) Figure 33 . Simulation and experimental beta at the instrumented strut slice Y = 11 in for UG3538 150 knots rime run.

Page 10 of 17 HTC vs. Z(in) for 150knts -17degC on Strut at Y = 11 in 200knts Main Body Junction Ice Shape Comparison at Centerline GlennICE 6.5 LEWICE3D K) 800 FENSAP-ICE GlennICE LEWICE3D 700 FENSAP-ICE UG3548 Z=0in UG3521 Z=0in 600 5.5 Clean Y (in) 4.5 3.5 Heat Transfer Coefficient, W/(m -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Z (in) 13 13.5 14 14.5 15 15.5 16 16.5 X (in) Figure 36 . Heat transfer coefficient vs. Z (in) simulation for UG3538 150 knots 0° AoA run at the strut LE upper outboard Y = 11 in slice. Figure 38 . Computational ice shapes vs. experimental UG3548 200 knots 0 ° AoA run at the main body slope and strut junction for a 7 min icing duration .

Ice Accretion Analysis at 0 ° AoA and 200 knots Ice comparisons on the instrumented strut can be seen in Fig s . 3 9 - 40 , where LEWICE3D got closer in prediction th a n GlennICE and Results for a single - shot approach and ice density held constant at FENSAP - ICE to the experimental ice shapes . Since this is a 0° AoA ⁄ 𝜌 = 917 𝑘𝑔 𝑚 are shown in Fig. 3 7 . GlennICE and FENSAP - ICE 𝑖𝑐𝑒 case slices on the non - instrument strut (Y = - 6 in and Y = - 11 in) are did not reach the maximum thickness of the experimental ice shapes.

plotted as well for comparison . LEWICE3D overpredicts on Fig. 3 9 LEWICE3D reached the maximum thickness but created a l arge divot on the lower root slice cut and underpredicts on the Fig. 40 upp er for a rime case. LEWICE3D uses an integral boundary layer method outboard slice cut. This seems to be due to LEWICE3D relying on a on the 2D slice cut. Perhaps LEWICE3D ’s edge velocity modification 2D slice approach, and its cut is not in line with the new flow angle on the small leading edge could be causing the htc calculation to be produced by the main body pressure field.

low enough at the stagnation point to create the large divot shown.

For GlennICE and FENSAP - ICE the maximum thickness is not reached in a sing le shot approach . This could be due to needing a 200knts Rime Main Body Nose Ice Shape Comparison multi - shot approach or the computational ice density needs to decrease below 917 kg/m even though this is a fully rime case. In addition, a pointed ice shape for the experiment is interesting for this rime condition. More investigation is needed by looking at the repeat 0.5 rime icing runs to narrow down the causes for the point in the ice GlennICE LEWICE3D shape , but it could be due to the thin stru t geometry.

FENSAP-ICE UG3548 Z=-4in UG3548 Z=4in Y (in) UG3521 Z=-4in 200knts Rime Lower Strut LE Ice Shape Comparison UG3521 Z=4in Clean GlennICE LEWICE3D 0.8 -0.5 FENSAP-ICE UG3548 Y=6in UG3548 Y=-6in 0.6 UG3521 Y=6in UG3521 Y=-6in Clean 0.4 -1 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 X (in) 0.2 Z (in) Figure 37 . Computational ice shape vs. experimental UG3548 200 knots 0 ° AoA run at the main body nose for a 7 min icing duration.

-0.2 Figure 3 8 shows a centerline line slice focusing on the main body and -0.4 instrumented strut root junction, where the strut leading edge is shown as a vertical black line. For the slope ice in Fig. 3 8 all three icing codes -0.6 15.2 15.4 15.6 15.8 16 16.2 16.4 16.6 16.8 produced comparable results at the slope, but LEWICE3D formed X (in) thicker ice on the strut . The experimental laser scan results of UG3548 Figure 39 . Computational ice shape vs. experimental UG3548 and UG3521 200 and repeat run of UG3521 are also plotted. The increase in knots run on the strut LE at Y= 6 in slice for a 7 min icing duration .

experimental ice thickness might be captured better in a multi - shot simulation approach at this area or a low er ice density . As the nose leading edge ice grows then it will a ffect droplet trajectories hitting the slope ice.

Page 11 of 17 200knts Rime Upper Strut LE Ice Shape Comparison 150knts Main Body Junction Ice Shape Comparison at Centerline 6.5 GlennICE LEWICE3D 0.8 GlennICE (917 kg/m^3) FENSAP-ICE UG3548 Y=11in LEWICE3D (917 kg/m^3) UG3548 Y=-11in UG3538 Z=0in 0.6 UG3521 Y=11in GlennICE (450 kg/m^3) 5.5 UG3521 Y=-11in LEWICE3D (Swept Density) Clean 0.4 Clean 0.2 Z (in) Y (in) 4.5 -0.2 -0.4 3.5 -0.6 15.2 15.4 15.6 15.8 16 16.2 16.4 16.6 16.8 X (in) 13 13.5 14 14.5 15 15.5 16 16.5 X (in) Figure 40 . Computational ice shape vs. experimental UG3548 and UG3521 200 knots 0 ° AoA run on the strut LE at Y= 11 in slice for a 7 min icing duration .

Figure 42 . Computational ice shapes vs. experimental UG3538 150 knots 0° AoA run at the main body slope and strut junction for a 5 min icing duration.

Ice Accretion Analysis at 0° AoA and 150 knots ⁄ Fig ures 4 3 - 4 4 show that the LEWICE3D at 𝜌 = 917 𝑘𝑔 𝑚 𝑖𝑐𝑒 Ice shape comparisons for UG3538 150 knots can be seen in Fig s .

matches the experimental maximum ice thickness better , but the icing 41 - 4 4 with the icing codes matching the experimental better than the duration is shorter for this 150 knots case . When switching to 𝜌 = 𝑖𝑐𝑒 UG3538 200 knots case. In addition, Fig s . 41 - 4 4 show what happens 450 𝑘𝑔 𝑚 ⁄ , a large ice shape is predicted by LEWICE3D for the ⁄ when changing 𝜌 = 917 𝑘𝑔 𝑚 (shown solid line) to 𝜌 = 𝑖𝑐𝑒 𝑖𝑐𝑒 lower strut cut (Y = 6 in). This seems to be due to a corre ction factor 450 𝑘𝑔 𝑚 ⁄ . GlennICE and FENSAP - ICE results were laying on top inside LEWICE3D at this swept flow location. Switching GlennICE of each other for these plots with almost no deviation , so the ⁄ to 𝜌 = 450 𝑘𝑔 𝑚 improves the computational prediction to the 𝑖𝑐𝑒 FENSAP - ICE simulations were hidden to make the plots easier to maximum ice thickness, but it still does not reach the max thickness.

read for this 150 knots case. Figure 41 shows changing the ice More investigation in the future will foc us on the experimental ice density to 450 kg/m at the nose prediction increases the maximum weights to narrow down the ice density on this thin strut.

ice thickness. LEWICE3D at 𝜌 = 917 𝑘𝑔 𝑚 ⁄ gets the best match, 𝑖𝑐𝑒 but still shows a divot.

150knts Rime Lower Strut LE Ice Shape Comparison 0.8 GlennICE (917 kg/m^3) 150knts Rime Main Body Nose Ice Shape Comparison LEWICE3D (917 kg/m^3) UG3538 Y=6in 0.6 UG3538 Y=-6in GlennICE (450 kg/m^3) LEWICE3D (Swept Density) Clean 0.4 0.5 GlennICE (917 kg/m^3) 0.2 LEWICE3D (917 kg/m^3) UG3538 Z=-4in Z (in) 0 UG3538 Z=4in GlennICE (450 kg/m^3) Y (in) LEWICE3D (Swept Density) Clean -0.5 -0.2 -0.4 15.5 15.6 15.7 15.8 15.9 16 16.1 16.2 16.3 16.4 16.5 16.6 16.7 -1 X (in) -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 X (in) Figure 43 . Computational ice shape vs. experimental UG3538 150 knots 0° Figure 41 . Computational ice shape vs. experimental UG3538 150 knots 0° AoA run on the strut LE at Y= 6 in slice for a 5 min icing duration.

AoA run at the main body nose for a 5 min icing duration.

Figure 42 shows how decreasing the ice density (GlennICE) or using the swept density model (LEWICE3D) puts the icing code slope ice thickness over the experimental ice laser scan result.

Page 12 of 17 (shown purple), the main body leading edge collection efficiency is captured well (where beta max is), but the slop e collection efficiency 150knts Rime Upper Strut LE Ice Shape Comparison 0.8 is not captured well when comparing to computational predictions .

GlennICE (917 kg/m^3) ⁄ Changing to the 𝜌 = 450 𝑘𝑔 𝑚 curve (shown orange) improves LEWICE3D (917 kg/m^3) 𝑖𝑐𝑒 UG3538 Y=11in the pressure side match to experimental very well, but the beta max is 0.6 UG3538 Y=-11in significantly reduced. This m ight show that ice density is not constant GlennICE (450 kg/m^3) LEWICE3D (Swept Density) at every chord location and the slope ice density might be between the 0.4 Clean two bounds used for constant ice density.

0.2 Z (in) -0.2 -0.4 15.5 15.6 15.7 15.8 15.9 16 16.1 16.2 16.3 16.4 16.5 16.6 16.7 X (in) Figure 44 . Computational ice shape vs. experimental UG3538 150 knots 0° AoA run on the strut LE at Y= 11 in slice for a 5 min icing duration.

Experimental Collection Efficiency Analysis for 4 ° AoA The UG3540 150 knots 4° AoA rime case is interesting to see the effect of collection efficiency with angle of attack. At an angle of attack the pressure side can increase in total water catch. The LWC contour for this case is shown in Fig. 4 5 , where the pressure side is increasing in Figure 46 . Simulation and experimental beta at the main body slice (Z = 4in) LWC above freestream ( 𝐿𝑊𝐶 = 0 . 45 𝑔 𝑚 ⁄ ) and the suction side sees for UG3540 150 knots rime run.

a large shadow zone (shown blue) coming off the nose. Above the shadow zone is a high concentration zone that strikes near the lower The dry air velocity contour for this case is shown in Fig. 4 7 and strut root on the non - inst rumented strut.

shows the swept flow change between the instrumented strut on the pressure side and non - instrumented strut on the suction side. In general, the velocit y near the strut leading edge is higher on the non - instrumented strut. This higher velocity increases the collection efficiency on the non - instrumented strut compared to the instrumented strut of interest.

Figure 45 . LWC contour at 150 knots at 4° AoA for composite to see high Figure 47 . Dry air velocity flow field for UG3540 150 knots at 4 ° AoA.

concentration zones.

The result in collection efficiency changes between both struts can be Like before, a constant freestream method was explored for Eq. 3 on seen in Fig s . 4 8 - 4 9 on the lower root strut slice s (Y = 6 in and Y = - 6 the main body and is shown in Fig. 4 6 . A constant freestream tunnel in ). Beta max increases on the non - instrumented strut leading edge in velocity ( 77.1 m/s for this UG354 0 case) and constant freestream 3 Fig. 4 9 almost to 1. The experimental beta curve using 𝜌 = 𝑖𝑐𝑒 tunnel 𝐿𝑊𝐶 (0.45 g/m for this case) were used. All three icing codes 917 𝑘𝑔 𝑚 ⁄ (shown red) also gets a better prediction to the icing match each other well on the pressure side and the nose section where codes results.

beta max is located . However, on the suction side three distinct spikes from the Lagrangian codes (GlennICE and LEWICE3D) are shown .

These spikes correspond to the bin runs, so the three spikes are due to three different bins or droplet sizes striking those distinct locations further down the chord . Fig. 4 6 shows the results for back calculating ⁄ experimental collection efficiency. In the 𝜌 = 917 𝑘𝑔 𝑚 curve 𝑖𝑐𝑒 Page 13 of 17 UG3540 Sim and Exp Composite Beta vs. Z (in) for 29.9 MVD on Inst. Strut at Y = 6in GlennICE UG3540 Sim and Exp Composite Beta vs. Z (in) for 29.9 MVD on Non-Inst. Strut at Y = -11in LEWICE3D 0.9 1 FENSAP-ICE GlennICE Exp Beta (917 kg/m^3 Ice Density) LEWICE3D 0.8 0.9 Exp Beta (450 kg/m^3 Ice Density) FENSAP-ICE Exp Beta (917 kg/m^3 Density) 0.7 0.8 Exp Beta (450 kg/m^3 Density) 0.6 0.7 0.5 0.6 0.4 0.5 0.3 0.4 0.2 0.3 Collection Efficiency, Beta 0.1 0.2 Collection Efficiency, Beta -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.1 Z (in) Figure 48 . Simulation and experimental beta at the instrumented strut slice -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 Z (in) Y=6in for UG3540 150 knots at 4 ° AoA rime run.

Figure 51 . Simulation and experimental beta at the non - instrumented strut slice Y= - 11in for UG3540 150 knots at 4° AoA rime run.

UG3540 Sim and Exp Composite Beta vs. Z (in) for 29.9 MVD on Non-Inst. Strut at Y = -6in GlennICE LEWICE3D Ice Accretion Analysis at 4 ° AoA and 150 knots 0.9 FENSAP-ICE Exp Beta (917 kg/m^3 Ice Density) 0.8 Exp Beta (450 kg/m^3 Ice Density) Results for a single - shot approach and ice density held constant to 0.7 𝜌 = 917 𝑘𝑔 𝑚 ⁄ results are shown in Fig. 52 . All the data for the 4° 𝑖𝑐𝑒 0.6 AoA ice shapes and clean geometry has been rotated in x, y, and z 0.5 coordinates so th at the clean geometry is displaying at the 0° AoA 0.4 position. GlennICE and FENSAP - ICE did not reach the maximum 0.3 thickness (angled towards the pressure side) of the experimental rime ice shape at the main body “nose” or leading edge . LEWICE3D 0.2 Collection Efficiency, Beta reached the maximum thickness without creating a divot like it did in 0.1 the 0° AoA cases.

-0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 Z (in) Figure 49 . Simulation and experimental beta at the non - instrumented strut 150knts at 4deg AoA Rime Main Body Nose Ice Shape Comparison slice Y= - 6in for UG3540 150 knots at 4 ° A oA rime run.

Fig ures 50 - 51 show the upper outboard strut slice s (Y = 11 in and Y = - 11 in ) ha ve a similar beta max now between the pressure side strut 0.5 and the suction side strut leading edge. The experimental beta curve GlennICE using 𝜌 = 917 𝑘𝑔 𝑚 ⁄ (shown red) matches the icing code better 𝑖𝑐𝑒 LEWICE3D than the 𝜌 = 450 𝑘𝑔 𝑚 ⁄ (shown blue) curve, but the match is FENSAP-ICE 𝑖𝑐𝑒 UG3540 Z=-4in worse than at the lower root slice cuts shown in Fig s . 4 8 - 4 9 . 0 UG3540 Z=4in Y (in) UG3547 Z=-4in UG3547 Z=4in Clean UG3540 Sim and Exp Composite Beta vs. Z (in) for 29.9 MVD on Inst. Strut at Y = 11in -0.5 GlennICE LEWICE3D 0.9 FENSAP-ICE Exp Beta (917 kg/m^3 Density) 0.8 Exp Beta (450 kg/m^3 Density) 0.7 -1 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 X (in) 0.6 0.5 Figure 52 . Computational ice shape vs. experimental UG3540 150 knots at 4 ° AoA run at the main body leading edge for a 5 min icing duration .

0.4 0.3 Figure 5 3 shows the icing code to experiment comparison at the 0.2 instrumented strut “pressure side” and slope junction at the centerline Collection Efficiency, Beta 0.1 (Z = 0). While Fig. 5 4 shows the comparison at the no n - instrumented strut “suction side” and slope junction. More experimental ice 0 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 thickness is on the pressure side slope, Fig. 5 3 , compared to the Z (in) suction side slope, Fig. 5 4 . This is due to the increase in collection Figure 50 . Simulation and experimental beta at the instrumented strut slice efficiency due to the angle of attack cha nge. However, more ice Y=11in for UG3540 150 knots at 4° AoA rime run.

exists on the non - instrumented strut leading edge in Fig. 5 4 due to the high LWC above the shadow zone impacting on the lower non - instrumented strut. The ice bumps on the root of the strut from GlennICE and FENSAP - ICE are artifacts from using limited discrete droplet distribution bins.

Page 14 of 17 LEWICE3D overpredicts the max ice thickness on Fig. 5 6 on the 150knts at 4deg AoA Main Body Junction Ice Shape Comparison suction side.

6.5 GlennICE 6 LEWICE3D 4deg AoA Lower Instr. Strut LE Ice Shape Comparison FENSAP-ICE 0.8 UG3540 Z=0in GlennICE 5.5 UG3547 Z=0in LEWICE3D Clean FENSAP-ICE 0.6 UG3540 Y=6in UG3547 Y=6in Clean 0.4 Y (in) 4.5 0.2 Z (in) 3.5 13 13.5 14 14.5 15 15.5 16 16.5 -0.2 X (in) Figure 53 . Computational ice shapes vs. experimental UG3540 150 knots at 4 ° -0.4 15.5 15.6 15.7 15.8 15.9 16 16.1 16.2 16.3 16.4 16.5 16.6 16.7 AoA run at the main body slope and strut junction on instrumented side for a X (in) 5 min icing duration .

Figure 55 . Computational ice shape vs. experimental UG3540 150 knots at 4 ° AoA run on the instrumented strut LE at Y= 6 in slice for a 5 min icing 150knts at 4deg AoA Main Body Junction Ice Shape Comparison duration .

-3 GlennICE LEWICE3D -3.5 FENSAP-ICE 4deg AoA Lower Non-Instr. Strut LE Ice Shape Comparison UG3540 Z=0in 0.8 UG3547 Z=0in GlennICE -4 Clean LEWICE3D FENSAP-ICE 0.6 UG3540 Y=-6in UG3547 Y=-6in -4.5 Clean 0.4 Y (in) -5 0.2 -5.5 Z (in) -6 -6.5 13 13.5 14 14.5 15 15.5 16 16.5 -0.2 X (in) Figure 54 . Computational ice shapes vs. experimental UG3540 150 knots at 4 ° -0.4 15.5 15.6 15.7 15.8 15.9 16 16.1 16.2 16.3 16.4 16.5 16.6 16.7 AoA run at the main body slope and strut junction on non - instrumented side X (in) for a 5 min icing duration .

Figure 56 . Computational ice shape vs. experimental UG3540 150 knots at 4 ° AoA run on the non - instrumented strut LE at Y= - 6 in slice for a 5 min icing Ice comparisons on the instrumented strut lower root slice (Y = 6 in) duration .

can be seen in Fig. 5 5 , and ice comparisons on the non - instrumented strut lower slice (Y = - 6 in) can be seen in Fig. 5 6 . LEWICE3D got Ice comparisons on the instrumented strut upper outboard slice (Y = closer in prediction than GlennICE and FENSAP - ICE to the 11 in) can be seen in Fig. 5 7 , and ice comparisons on the non - experimental ice shapes on the instrumented strut in Fig. 5 5 .

instrumented strut lower slice (Y = - 11 in) can be seen in Fig. 5 8 .

Here the icing codes match the experimental ice shape better, but none of the three meet the maximum ice thickness. This might be due to using a single shot approach and using a constant ice density.

Page 15 of 17 in good agreement at the different tunnel airspeeds , and between the 150knts Rime Upper Strut LE Ice Shape Comparison 0° and 4° angle of attack cases. The aerothermal thermocouple 0.8 comparisons were closer in matching exp erimental than the 30 sec ond GlennICE LEWICE3D dry air result s before the rime icing cloud hit the test article . This might FENSAP-ICE 0.6 be because the test article ’s skin temperature does not reach full UG3538 Y=11in UG3538 Y=-11in thermal equilibrium between the rime tunnel run s. Future analysis will Clean focus on the thermocouples and the heat flux measurements to 0.4 simulations, and that might narrow down this thermal equilibrium issue.

0.2 Z (in) The ice shape at SIDRM’s leading edge nose is captured using all three icing codes when using a constant ice density of 𝜌 = 917 𝑘𝑔 𝑚 ⁄ , 𝑖𝑐𝑒 but this area could benefit more with switching from a single - shot to -0.2 multi - shot approach in the icing setup. T he experimental beta calculat ion matches the predicted droplet code collection efficiency results at the nose (beta max) when assuming 𝜌 = 917 𝑘𝑔 𝑚 ⁄ in -0.4 𝑖𝑐𝑒 15.5 15.6 15.7 15.8 15.9 16 16.1 16.2 16.3 16.4 16.5 16.6 16.7 Eq. 3.

X (in) Figure 57 . Computational ice shape vs. experimental UG3540 150 knots at 4 ° Changing the constant ice density to 𝜌 = 450 𝑘𝑔 𝑚 ⁄ , improves the AoA run on the instrumented strut LE at Y= 11 in slice for a 5 min icing 𝑖𝑐𝑒 experimental beta calculation match on the main body slope to the duration .

droplet composite (7 - bin) predictions from the three icing codes. To reduce the ice feathers complication, it is suggested to use shorter spray 4deg AoA Upper Non-Instr. Strut LE Ice Shape Comparison time runs to calculate collection efficiency from a rime 3D laser scan 0.8 ice shape. However, for the ice shape comparisons , the computational GlennICE LEWICE3D ⁄ single - shot ice shape overpredicts when changing 𝜌 = 450 𝑘𝑔 𝑚 𝑖𝑐𝑒 FENSAP-ICE 0.6 on GlennICE compared to the experimental laser scan ice thickness.

UG3540 Y=-11in UG3547 Y=-11in This suggests LWC , flow angle, droplet impact angle, or airspeed Clean changes have more of an effect on the main body slope ice . Alternately, 0.4 this may suggest that a single - shot approach is not as accurate for this sloped geometry ice accretion. The main body nose ice shape could 0.2 have a greater effect on droplet trajectories impacting the ramp slo pe.

Z (in) Improvements on gathering experimental ice mass measurements in just this slope area are being implemented.

-0.2 For the strut leading edge, GlennICE and FENSAP - ICE consistently underpredicted the ice shape no matter the chosen density. This strut leading edge area would benefit from a multi - shot approach -0.4 15.5 15.6 15.7 15.8 15.9 16 16.1 16.2 16.3 16.4 16.5 16.6 16.7 comparison at the same constant ice densities. LEWICE3D matched X (in) the experimental lower strut slice (Y = 6 in) ice shape better when Figure 58 . Computational ice shape vs. experimental UG3540 150 knots at 4 ° 3 3 ⁄ ⁄ assuming 𝜌 = 917 𝑘𝑔 𝑚 on the 0° cases , and 𝜌 = 917 𝑘𝑔 𝑚 𝑖𝑐𝑒 𝑖𝑐𝑒 AoA run on the non - instrumented strut LE at Y= - 11 in slice for a 5 min icing on the instrumented strut “pressure side” on the 4° case . LEWICE3D’s duration .

swept wing ice density model gave mixed results on the strut leading edge as well on the 150 knots and 0° case.

Conclusion

The test article angle of attack impacted the location of ice accretion Computational aerothermal dry air and icing simulations were on the main body slope and strut leading edges due to high LWC conducted on the SIDRM test article and compared to the experimental concentrations and airspeed flow changes. FENSAP - ICE and measurements collected from tests performed at the NASA IRT icing GlennICE single - shot results on the stru t slices were identical, but wind tunnel in early 2022. This paper covered an icing study analyzing underpredicted the maximum ice thickness . Changing the ice density different icing codes involving FENSAP - ICE (Eulerian approach), 3 𝜌 = 450 𝑘𝑔 𝑚 ⁄ improves the FENSAP - ICE and GlennICE 𝑖𝑐𝑒 LEWICE 3D (Lagrangian approach), and GlennICE (Lagrangian prediction, but they still underpredicted the maximum thickness. A approach) on SIDRM’s complex body flow - field and results were multi - shot approach might have better matching to the strut ice profile compared to experimental aerothermal results and rime icing tunnel shapes. Utilizing the two constant ice densities bounded the results for run data . Rime icing runs were chosen due to the objective of much of the analysis, and more focus in future tests will be to measure calculating collection efficiency through an empirical equation from the experimen tal ice density on separate locations of this test article.

the laser ice shapes. In addition, single - shot approaches in the icing codes d o better in capturing rime ice shapes compared to glaze due to

References

glaz e icing having runback icing and a varying ice density as the ice 1. Mason, J. G., Strapp, J. W., and Chow, P., “The Ice Particle grows.

Threat to Engines in Flight,” 44th AIAA Aerospace Sciences Meeting and Exhibit , AIAA, Reno, NV, 2006, AIAA - 2006 - 206 Thermocouple and pressure taps results measured on the t est article https://doi.org/10.2514/6.2006 - 206 .

were compared to the computational dry air CFD simulations for 2. Struk, P. M., Bartkus, T. P., Bencic, T. J., King, M. C., surface pressure and temperature validation of the s imulations before Ratvasky, T. P., Van Zante, J. F., and Tsao, J. C., “An initial the icing cloud was turned on. The experimental to CFD results were Page 16 of 17 study of the fundamen tals of ice crystal icing physics in the 20. Tsao, J en - Ching and Lee, Sam. Evaluation of Icing Scaling on NASA propulsion systems laboratory,” 9th AIAA Atmospheric Swept NACA 0012 Airfoil Models NASA/CR - 2012 - 217419, and Space Environments Conference, June, 2017, pp. 1 – 40. 2012 https://ntrs.nasa.gov/citations/20120009186 .

https://doi.org/10.2514/6.2017 - 4242 . 21. Tsao, J en - Ching, and Christopher E. Porter. “Characterization of 3. Struk, P. M., Agui, J., Ratvasky, T., King, M., Bartkus, T., and Collection Efficiency of the Common Research Model Midspan Tsao, J. C., “Ice - Crystal Icing Accretion Studies at the NASA Wing Section in the IRT.” AIAA AVIATION 2021 FORUM, 28 Propulsion Systems Laboratory,” SAE Technical Papers, Vol. Aug. 2021, https://doi.org/10.2514 /6.2021 - 2681 .

2019 - June, No. June 2019, pp. 1 – 12.

https://doi.org/10.4271/2019 - 01 - 1921 .

C ontact Information

4. Bartkus, T. P., Tsao, J. C., and Struk, P. M., “Analysis of Experimental Ice Accretion Data and Assessment of a Eric Stewart Thermodynamic Model during Ice Crystal Icing,” SAE Work phone: (2 10 ) 833 - 4647 Technical Papers, Vol. 2019 - June, No. June 2019.

E - mail: eric.a.stewart@nasa.gov https://doi.org/10.4271/2019 - 01 - 2016 .

Affiliation: Naval Air Warfare Center Aircraft Division (NAWCAD) 5. Bartkus, T. P., Lee, S., Potapczuk, M. G., and Flack. C. A., “D escription of Cloud Characterization and Icing Tests for a 3D Heated Test Article at the NASA Icing Research Tunnel,” AIAA AVIATION 2022 FORUM, AIAA, Chicago, IL, 2022, AIAA - 2022 - 3700 https://doi.org/10.2514/6.2022 - 3700 .

6. Bartkus, T., Potapczuk, M., Lee, S., Stewart, E., and Chen, R. - C., “Plans for Ice Crystal Icing Tests Using a 3D Heated Test Article at the NASA Icing Research Tunnel,” AIAA AVIATION 2021 FORUM, AIAA, Virtual Event, 2021, Oral Presentation.

7. William, Mathias. “General Electric GE90 - 94B Turbofan Engine”. https://grabcad.com/library/general - electric - ge90 - 94b - turbofan - engine - 1. GrabCAD. Nov. 2016.

8. SAE Aerospace. Aircraft Inflight Icing Terminology . SAE International, April 2013. ARP5624.

https://doi.org/10.4721/ARP5624 .

9. “ANSYS FENSAP - ICE User Manual , ” Release 202 1 R1, ANSYS, Inc., 2021.

10. Bidwell, C.S., and Potapczuk, M.G., “ User’s Manual for the NASA Lewis Three - Dimensional Ice Accretion Code (LEWICE3D),” NASA/TM 105974, 1993.

https://ntrs.nasa.gov/citations /19940017117 .

11. “GlennICE User Manual Software Version 2.2.0”, NASA, June 2022 .

12. Bartkus, T. P. and Stewart. E. A., “Icing Physics Studies using the 3D SIDRM Test Article: Aerodynamic and Supercooled Liquid Icing Analysis,” SAE International Conference on Icing of Aircraft, Engines, and Structures , SAE, Vienna, Austria, 2023 (submitted for consideration).

13. “ANSYS Fluent Theory Guide , ” Release 2021 R1, ANSYS, Inc., 2021.

14. Bourgault, Y., Beaugendre, H., and Habashi, W.G., “Development of a Shallow Water Icing Model in FENSAP - ICE”, AIAA Journal of Aircraft , Vol. 37, No. 4, 2000, pp. 640 - 646. https://doi.org/10.2514/2.2646 .

15. Bidwell, Colin. “Icing Analysis of the NASA S3 Icing Research Aircraft Using LEWICE3D Version 2.” SAE Technical Paper Series , 2007, https://doi.org/10.4271/2007 - 01 - 3324 .

16. Wright, W., Porter, C., Galloway, E., and Rigby, D. “An Automated Refinement Process for Particle Trajectory Methods in GlennICE,” AIAA AVIATION 2021 FORUM, AIAA, Virtual Event, 2021, https://doi.org/10.2514/6.2021 - 2631 .

17. Porter, C. E., “A Comparison of Trajectory Refinement Schemes for GlennICE,” AIAA AVIATION 2022 FORUM, AIAA, Chicago, IL, 2022, ht tps://doi.org/10.2514/6.2022 - 3692 .

18. Timko, Emily N, et al. NASA Glenn Icing Research Tunnel: 2019 Cloud Calibration Procedure and Results NASA/TM - 20205009045, 2021 https://ntrs.nasa.gov/citations/20205009045 .

19. Anderson, D. and Tsao, J. Overview of Icing Physics Relevant to Scaling NASA CR 2008 - 213851, 2008.

https://ntrs.nasa.gov/citations/20050215167 .

Page 17 of 17