Kinetic and thermal simulation of water droplets in icing wind tunnels

Within the certification process of aircraft, tests under specific icing conditions are required. For such safety relevant tests—which are performed under defined and repeatable test conditions—specially equipped Icing Wind Tunnels (IWT) are required. In such IWTs, supercooled water droplets are created with the aid of a spray system injecting pre-tempered water droplets of specific diameters into the free stream air flow. Especially tests with a droplet size up to 2mm (Supercooled Large Droplets - SLDs) are of great importance. SLDs are difficult to generate under laboratory conditions in IWT since usually the available droplet flight time from the injection location to the impact position on the test object is insufficient to reliably cool down a droplet at least to freezing temperature. To investigate the limitations associated with the application of SLD, the current work provides a method to allow detailed insight into the behavior of droplets on the path from the injection spray nozzle to the test section. In this work a state space model of a single droplet is derived that combines the kinetic aspects, thermal properties as well as the governing differential equations for motion, convective heat transfer at the droplet surface and heat conduction inside the droplet. Beside the states for the droplet’s position and velocity in space, the state space vector comprises various fluid and thermodynamic parameters. The droplet-internal temperature distribution is modelled by a discrete one-dimensional spherical shell model that also incorporates the aggregate phase (freezing mass fraction) at each shell node. This approach allows, therefore, the simulation of potential droplet phase change processes (freezing/melting) as well. With the model at hand, the influence of various boundary conditions (initial droplet temperature, flow field, ambient air temperature, etc.) can be determined and evaluated. As a result, concrete measures to achieve a desired operating condition (e.g. droplet temperature at the test object) for various model assumptions can be derived. In addition, the simulation model facilitates the prediction of the droplet diameter threshold for ensuring a supercooled state upon the impact on the test object. The governing theoretical influences are described, and various simulation results for representative test conditions that occur at the Rail-Tec-Arsenal (RTA) in Vienna are presented.


Introduction
The problem arising due to icing of aircraft components is inevitable, and various fields of research and experiments are dealing with an exact description of the complex processes involved. There exist several approaches that try to tackle the process of ice accretion in a numerical way [1].
For the application of SLD testing within Icing Wind tunnels, two major concerns arise: Firstly, the overall heat transfer to reach a supercooled state of the droplets must be achieved, and secondly, gravitational influences and momentum of the droplets may restrict the usable test section area [2]. To allow for the investigation of the temperature of a water droplet within an IWT such as the RTA icing wind tunnel and for the general investigation in the layout of IWT, the described simulation model was set up. With the aid of this model, the actual temperature of a water drop can be predicted. It can be checked, whether it is theoretically possible that a water droplet reaches a "supercooled" state at the point where it impacts on a test specimen. To achieve a simulated time history of the droplet position in space, the temperature distribution within the droplet shell model and the freezing mass fraction, the governing differential equations are numerically solved. The flow field itself may be of arbitrary origin such as CFD calculations, analytical calculations or measurements. For the presentation of exemplary results within this work, the flow field is determined by means of CFD calculations. As the underlying phenomena of the cool down history as the main contribution of this work can be transferred to arbitrary layout of future IWT to be designed, any quantitative assessments of the flow field properties itself are not within the scope of this work.
A common model is documented within the work of Willbanks and Schulz [3], describing an one-dimensional, multiphase flow code (the "AEDC-1DMP" code) also incorporating a partial freezing model that is used in various studies by Al-Khalil et al. [2] and Orchard et al. [4]. The AEDC-1DMP also includes model assumptions reducing its applicability to the simulation of temperature and trajectories of SLD. The work of Orchard et al. includes an overview of simulation methods documented in literature. Various approaches incorporate evaporation and partial freezing models while having the scope of e.g. predicting the LWC content at a given test section [5] or investigating on droplets trajectories and the effects of the flow field on the cloud density variation within axisymmetric wind tunnels [6] or for predicting ice and water deposition in aeronautical mixing manifolds [7]. While basically following the same approach for the kinetic simulation of droplets exposed to a flow field, the current work is not limited to an analytic expression for the flow field. Extending the applicability to arbitrary flow fields from various sources greatly enhances the field of application to arbitrary wind tunnel layouts. Most importantly, the aforementioned models assume a lumped thermal drop model and to the best knowledge of the authors, no comparable study of inner heat transfer of large droplets within an IWT as described in the current work is available.
The approach chosen here differs from the previously methods and the current work extends the scope of the simulation and the model assumptions to better fulfill the specific requirements of SLD simulation for IWT design tasks: • The trajectory, relative velocity and subsequently heat transfer calculations are based on a separately performed CFD calculation and are capable of handling 3-dimensional flows. Separately calculating the flow field allows more elaborated calculations including complex wind tunnel geometry, varying cross sections, fixtures etc. as well as sophisticated turbulence models. • Gravity is of importance for the trajectory and impingement angle of SLD and is thus incorporated in the model at hand. Studies of trajectories of droplets within a 3-dimensional flow field emphasize the importance of considering gravity, especially for SLD [4,6]. • Internal heat transfer has to be considered for larger droplets having a Biot number Bi greater than 0.1. • The injection position, velocity and angle can be varied arbitrarily to simulate a nozzle spray cone. • The simulation capabilities are extended to arbitrary flow fields.
Another motivation for a deeper investigation of the proposed shell model is to estimate boundaries for droplet temperatures within an IWT: The lumped model approach assumes a perfect heat transfer within the drop, resulting in a higher interface temperature and subsequently higher heat loss of the drop. Thus the cooldown potential of a drop traveling through the wind tunnel tends to be overestimated. While apparently underestimating the cooldown potential with the proposed shell model, incorporating heat conduction but no convection, an upper bound of the droplet temperature is obtained.

Methods
Within this section the determining parameters of the model for a single droplet are described. The model consolidates the kinetics and thermal behavior to a single state space model.

State space model of the droplet
For the simulation of the time history of a droplet, various properties are consolidated to a state space model [8]. The state vector consists of the following parameters: With x,y and z the coordinates of the actual position are denoted and u,v and w denote the corresponding velocity components. The state space contains 3-dimensional position and velocity components to maintain data integrity for a future use of a 3-dimensional vector field that is not implemented yet. The droplet diameter and the mass are saved in the parameters d drop and m drop . They may change with evaporation and water density change during simulation. The relevant fluid and thermodynamical parameters are the Reynolds number (Re), the Nusselt number (Nu), the Biot number (Bi), the Weber number (We), the Ohnesorge number (Oh) and the Prandtl number (Pr) as well as the local convective heat transfer coefficient ( h c ) . These states are required in the computation of the the governing differential equations during computation and are stored in the time series of the results. Since they mostly depend on local relative velocities, a recalculation using the flow field would be required otherwise.

Position
The change of position results from a pure integration of the actual velocities. This yields: No rotational motion of the droplet is considered and thus no angular position is required.

Velocity
The velocity components of the droplet are given by: The drag force is given by represents the magnitude of the drag force and − rel | rel | determines the direction of the drag force opposing the actual relative velocity vector rel . The slip velocity vector rel is calculated as the difference of the actual velocity and the local velocity of the airflow field that is known.
In addition, a constant gravitational force is added in vertical direction. The change of velocity ̇ (i.e. the acceleration) is given by the attacking force divided by the mass which yields: For coefficient of drag ( C d ) calculations of smooth spheres, various formulas can be found in literature [1,[9][10][11]. While for low Reynolds number flow ( Re < 1 ), the correlation C d ∼ 24 Re (Stokes Law) is common among the sources, the values vary for higher Reynolds numbers. Given the stated uncertainty of single-expression curve-fit errors of up to 5% [9], up to 10% [10] or discontinuities [1], a mixture of the provided formulas seems a reasonable approach.
Munson et al. offer a graph that merges various sources over the relevant range of Reynolds numbers [12]. Finally, C d is calculated as a function of the local Reynolds number by interpolating a lookup table based on the curve depicted in Fig. 1.
To compensate for the small diameter of the droplet relative to the mean distance between air molecules, the Cunningham correction factor f cun is taken into consideration [13]. The Cunningham correction factor is given by: Fig. 1 The drag coefficient as a function on the Reynolds number (modified from [12]) For small droplets of 1 μ m the influence is roughly 15%, for droplet size over 20μ m the effect vanishes.
Therefore the drag coefficient considering the slip phenomenon is altered to: For sake of completeness, the buoyancy force is given: The buoyancy force of a water droplet within air is negligible.

Temperature
The temperature of the droplet changes due to convective heat transfer to the surrounding air flow, and the Reynolds number (based on the relative velocity of the droplet to the surrounding air) is used to calculate the Nusselt number. The Reynolds number is given by: The Nusselt number is calculated according to the formula of Brauer/Sucker for a sphere [14,15] valid for a range of 0 ≤ Re ≤ 3 ⋅ 10 5 : With the calculated Nusselt number, the heat transfer coefficient h c is calculated by The Prandtl number is calculated using and is considered as constant for the whole simulation, since the temperature of the surrounding flow field is considered The thermal conductivity air is calculated using the following relation [16, p.524] and is kept constant ( T air = const ) over the whole simulation.
The actual temperature change for the water droplet in the lumped model is then given by: where the heat capacity of water is given as a function of temperature according to Fig. 2.
A change in temperature results in a change of the density, which in turn yields a change of the diameter. While the diameter change due to the change of density may be negligible, the mass of the drop over the simulation is included in the state space to enable consideration of evaporation in the extended model.
The density of water is implemented as an interpolation function of the graph depicted in Fig. 3.
Additional numbers are included into the state vector to enable further examination of the results even if they are not explicitly required for solving the differential equations: The Biot number Bi to emphasize applicability of the shell model versus the lumped model approach, the Weber number We to observe limitations for droplet break up [18] in further research focusing on the layout of a wind tunnel specifically designed for SLD [19].
(13) air = 7.655 ⋅ 10 −5 ⋅ T air + 0.024139 The Biot number is defined in the current context by [20]: The Weber number is defined by and the Ohnesorge number-the ratio of viscous forces to surface tension and inertia forces-is given by where water is the dynamic viscosity of water at the mixed mean temperature T drop . The Weber number and the Ohnesorge number are implemented for future use and additional analysis of the droplet.

Initial values
The initial values for the numerical simulation are set according to Table 1. Fig. 3 The density of water as a function of the temperature for various values of the ambient pressure P [17]. The figure presents experimental results (rectangles and circles) and values from the IAPWS-95 formulation (continuous lines) . The line at P = 0.1 MPa is used as a basis for interpolation

Flow field
The governing flow field and the surrounding temperature is taken from a precomputed CFD solution given at discrete data points in conjunction with a spatial interpolation. Based on a 3D CAD-model, a mesh was created with the commercial software package ANSYS ICEM CFD.
For simple geometries a blockstructured hexahedral mesh and for more complex geometries an unstructured tetrahedral mesh with prism layers in wall vicinity was used. For both approaches the thickness of the first element adjacent to the wall was adjusted to yield (depending on the flow velocities of the operating condition) a Y+ value between 0.1 and 10. With the help of mesh smoothing steps aspect ratios <100, face angles >20 • and determinants >0.3 were aimed to ensure high quality meshes.
Then, the convective flow field was determined by a 2D or 3D steady state CFD simulation, performed with the commercial software package ANSYS CFX.
For modelling the turbulent flow structure, the Reynoldsaveraged Navier-Stokes (RANS) approach using the Shear Stress Transport (SST) eddy-viscosity turbulence model from Menter [22] was used. The model works by solving a turbulence/frequency-based model ( k − ) at the wall and k − in the bulk flow. Menter developed advanced nearwall treatments and blending functions, that ensure a smooth transition between the two models within the ANSYS CFX software. The SST model performance has been studied in a large number of cases. In a NASA Technical Memorandum [23], SST was rated the most accurate model for aerodynamic applications.
In the case of design studies, an inlet mass flux was set as inlet boundary condition (BC). All walls are treated as hydraulically smooth with zero slip conditions. For all outlets zero relative static pressure BCs were used.
Initial conditions were adopted for various operating cases to ensure a robust and fast convergence.
All CFD simulations were performed with a continuous fluid air (treated as ideal gas) including buoyancy forces with the total energy heat model (including viscous work) to take into account compressibility effects.
The advection scheme and turbulence numerics are set to a high resolution discretization scheme. All simulations were performed with double precision and the target of a maximum Root Mean square (RMS) for the residuals for mass, momentum, energy and turbulence of 10 −4 . In addition, various monitor points inside the fluid domain reported a stable behaviour of the velocity components to approve the results as a steady state converged solution.
Calculated flow fields are available for measurement section velocities having the following values (in m/s): and for measurement section static temperatures An example flow field that is used for calculations within this work is depicted in Fig. 4. For more detailed information on the flow field's properties, interested readers are referred to to the work of Breitfuss [24].
The initial value of the z-coordinate of the droplet corresponds to the elevation of the nozzles on the spray bar grid and ranges from -1.43m to 2.41m with a vertical bar distance of 0.11m.

Droplet size
The occurring droplet size spectra for the relevant certification specification scenarios given in CS-25 [25] were measured at the RTA Icing Wind Tunnel using a Malvern Instruments Spraytec system. The droplet spectra for the following scenarios were available for the current simulations [  The measured droplet spectra are depicted in Fig. 5.
Since the droplet simulation is based on a single-droplet formulation, the actual mass fraction distribution of a specific droplet size is of no influence for now. Nevertheless, it may be used for assessing the relevance of the results obtained for a certain droplet size.

Lumped thermal droplet model
The lumped thermal droplet model assumes an uniform temperature distribution within the droplet as a result of instantaneous heat conduction. Such an assumption is only valid for droplets with a Biot number satisfying Bi < 0.1 . At greater Bi numbers, the lumped thermal model is not valid anymore [20].
It can be shown (cf. (15)) that for a relevant initial slip velocity of up to 40 m/s a Biot number greater than 0.1 occurs for droplets greater than d 0 ≈ 500 μm . For SLD, therefore, an extended thermal model has to be introduced.

Extended thermal model of the droplet
As stated, the internal heat conduction of the droplet has to be taken into account for larger droplets within the relevant droplet size spectrum.
While the kinetic equations and the heat transfer on the droplet surface remain unchanged also for the extended thermal model of the droplet, the internal temperature distribution of the droplet is modeled as a spherical shell model. The derivation of the shell model is explained in the following sections.

Conservation of energy
The energy balance equation for a control volume V of one single ice water shell is given by the expression Here e is the specific enthalpy within V, q cond is the conductive heat flux and iw is the ice water shell density defined by: Having the thermal conductivity of the ice water shell k iw and the temperature gradient ∇T between adjacent volumes (i.e. neighbour shells) and considering q cond = −k iw ∇T , one subsequently obtains: or In analogy to 21, the thermal conductivity k iw of the ice water shell is given by: The grid size is assumed to be uniform over d 2 with a radial cell size of Δr = d

2⋅N iw
, where d is the droplet's diameter and N iw is the node of the most outward spherical shell as depicted in Fig. 6. The reference state of the water species having e = 0 is considered to be liquid water at T melt .
The following derivation of the formulas vary for three cases: • A regular inner shell with two adjacent cells. • The innermost shell with the node containing the droplet center • The outermost shell with one regular adjacent cell to the inside and the air/droplet interface at the outside.
Consolidating the formulas given above, the discretized transient energy balance equation for an inner spherical shell can be derived: The full explicit transient discretization of (27) for an inner spherical shell at the timestep k is given by: The parameter is defined by: In (28) k k iw,j+1∕2 and k k iw,j−1∕2 are evaluated by (27) iw e k+1 j − e k j = Δt r 2 j Δr The energy balance for the innermost spherical shell is given by: Since the innermost spherical shell exhibits a symmetry around the droplet center, the heat flux at the center is set to q k 1∕2 = 0 (Adiabatic boundary condition). Thus, (31) can be reduced to: The energy balance for the outermost spherical shell (denoted by the shell node index N iw ) is then given by: Again, k k iw,j−1∕2 can be evaluated according to (30b), and is given by (29).
The parameter q k N iw +1∕2 is the convective/evaporative wall heat flux at the air/droplet interface which can be broken up into two parts: The convective heat flux at the air/droplet interface is denoted by q conv and the evaporative heat flux by q lat .
The wall heat flux directly influences the air droplet interface temperature and thus a Dirichlet boundary condition is introduced: The convective wall heat flux q conv requires the heat transfer coefficient h c given in (11) and the temperature gradient at the surface as the difference between the ambient air temperature T a and the surface Temperature T N iw +1∕2 and is calculated by: The evaporative heat flux q lat is given by: with a,humid as the density of humid air and L e,s as the latent specific heat of evaporation and sublimation. By use of a limiting factor that depends on the relative humidity of the surrounding air and by implementation of a as well as saturation of evaporation is also considered. A thorough description of the evaporation process is available in [20] and a description of the implementation is available in [28]. A detailed recap here would exceed the scope of the current work.

Conservation of mass
The change of the droplet's mass with time is given by the expression with the right hand side similar to the wall heat flux as stated in (37).
The corresponding mass change of the surrounding air of the droplet is consequently given by:

Stability criterion
The implemented extended thermal model is based on an explicit finite-difference representation and thus vulnerable to numeric stability issues while solving the implemented differential equations. Therefore a specific stability criterion has to be fulfilled as stated in [27].
with Δx as the grid spacing (or relevant length i.e. Δx ∼ d droplet 2N shells ) and the thermal diffusivity p of the particle given by: Considering a realistic scenario with an extreme value of a droplet diameter of d = 1μm this results in (37) q lat = h m a,humid L e,s a − The maximum calculated time step has to be considered in the implementation of the solver or the parametrization of adaptive time step routines.

Numerical solver
For numerically solving the given differential equations, the MATLAB-implemented solver "ode23t" is used since it can be set to handle "state variables" that shall not be integrated during the solution procedure but are rather just updated. As a basis for the set of differential equations various numbers of interest such as Weber number, Biot number, Nusselt number etc. are calculated, and no additional computational effort is required after calling the numerical solver routine. For the time step adaption, the absolute tolerance has been set to values of d drop ∕100 for all components of the droplets position, 10 −8 for all components of the droplets velocity, 10 −4 for the temperature and 10 −6 for the mass. An overall relative tolerance of 10 −4 has been set for the solver. The accuracy is checked during numerical integration by comparing the achieved results with the results of a more accurate solver of higher order at specified time intervals. Whichever of the set tolerances is to be violated triggers a time step refinement of the solver. Conversely, if all criteria are met, then a relaxation of the time step size is performed [29].

Results
The simulations where conducted for the relevant scenarios as stated before (Sect. 2.4). For a further detailed discussion the results for a CS-25, Appendix C simulation and a CS-25, Appendix O (SLD) simulation are given below. The used parameters for the simulation are given in Table 2.

CS-25 Appendix C simulation results
For the scenario " MVD = 20 μm " droplets with a size ranging from 7.5 μ m to 135μ m are used. This range includes 99% of the measured droplets (see Fig. 5a. The majority of the droplets reaches a temperature below 0 • C within a flight distance of around 4m after the injection point (Fig. 7a). Generally, the velocity of the droplets follows the free stream velocity, as shown in Fig. 7b. The effect of evaporation and density change on the droplet diameter amounts to a rather negligible reduction of 2.5% compared to the initial diameter as depicted in Fig. 7c. The equidistant distribution of the droplets over the injection array is mainly preserved, that is, no major deviation from cloud uniformity is predicted by the present simulations (Fig. 7d). The increased thermal inertia of the droplets from the SLD spectrum (see Fig. 5d) causes a shallow cool down behavior and therefore requires pre-cooled water at considerable low temperatures. By pre-cooling the water to 2 • C and at a static air temperature of -13 • C within the wind tunnel, the mean droplet temperature of all droplets of the considered diameter spectrum cool down at least to freezing temperature before reaching the red marked test section (see Fig. 8a). However, the larger droplets do not acclimate to the ambient air temperature to such an extent as observed for a CS-25, Appendix C scenario which can clearly be seen by observing the position where the droplet center reaches 0 • C in Fig. 8a. The velocity of the droplets does not completely reach the flow field's velocity (Fig. 8b).

CS-25 Appendix O simulation results
The effect of evaporation and density change on the droplet diameter amounts to negligible reduction of less than 0.02% compared to the initial diameter as depicted in Fig. 8c.
Two effects are observable at the droplets' trajectories for CS25, Appendix O simulations. Firstly, the larger droplets do not completely follow the flow field at the contraction vane due to higher inertia and thus impinge on the icing nozzle (red line) at the position of the airflow contraction ( Fig. 8d at x ≈ 4.5m, y ≈ 2m ). Secondly, the droplets trajectories tend to divert downwards close to the test section.   Although this decline may be partly explainable by the slightly downward oriented velocity vectors at the test section, gravitational effects influence the trajectories of larger droplets too. An detailed view of the flow field in the vicinity of the test section is given in Fig.9. The deviation of the flow from horizontal is compensated by adding an offset to the tilt of a test specimen to achieve the desired angle of attack, i.e. conducting a zero angle of attack test not automatically implies a horizontal alignment of the profile itself. Such phenomena are also reported within the work of Breitfuss et al., where the effectively usable area with an acceptable cloud uniformity is stated to be approximately 0.35 m in height and 2 m in width [24]. Precise knowledge of this region is required for placement of test specimen but the area is sufficiently large for wing profile icing tests. However, the droplet trajectories (Fig. 8d) confirm RTA's capability [26, p. 8] to transport the larger droplets from the injection to certain regions of the test section.
The droplets' mean temperatures at the test section entry are consistent over the origin of the droplets from various spraybars as depicted in Fig. 10a. Regarding the velocity of larger droplets, a dependence on the origin of the droplet can be observed, whereas the smaller droplets practically follow the free stream flow field (see Fig. 10b). While these results are given only for a CS25-O simulation, the results for relative velocity are also valid for other scenarios.

Comparison of the lumped thermal model and the extended thermal model
As expected, the differences between the two models considered is relevant for larger droplets of the relevant spectra. In Fig. 11a the Biot numbers along the distance travelled for various sized droplets originating from the nozzle at z = 0.106 m are depicted. Droplet diameters of up to 500 μ m show Biot numbers below the horizontally red marked threshold of 0.1 (see Sect. 2.6) and confirm, therefore, the validity of the lumped approach. As expected, both models indicate matching simulation results in this droplet diameter domain. For diameters larger than 500 μ m, however, the critical Biot number of 0.1 is exceeded and the dropletinternal heat conduction included in the extended approach results in an increased cool down distance compared to the lumped approach. The difference in droplet temperature at the test section is approximately 0.2 • C, the position where the value for the center shell reaches 0 • C is shifted by ≈ 1m for 500 μ m, by ≈ 5 m for 1000 μ m and even more for larger droplets (see Fig. 12).
For the more critical case "CS25, Appendix O" with simulated droplets sizes up to 1143μ m (see Fig. 11a) the observed Weber numbers are below 7. Critical locations for   Fig. 11a, at position x ≈ 7m ). A synopsis of various studies consolidates no breakup for a Weber number below 13 [30]. Even though these studies eventually address different ranges of density ratio = gas liquid , it is reported that the fundamental behaviour and underlying mechanisms of droplet break-up are valid for a wider range of [31] and results are applicable for the problem at hand. Based on the work of King et. al. [18], the Weber number should not exceed 15 to prevent droplet breakup.
A perfect sphere droplet shape is assumed within the model of the kinetic simulation and its applicability may be assessed as follows. An oscillatory deformation may be occurring temporarily in a phase of the droplet's trajectory approaching the maximum Weber numbers. For a Weber Number of 2.1 a deformation (i.e. the maximum ratio of the droplets cross stream diameter to its initial diameter) of 1.2 (20%) is reported [32]. For a deformation of 20%, the drag coefficient increases from 0.44 to approximately 0.6 for 1000 < Re < 2500 , meaning an increase by approximately 35% [33]. The Reynolds number as a function of the distance traveled from injection for a CS25, Appendix O simulation is depicted in Fig. 11b. The simulations within this work consider the transport of a droplet from the injection nozzle to the specimen under investigation. It is known, that in the vicinity of an airfoil major droplet deformation and even breakup may occur due to high accelerations [34]. These processes take place near experimental targets and need not be considered in simulations used for the overall design and layout of an IWT.
The observed temperature differences between the lumped and the extended shell model prompts the conclusion that the influence of a droplet's inner heat transfer on the resulting temperature is not the main contributing factor to droplet cooldown.
From a thermodynamic view, the potentially occurring oscillatory deformation [30] and the associated internal flow leads to an increased internal droplet heat exchange (forced convection) which results in approaching the lumped model assumption. An increased droplet surface may lead to enhanced heat transfer, too. The conclusion that can be drawn from this is that the extended model is likely to underpredict the cooling rate of droplets.
It is reported that the droplet's actual temperature has little impact on the resulting shape of the accumulated ice [2] which is in close agreement with preliminary results of ongoing research of the authors. Simulations with the ICEAC2Dv2 Icing Code [35] have not proven major shape variations for droplet temperature differences within the observed variability as in Fig. 12.
This can be attributed to the fact that convective cooling of the airfoil and the just accumulated ice layer is a major contributor to the overall icing process. These assumptions of course are limited to an extent, up to when droplets are (partially) frozen. For the analysis of ice accretion on a surface under investigation, the freezing fraction has significant impact on the resulting ice shape due to completely different processes involved at impingement [36]. However, partial or complete freezing of water droplets is a frequently observed phenomenon in icing wind tunnel experiments. This the primary reason why phase transitions were considered in the current model and the ice mass fraction of the droplet approaching a surface are calculated for future utilization. In certain types of tests it is even desired to produce frozen particles rather than supercooled droplets (e.g. CS-25, Appendix P, "Mixed phase and ice crystal icing").

Conclusion
Experimental icing tests are essential for the certification of aircraft based on the icing-related Appendices C and O of the CS-25 published by the European Aviation Safety Agency (EASA) [25]. The regulations shall prompt aircraft manufacturers to implement reliable and effective protective measures against hazardous icing conditions.
For demonstrating compliance to both CS-25 Appendices, manufacturers perform expensive flight experiments as well as experimental icing tests at appropriate facilities. Especially icing caused by Supercooled Large Droplets (SLD) with diameters of up to 2mm (CS-25-O) requires a profound knowledge of the kinetics and the thermal aspects a droplet experiences between its injection into an Icing Wind Tunnel (IWT) and the impact on the affected structure (e.g. a wing). SLDs are critical for experimental icing tests at IWTs since the cool down process to freezing temperature takes place relatively slowly. This paper presents a simplified (lumped) and an extended thermal model aiming at the prediction of a droplet's temperature upon reaching the test section within an IWT. The lumped approach assumes infinite droplet-internal heat conduction processes and is thus limited to the simulation of droplets smaller than approximately 500 μm.
The extended thermal approach, however, accounts for the heat conduction inside the droplet and allows the thermal simulation of SLDs of the CS-25-O spectrum as well. In addition, mass transfer aspects between air and droplet as well as potential phase change processes (freezing/melting) are implemented. The extended thermal model inevitably requires additional computational effort. It has been found that execution time is increased by roughly a factor of 20. Even though a single trajectory is simulated within several seconds, the required simulation time can play a significant role when using the model for a multi-dimensional optimization problem, e.g., the overall design of an IWT layout in the planning phase. The comparison of both thermal droplet models indicated matching results for small diameters but also justifies the higher computational effort using the extended approach for larger droplets.
The results have shown differences for the droplet's temperature using both of the models at hand. These differences in calculated temperature subsequently only have little impact on the simulated ice shapes. Even if resulting ice shapes are similar from a geometric point of view, hence having no major impact on aerodynamic properties of the device, the mass averaged enthalpy of a droplet cloud influences the course of subsequent phenomena. For instance, the energy level directly influences the required heat flux and the effectiveness of deicing systems and their overall layout.

Outlook
For the future, additional validation efforts are contemplated. The Department of Aerospace Engineering at the Pennsylvania State University is currently investigating the freezing and melting behavior of single droplets by means of luminescent imaging [37] which might result in meaningful validation data for the phase change routine of the extended thermal droplet model [38]. It has to be further investigated, whether the results obtained may carry over to the application of the current model regarding the different range of Reynolds number for the cases considered.
It was shown that the model at hand provides a unified investigation approach of the influence of the flow field on the flow path as well as the thermodynamic properties of the droplets. In conjunction with investigations on the real flow field [24], further improvements on the overall capabilities may be achieved while possibly requiring only minor modifications (e.g. vanes, elongation of floor coverage at the orifice).
Moreover, the already investigated validation approaches included tests with even lower temperature within the IWT. If the critical temperature [11] of a droplet is reached, no ice accretion should be observed at a test specimen at a given distance since the then solid drops bounce off the surface. However, such tests are a precarious method and have not given fully satisfactory results so far and thus may need to be refined in a future research.
The model developed in the current research is used for the overall design and layout of a planned IWT specialized on SLD capabilities where the model has shown its potential to efficiently simulate diverse flow field properties within a multidimensional optimization [19]. The relatively low computational effort (when compared to a full 3D finite element method of a droplet) while maintaining an enhanced thermal model of a single drop also opens up possibilities to perform trajectory and thermal simulation for a whole regime of droplets. Once the time consuming calculation of a transient flow field solution is performed, even such a flow field may serve as a future basis for droplet trajectory simulation. Such investigations require only minor adaptions to the underlying data structure. By incorporating an assessment of the droplets impingement points and binning, assumptions on cloud uniformity of planned IWT facilities or measures for further improvement of existing facilities may be derived in future work.
Funding Open access funding provided by FH Joanneum -University of Applied Sciences.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.