On the influence of ice accretion over the aerodynamic performance of airfoils

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Máster Máster Universitario en Ingenería Aeronáutica On the influence of ice accretion over the aerodynamic performance of airfoils Autor: Antonio José Díaz García Tutor: Dr. Francisco José Blanco Rodríguez Tutor: Dr. Francisco Cruz Mazo Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 Trabajo Fin de Máster Máster Universitario en Ingenería Aeronáutica On the influence of ice accretion over the aerodynamic performance of airfoils Autor: Antonio José Díaz García Tutores: Dr. Francisco José Blanco Rodríguez Profesor Titular de Universidad, Universidad de Sevilla Dr. Francisco Cruz Mazo Profesor Ayudante Doctor, Universidad Politécnica de Madrid Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 Trabajo Fin de Máster: On the influence of ice accretion over the aerodynamic performance of airfoils Autor: Antonio José Díaz García Tutores: Dr. Francisco José Blanco Rodríguez Dr. Francisco Cruz Mazo El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos A mi familia, que gracias a ellos no sería lo que soy hoy, y me han proporcionado desde el día uno un apoyo y una ayuda que no se puede medir por palabras, sacrificándose por mí en todo momento y haciendo que nunca me faltara de nada. A mi gente que durante esta andadura universitaria me han ayudado a seguir adelante, especialmente a los que llevan conmigo desde el principio, tanto en los momentos buenos como en los malos, llevándome a ser la persona que soy hoy en día. Y por último, a mis tutores Fran y Paco, que desde el día uno con su pasión y atención, han conseguido que haya disfrutado y aprendido mucho durante la realización del trabajo y consecuentemente éste haya llegado a buen puerto. Antonio José Díaz García Sevilla, 2023 I Resumen La formación de hielo sobre superficies, concretamente, sobre las superficies aerodinámicas, es un problema al que se ha enfrentado la aviación desde prácticamente sus inicios. Este proceso es complejo y difícil de estimar, puesto que depende de multitud de factores como las características del vuelo, la geometría del perfil o el ala, las condiciones atmosféricas, etc. El desarrollo en los últimos años de la Mecánica de Fluidos Computacional (CFD por sus siglas en inglés) ha permitido abordar de una forma más precisa este fenómeno. En este documento se tiene como principal objetivo la realización de un análisis exhaustivo de la formación de hielo sobre superficies, tanto aerodinámicas como los perfiles, como de superficies (canónicas) más simples, como la placa plana o el cilindro, con el objetivo de caracterizar el comportamiento de la formación de hielo sobre estas superficies en base a la variación de parámetros cruciales como la velocidad del flujo, temperatura, geometría del perfil o de la placa plana, etc. Para ello, se ha empleado Ansys FENSAP ICE . Este software es un programa que nos permite simular la formación de hielo sobre distintas superficies, bastando con introducir la geometría de nuestro problema, malla y las condiciones de contorno, tanto de aire incidente como del agua (que definen las gotas). Este documento, además, pretende servir como base para futuros estudios relacionados con la formación de hielo, y es por ello que se ha realizado una breve explicación general sobre los distintos módulos que conforman Ansys FENSAP ICE , así como de sus principales características y funcionalidades. Con estos dos objetivos en mente, se quiere enfatizar la importancia que tiene en la aviación la formación de hielo sobre las superficies aerodinámicas, puesto que pueden afectar negativamente a la actuación (sustentación y arrastre) de la aeronave, y cómo afecta a la formación de hielo los parámetros principales de los que dependen este fenómeno. III XChapter 0. Notation Subscript a Air d Droplet f Fluid i Ice S Solid w Wall ∞Reference 1 INTRODUCTION In this section, a revision of the ice accretion study history will be made, starting with the beginning of the 30s with the experimental studies and the initial stages of the Computational Fluid Dynamics (CFD) until recent years and state of art of this phenomena. Then, the objectives of the study will be stated, as well as the layout of this document, serving as a brief summary of the chapters presented in this project. 1.1 Ice Accretion Background Ice accretion problem on aerodynamic surfaces is a problem that has been known since the early years of aviation. This paragraph will summarize ice accretion research over the years, from the early studies made since 1930 until recent and current research in this field. Bragg et al. ([ 2 ]) divide ice accretion research into three main stages since two main events change entirely the focus of ice accretion research, which are the NASA workshop (1978) and the ATR-72 accident due to ice accretion (1994). Those stages are the following: •Early years: Performance effect and NASA Workshop (until 1978) •Progress in Ice Accretion Research: ATR-72 accident (until 1994) •Recent years and actual years (from 1995 to up to date) 1.1.1 Early studies: NASA workshop Ice accretion problem research started in the 30s, with experimental studies on different airfoils. During those years, ice accretion was seen as an operational problem and one of the most important ones. Thus, the main focus of the research was the study of the effect of ice accretion on the aerodynamics performance of airfoils, significantly lift and drag. This led to almost exclusively experimental research during these years since other critical future techniques, such as CFD, still needed to be fully developed. The construction of wind tunnels to study this phenomenon became critical to take data from aerodynamic degradation due to ice accretion. Therefore, during those years (the 50s-60s), the first data about performance degradation was obtained, even under different conditions. Even though fundamental data was collected from icing tests during those years, those tests and data were collected from a unique type of airfoil. Consequently, since the ice accretion phenomenon is complex, it could have been more useful for other airfoil types. However, in the 70s, aviation started to grow a lot. Thus, more problems related to ice accretion occurred since a wide variety of aircraft were operating (different performance, wings, etc). It was 1 2Chapter 1. INTRODUCTION especially noticeable in rotorcraft and general aviation. This led to an increasing interest in ice accretion research, culminating in a workshop sponsored by the National Aeronautics and Space Administration (NASA) and Federal Aviation Administration (FAA) in 1978 . This workshop was a meeting in which the most influential personalities of this field engaged in an extensive discussion of ice accretion phenomena. Therefore, a new horizon and focus on ice accretion research was made, especially since new technologies such as Computational Fluid Dynamics (CFD) became available and developed. 1.1.2 Progress in Ice Accretion research: ATR-72 accident From 1979, ice accretion research started to grow interested in other regions, such as Canada or Europe, thanks to the workshop held in 1978 in the United States. As said previously, CFD methods became fully developed, and thus a new field of studies became available. NASA (and other researchers) focused on aviation-related areas, such as ground facilities, flight tests, ice protection devices, and icing aerodynamics. Consequently, the new focus was developing new CFD methods and validating them using experimental measurements of airfoils and wings in icing tunnel tests. It was mainly focused on performance degradation and the flow field around the airfoil and wing During those years, numerous ice tests were made to get the most significant amount of data possible to validate and compare it with the numerical data and solution from the new CFD methods. It was mainly done in the NACA 0012 airfoil, and it was studied under different conditions of LWC, droplet diameter, temperature, velocity, etc., to have enough data. However, the research was done almost exclusively on the symmetric NACA 0012 airfoil. Figure 1.1 ATR-72 Aircraft. In 1994, the ATR-72 accident led to a new focus on the ice accretion problem. This accident was due to the severe icing conditions that the aircraft experimented with during its flight, leading to losing control and crashing into the ground. Concretely, investigation was concluded that slight unevenness in the ice accretion on the wing upper surfaces aft of the leading edge and ahead of ailerons towards each wing tip had, when the flaps were retracted and angle of attack had consequently increased, caused differential airflow separation and disturbed the normally balanced aileron hinge moments so as to cause sudden self-deflection. Hence, this accident was due to Super Large Droplet (SLD) conditions, which were new at that moment. Consequently, focus on studying different ice accretion shapes and geometries and critical ice accretion conditions such as SLD 1.1 Ice Accretion Background 3 conditions started to develop. Moreover, companies started to study also ice accretion problems under a wide variety of airfoils geometries. This accident was so crucial that not only did it completely change the focus on ice accretion research, but it also led to a revision of operating procedures on icing conditions, the issue of 18 Airworthiness Directives (ADs) related to in-flight aircraft icing or even some companies transferred their aircraft fleets to southern airports. 1.1.3 Recent years in Ice Accretion research As said, there were two main problems in icing research until 1994: the majority of experimental data and CFD analysis was almost exclusively in NACA 0012. Thus few airfoil shapes were studied, and the knowledge of a new type of ice accretion condition (SLD) was detected. Therefore, since the ATR-72 accident, there has been an interest in testing a new wide variety of wing and airfoil geometries and determining what type of ice shapes are critical for the wing or airfoil performance. Ice accretion research in recent years has focused on improving the understanding of the physical processes that lead to ice formation on surfaces in aircraft and other related areas such as wind turbines or power lines. Those studies have investigated the impact of ice accretion on aircraft performance, especially in lift and drag and other aircraft characteristics. Moreover, it has also been studied the effect of ice accretion in other components or other fields of aviation, such as loads on the structure. However, the main thing that developed during those years was the vast development that numerical methods based on Computational Fluid Dynamics (CFD) have experimented with ice accretion phenomena. Those methods try to simulate and predict the ice accretion over the aerodynamic surface and its implication in aircraft performance. NASA, for example, has developed its analytical ice accretion model that evaluates the ice accretion process in a surface ([ 16 ]) called Lewice. Another analytical ice accretion model example is PoliMIce ([ 5 ]). This two-dimensional model can simulate ice accretion over solid bodies in wet airflow made by Politecnico di Milano. Those models and software have been the base for various research and are only examples of how ice accretion research has increased over the years. For example, research has focused on geometry and flow conditions. Price Raj et al. ([ 12 ]) can be seen in an extensive study of ice accretion and aerodynamic effects of airfoil under SLD conditions, which was the critical condition that triggered the ATR-72 accident. For example, S.G. Pouryoussefi et al. ([ 11 ]) have focused on the influence of airfoil geometry and flow conditions variation. In this study, it has been studied the airfoil NACA 23012 under different conditions of ice shapes (clean airfoil, horn ice, runback ice accretion, etc.) under different conditions and its implication in the airfoil performance, or Johansen et al. ([ 7 ]), focusing on the influence of airspeed and chord length. However, other research is focused on creating an analytical model for ice accretion phenomena such as Lewice or PoliMIce. For example, Chen et al. ([ 3 ]) focus on creating a mathematical model for ice accretion based on unstructured mesh in a 3D model. Dai et al. ([ 4 ]) present a similar way of studying ice accretion phenomena, but this model focuses on the phase-field method. Other investigators have focused on one of the main goals of ice accretion research, which is the improvement of aircraft safety by developing a better method of detecting and preventing ice accretion, using different anti-icing and de-icing methods, such as heating systems or electrothermal systems. Examples of those studies are Arizmendi et al. ([ 6 ]), which made a simulation of a thermal ice protection system in aircraft, using PoliMIce as a numerical model, Xin et al. ([ 17 ]) with the implementation of a coupled moving-boundary method for the analysis of an airfoil electrothermal de-icing system, with a comparison between his method and experimental data. 4Chapter 1. INTRODUCTION 1.2 Objectives and Layout The main objectives of this project are the following: • Make an introduction of ice accretion phenomena and problems related, as well as its main characteristics, types, geometries, etc. • Make an introduction and development of the software used in this project, ANSYS FENSAP ICE, in order to use document as a guide of the main characteristics, modules, problems and tools of this software, so that everybody interested in the study of ice accretion has a guide to start taking their first step. • Prove various methodologies to study this problem in FENSAP ICE, in order to take the optimal choice to approach those analysis. • Make various analysis of the ice accretion phenomena changing various key parameters that defines and has a key role in ice accretion, such as velocity, temperature, Reynolds number, chord length, etc. •Study the importance of the selection of a correct time step in the ice accretion simulations. In order to achieve this goal, this document has been divided into those following chapters apart from this introduction one: • Chapter 2: Ice Accretion. A summary of the ice accretion physics, principal geometries and key parameters. • Chapter 3: Methodology. An explanation of the followed methodology, starting with an introduction of Ansys FENSAP ICE, its principal modules and functionalities. • Chapter 4: Simulations. A summary of all the simulations performed, including its geometries, mesh, set-up and boundary conditions of every case. • Chapter 5: Results. Exposition of the results obtained on the performed simulation, as well as an analysis and discussion of those results. •Chapter 6: Conclusions. 2 ICE ACCRETION This section is devoted to a brief explanation of the physics related to the ice accretion phenomenon. First, a summary of the phenomena will be given, starting with the essential types of atmospheric icing and their main characteristics, as the necessary conditions for those types, such as temperature or geometry. The main ice-shape types and their characteristics, such as geometry or adverse effects on aerodynamics, will be presented. Finally, the most critical parameters that intervene in ice accretion will be explained, as it is needed for the future simulations that will be carried out. 2.1 Types of atmospheric icing Atmospherics conditions for ice accretion can vary depending on temperature, liquid water content (LWC) for clouds, or air velocity magnitude. However, commonly in the literature, it is differentiated into two main types of atmospherics icing conditions, which are rime ice and glaze ice ([8]). Rime ice conditions tend to form at low temperatures, in the order of -10 ◦ C/-15 ◦ C or below, with low liquid water contents at clouds. This produces a high rate of heat loss of the impinging water droplets. Thus, droplets freeze entirely when they contact the aerodynamic surface. Therefore, as droplets freeze upon contact with the surface, no downstream water is over the surface. Ice usually accumulates on the leading edge of the surface or the stagnation point. It generates more regular shapes, following the aerodynamics surface shape, leading to smooth surfaces [15]. Figure 2.1 Rime ice example. 5 6Chapter 2. ICE ACCRETION On the other hand, glaze ice conditions are given in higher temperatures than rime ice, between -3 ◦ C and -10 ◦ C/-15 ◦ C and higher LWCs. In those conditions, there is no high rate of latent heat loss of the water. Thus impinging water droplets do not freeze upon contact with the surface. It leads to a remaining mass of water over the surface, which spreads along the surface downstream before freezing, or even remains unfrozen, creating a film of liquid water over the surface. Those conditions produce larger ice accretions and geometries more different than in rime ice conditions. There is a remaining unfrozen mass of water over the surface, leading to rough surfaces of ice or horn shapes. Figure 2.2 Glaze ice example. However, those conditions are only sometimes fulfilled. Thus, another type of icing condition lies between rime and glaze conditions, called mixed ice. It is just an icing condition in which rime and glaze ice coexists. Thus it has mixed characteristics between both of them. Studying the ice accretion phenomenon can be challenging; even if the simulation methods are the best, there are many variations and casuistry around ice accretion. A wide range of possible temperatures, surface geometries, LWC, or droplets diameter exist. Figure 2.3 Increasing aerodynamic effects vs. Icing conditions and geometry type [2]. 2.2 Ice Accretion geometries In the previous paragraph, it has been seen that atmospheric conditions, even though there are two dominants (rime and glaze ice), can be very different as the range of values for each parameter 2.2 Ice Accretion geometries 7 is extensive. Therefore, it leads also to a wide variety of possible geometries over the surface. However, those types of ice geometries are usually divided into four main ones, as the other possible geometries are usually a mix of those main types [ 2 ]. Fig. 2.4 shows the four main types which are divided in the literature different ice geometries shapes: roughness (Fig. 2.4a), horn ice (Fig. 2.4b), streamwise ice ((Fig. 2.4c) and spanwise-ridge ice (Fig. 2.4d). (a) Ice Roughness. (b) Horn Ice. (c) Streamwise Ice. (d) Spanwise-Ridge Ice. Figure 2.4 Main ice geometry shapes. In the following sub-paragraphs, a brief explanation of the main characteristics of these main types will be explained. For further information, Bragg et al. [ 2 ] describe those types of ice geometry shapes more profoundly and their impact on aerodynamics, as it can be seen qualitatively, for example, in Fig. 2.5. 2.2.1 Ice Roughness Ice roughness is the most straightforward ice shape that can be found, as it occurs at the beginning of the ice accretion process or when the icing exposure is short-lived. Later, during the ice accretion, this shape forms another ice shape, such as ice horn or streamwise ice. Even at the final stages, ice roughness tends to form another ice shape. It is studied separately because, during this initial stage, ice roughness has already altered the surface geometry. Moreover, this ice shape can be found in both glaze and rime ice conditions, as it is a shape that occurs at the initial stages. In this type of shape, three different regions are found near the leading edge or stagnation point: smooth zone, rough zone, and feather region (Fig.2.6). Usually, those ice roughness elements are separated based on the roughness size. This ice shape is considered a 3D shape, as seen in Fig.2.5. Relative to its aerodynamic effect, it will depend primarily on its height, density, and surface location. Thus in the function of those parameters, the ice roughness effect on an airfoil (surface) 8Chapter 2. ICE ACCRETION Figure 2.5 Increasing aerodynamic effects vs. Icing geometry shape. Figure 2.6 Ice Roughness central regions. performance will be higher or lower. The main effect of ice roughness is that it increases noticeable skin friction on the surface. Therefore it directly affects surface performance. Furthermore, it has been studied that higher ice roughness affects boundary later transition and separation, promoting variations in pressure distributions and leading to aerodynamic degradation. 2.2.2 Horn Ice Horn’s ice shape is a geometry usually formed due to glazed and mixed ice generally. It is characterized by its height, its location on the surface, and the angle that respect the airfoil’s chord line, as seen in Fig. 2.7. Moreover, as this ice shape generally has a dominant direction of accretion, which is its height, it is said that it is more a 2D shape than a 3D, unlike ice roughness. Of course, as its height is the most dominant direction of ice accretion, the main aerodynamic effect will be due to the formation of this ice horn. The most critical effect encountered when an ice horn shape appears the formation of bubble zones that appear downstream of the horn, which is a problem related to the separation of the boundary layer that mainly affects the pressure distribution. This is caused because the airflow encounters an irregular and discontinued surface (as it is a horn-ice shape), and its leads to a formation of an adverse pressure gradient near the ice horn shape, leading to a separation of the boundary layer on the top of the horn (Fig. 2.8) 2.2 Ice Accretion geometries 9 Figure 2.7 Horn-ice shape geometry. Figure 2.8 Bubble zone created on an airfoil surface due to the ice horn shape. Moreover, the airflow can transition from laminar flow to turbulent flow in this zone. The direct effect of this transition is the modification of the pressure distribution in the bubble, which becomes virtually constant. However, there is a point after the transition when the reverse airflow velocity magnitude increases, creating vortex flow in the bubble zone. This is critical because it can lead to a pressure recovery and an airflow reattachment on the iced airfoil. Bragg et al. can be seen a study of ice horn shape effect on a NACA 0012, as shown in Fig. 2.9a with the separation streamlines depending of the angle of attack of the incoming airflow and in Fig. 2.9b with the variation of the pressure distribution on a NACA 0012 due to the ice horn shape (a) Pressure distribution on an iced and clean airfoil (NACA 0012). (b) Separation Streamlines. Figure 2.9 Effects of Ice horn shape on a NACA 0012. 16 Chapter 3. METHODOLOGY Following sections will provide an introduction to those modules used, introducing inputs, outputs, governing equations, tips, built-in models and formulas, capabilities, etc. Moreover, Ansys FENSAPICE has more built-in modules, specifically related with turbomachinery, heat transfer problems, remeshing, etc. More info about those modules, which have not been used in this paper, can be found in Ansys FENSAP-ICE’s user guide. Figure 3.4 FENSAP modules. 3.2.1 FENSAP FENSAP is the CFD flow solver module that has FENSAP-ICE implemented. This module is similar to Ansys Fluent and is able to solve 3D complex, compressible, steady and unsteady flows, etc. It has two inputs, geometry and mesh grid, and a output, which is a file that contains flow field solution. The following paragraph will explain governing equations used in FENSAP and its input parameters in order to set-up properly the simulation, as the conditions, boundaries or type of solver selected. Model In this module, the physical model of the problem should be selected, such as the governing equations and the turbulence model. For viscous flow, the option of Navier-Stokes should be selected in the physical model. The governing equations used in ANSYS Fensap for the flow field is modeled by partial differential equations for the conservation of mass (also known as continuity equation), momentum and energy. 3.2 Ansys FENSAP-ICE 17 Figure 3.5 Model section in FENSAP Panel. Starting with the conservation of mass equation, for compressible flows it can be written as: ∂ρa ∂t+ ∇·(ρa Va) = 0(3.1) where: •ρais the fluid density. • Vais the velocity vector field. •tis time. This equation states that the rate of mass entering a system equals the rate of mass leaving the system plus the accumulation or loss of mass within the system. Following with conservation of momentum equation for Newtonian fluid, also known as Navier-Stokes equations, it can be written as: ∂ρa Va ∂t+ ∇·(ρa Va Va) =  ∇·σi j +ρag,(3.2) where: •σi j is the stress tensor. •µis the fluid dynamic viscosity. Those are a set of 3 non-linear equations written in 3D, describing the movement of the viscous fluid. The left side of those equations expresses the acceleration of the fluid. It comprises two components: local acceleration and convective acceleration, respectively. However, the right side of the equations is the summation of effects as hydrostatic pressure, deviatoric stress (both included in 18 Chapter 3. METHODOLOGY the stress tensor component), and body forces. Related with the stress tensor, it can be written as: σi j =−δi j pa+µaδjk∇kvi+δik∇kvj−2 3δjk∇kvk=−δi j pa+τi j,(3.3) where δi j is the Kronecker delta and •pais the static pressure of the fluid. It can bee seen in the Eq. 3.3 the pressure (compression) effect ( −δi j pa ) and the divergence of deviatoric of the stress tensor, τi j. Finally, the final governing equation is the conservation of energy. It states that the energy of the system must be conserved. It can be written as: ∂ρaEa ∂t+ ∇·(ρa VaHa) =  ∇·(κa( ∇Ta)+ viτi j)+ρag· Va(3.4) Where: •Eaequals the total internal energy of the fluid. •Ha=Ea+(pa/ρa)is equal to the enthalpy. •κais the thermal conductivity of the fluid. It is the FULL PDE option integrated into Ansys FENSAP. However, Ansys FENSAP can simplify this equation, assuming the flow is adiabatic. If the assumptions of adiabatic flow are accurate, then the option Constant Enthalpy can be selected. It enables to simplify of the entire energy equation into an algebraic equation, as the following: Ha=H∞=1 2|| Va||2+γ γ−1 pa ρa (3.5) which states that the stagnation enthalpy remains constant along streamlines. In addition, γ=cp/cv is the ratio of specific heats. This option reduce noticeably the computational cost of simulation, but it is necessary to make sure that flow is completely adiabatic in order to have a proper solution. Figure 3.6 Turbulence model integrated in FENSAP Panel. The next option that should be selected in the model is the turbulence model, in case that turbulent flows are studied. FENSAP has implemented a variety of turbulence models, such as: • The Spalart-Allmaras Model: This is a simple model for turbulent flows that only has one equation. • The Low Reynolds kω Model: This is a model that has FENSAP implemented that is made specifically for cases in which Reynolds number is low. 3.2 Ansys FENSAP-ICE 19 • The kω Model: It is a more complex turbulence model which has become popular. It is a two-equation eddy-viscosity model. For further information about those turbulence models implemented in FENSAP, ANSYS FENSAP-ICE User Manual should be reviewed. [1] Figure 3.7 Roughness Model in FENSAP Panel. In the special case for ice accretion, the formation of ice in the leading edge of the airfoil (or surface in general) and therefore, surface roughness, is enough to trigger the transition of laminar flow to turbulent flow. This is why a turbulent model for ice accretion is going to be considered, in order to have a more accurate solution. Next thing to consider in FENSAP model is the surface roughness. It can be selected the option of "No Roughness" for cases of smooth wall. However, as ice accretion problems are studied, there will be a surface roughness due to this ice formation. For that reason, a surface roughness model should be selected. FENSAP has various roughness model to select, as: The next thing to consider in the FENSAP model is surface roughness. It can be selected the option of "No Roughness" for cases of smooth walls. However, as ice accretion problems are studied, the surface roughness will be due to this ice formation. For that reason, a surface roughness model should be selected. FENSAP has various roughness models to select from, such as: •Sand-Grain Roughness: A sand-grain roughness height is needed as a input. • The NASA Roughness Model: The NASA has an empirical correlation for icing that calculates the sand-grain roughness using as a input the Liquid Water Content (LWC). • The Shin et al. Roughness Model: It is a model derived from the NASA Roughness Model. It is just a modification from the empirical correlation used in NASA model. The main different is that in this model, it is need not only LWC, but also the Mean Volumetric Diameter (from now on, MVD) of the droplet in microns. • Personalized Variable Roughness: FENSAP also has the option to introduce an variable roughness based on a input file introduced by the used, and to assign this roughness to all wall or to an specific wall of our problem. As said before, a roughness model should be selected to have a more realistic model for the ice accretion phenomena. In this case, as LWC and MVD are parameters that should be selected as input, and those are boundary conditions, the Shin et al. Roughness model will be selected in this study. Finally, the last thing to complete our model in FENSAP is the selection of body forces that apply to the simulations. In this case, FENSAP permits us to introduce a gravity value or a rotating frame of reference as a type of body force. For our purpose, this will not be considered. 20 Chapter 3. METHODOLOGY Conditions In this section in FENSAP, the reference conditions (flow conditions) and a value of the initial solution should be entered. This section is a crucial factor in the simulation because FENSAP and DROP3D solve equations in a non-dimensional form. Moreover, those reference and flow conditions should be the same for all modules (FENSAP, DROP3D, and ICE3D). Reference values that should be introduced in FENSAP are similar to the ones that should be introduced in Fluent, and those are the following: •The characteristic length (L∞) •The static pressure (p∞) •The velocity magnitude (V∞) •The static temperature (T∞) Figure 3.8 Reference and flow conditions in FENSAP Panel. With those values, FENSAP also calculates other non-dimensional variables or dimensional variables as the following: •Reynolds Number, Re: Re =ρ∞L∞V∞ µ∞ (3.6) •Mach number, M: M∞=|| V∞|| γRgT∞ (3.7) •Adiabatic stagnation temperature, To,∞: To,∞=Ts,∞1+γ−1 2M2 ∞(3.8) To end this section, an initial solution should be entered. This is used as an initial guess to solve the discretized system of equations (remember that Navier-Stokes equations are non-linear, so those are discretized). Therefore the better those values are, the better and more accurate solution will be obtained, and the iterative process will converge faster. 3.2 Ansys FENSAP-ICE 21 Figure 3.9 Initial solution. Velocity Components.. FENSAP can be selected as an initial solution for velocity angles or velocity components imposed in the domain. Good practice in this option is the introduction of boundary conditions as an initial solution. For example, if the Velocity component at the inlet is 10 m/s in the x-component, use this value as an initial solution. It will give a correct first approximation to the problem solution. Boundary Conditions and Domain In this section, the boundary conditions of the problem should be entered. In this panel, an overview of the boundaries of the problem can be seen. Boundaries are the same as the ones introduced in the geometry and mesh (input of FENSAP). However, those boundaries can be changed in this panel. Figure 3.10 Boundary conditions panel.. Furthermore, once boundaries have been checked, boundary conditions should be introduced in every boundary. For this project, the most important ones are inlet, outlet, and wall boundaries, but FENSAP permits us to introduce other boundary conditions such as periodic boundaries, rotors, symmetry planes, etc. Many inlet types, such as subsonic inlet, supersonic, far-field, or Mass flow, can be introduced. In general, parameters that should be introduced are velocity component and temperature, but depending on the inlet type, other parameters such as pressure or mass flow rate can be introduced. There are two things to configure, surface type and heat flux (or temperature). A slip wall (for inviscid flow) or no-slip wall (for viscous flow) should be selected in surface type. For our purpose, the wall type will be no-slip. Then, a heat flux or a different temperature for the wall can be introduced. For the domain panel, this section is only interesting if our problem has more than one domain. Therefore, this panel have not been used in this study. Solver In this panel, the parameters of the solver used should be introduced. The key thing here is the selection of the flow state: steady or unsteady. Depending on this selection, a set of parameters should be introduced. The maximum number of time steps and the CFL number should be set for steady flows. The first one is just an upper limit for the simulation. If the residual convergence level has been reached, doing the maximum number of timestep is unnecessary. However, the CFL number (Courant-Friedrichs-Lewy number) has a deeper meaning. This number is a ratio of how 22 Chapter 3. METHODOLOGY much "information" travels across a computational grid cell in a time step. If the CFL number is large, a convergence problem can be encountered. In order to reduce those problems, a reduction of CFL number should be applied. In FENSAP, a value of 50 is used as a default value. Figure 3.11 Steady Flow panel. The total simulation time and the time step of every iteration should be introduced for unsteady flows. The main condition that FENSAP used to end a time step and to start the following is if the residual values of convergence have been reached, as in the case of steady flows. In addition, FENSAP proposes a mechanism to reduce the instability of the solution found: the introduction of artificial dissipation. This consists of introducing an artificial viscosity to the equations that reduce the oscillations of the solution. FENSAP has many schemes that introduce this concept, such as: •The Streamline Upwind Scheme •Streamline Upwind for Shocks •Central schemes (1st order or 2nd order) Figure 3.12 Example of Streamline upwind scheme.. Output This panel has the option to personalise FENSAP’s output and logs. For example, the solution file can give us the final solution of the simulation, or can save various outputs solutions every "N" iterations. 3.2 Ansys FENSAP-ICE 23 Figure 3.13 Output panel.. The most exciting thing here is that FENSAP can have as an output the resultant force in every surface (solid) of the problem, such as lift, drag, or moment. The critical thing here is the selection of the desired direction for the calculation of lift and drag and the moment reference point if a moment calculation is desired. 3.2.2 DROP3D DROP3D is the droplet and crystal impingement 3D solver module with ANSYS FENSAP integration. As seen before, it works with other CFD solvers, as DROP3D uses the flow solution file obtained from a CFD solver as input. Although any flow data file can be used as an input, working with FENSAP or Ansys Fluent as a main flow solver is easier because Ansys FENSAP has implemented capabilities (such as Multishot) that fit well with Fluent and FENSAP. DROP3D solves differential equations for particle velocity and water concentration. It gives us various parameters as a solution, such as water concentration, droplet velocity, efficiency distributions, impingement zones, etc. Even in this project is used mainly for droplets falling into airfoil or flat plate surfaces, it can be used for other applications such as screens, pollutants, suspender droplets in carrier fluids, etc. As the section before with FENSAP, this section will provide a summary of governing equations and central capabilities of DROP3D, as well as boundary conditions and reference values that should be introduced. Model This section will present the physical model and the governing equations for droplets integrated into DROP3D. Even though more capabilities can be introduced and considered in the physical model, such as the introduction of a vapor model or Supercooled Large Droplet (SLD Conditions), those 24 Chapter 3. METHODOLOGY are outside the project’s scope; thus, it will not be developed in this section. Figure 3.14 Model Section in DROP3D Panel. Starting with governing equations for droplets are very similar to the FENSAP’s equations. Those are the same equations but are augmented by droplets or crystal continuity and momentum equations. Taking into account that subscript "a" refers to the fluid, and "d" subscript refers to droplet, continuity and momentum equations for droplets can be written as follows: ∂(α Vd) ∂t+ ∇·(α Vd) = 0(3.9) where: •αis the mean field value of the water volume fraction. •Vaand Vdare the mean field value of flow and droplet velocity, respectively. ∂(α Vd) ∂t+ ∇[α Vd⊗ Vd] = CDRed 24 Kα( Va− Vd)+α1−ρa ρd1 Fr2g(3.10) where: •dis the mean droplet diameter •Redis the droplet Reynolds number, defined as a follow: Red=ρadVa,∞ µa (3.11) •K is the inertial parameter, defined as a follow: K=ρdd2Va,∞ 18 L∞µa (3.12) 3.2 Ansys FENSAP-ICE 25 •Fr is the local Froude number, defined as a follow: Fr =Va,∞ √L∞g(3.13) As can be seen in the momentum equation, the first term on the right-hand side represents the drag acting on droplets. This term is proportional to Reynolds number and relative droplet velocity. However, it decreases with the inertial parameter. This term is critical in the setup of the physical model in DROP3D. The next thing to introduce in the model section is the droplet drag model, essentially the correlation selection defining the droplet coefficient. There are four droplet drag model correlations that can be applied in DROP3D. Those are: • Water - Default: It is optimal for low-medium values of Reynolds, because is based on an correlation for flow around spherical droplet, but it is known that at higher number of Reynolds, droplet can deform and lose their spherical shape. • Water - Stokes Law:It is optimal for low-medium values of Reynolds, because is based on an correlation for flow around spherical droplet, but it is known that at higher number of Reynolds, droplet can deform and lose their spherical shape.. • Water - Extended Reynolds: It is similar to Water - Default drag model, but it is extended to higher Reynolds number, taking into account the deformation of droplet due to the increase of "Re". •Snowflakes: It is a drag model which is more useful for calculation involving snowflakes For ice accretion problems, it is more accurate to employ the Water - Extended Reynolds, because Reynolds number in aerodynamics applications can reach very high Reynolds numbers. Figure 3.15 Body Force in DROP3D. Finally, the last input in the physical model is the introduction of body forces. It works similarly to introducing body forces in the FENSAP physical model. In ice accretion problems, gravity is the most common body force that should be introduced. The introduction of this body force will depend mainly on MVD. For example, if droplets have a small MVD, it could be assumed that they remain in static atmospheric suspension so that the gravity effect could be neglected. However, not all droplets have the same size. For example, SLD is a case in which droplet MVD is very high. Consequently, gravity will be used and introduced in the physical model as a body force due to the disparity of droplet size found in the atmosphere and the assumption of an uneven distribution of droplet size in this study. Conditions This section is very similar to the conditions panel used in FENSAP, as reference condiitons and a droplet initial solution should be introduced. 32 Chapter 3. METHODOLOGY Figure 3.24 Condition Section in iCE3D Panel. References conditions should be the same as those imposed in FENSAP and DROP3D to guarantee a correct simulation of the ice accretion problem. Thus it is fundamental to make sure that those have not changed. In terms of ice model parameters, as an input must be included parameters such as icing air temperature, which is the temperature at which ice accretion is simulated (in most cases, is the same as the reference air temperature), fluid properties, relative humidity, radiation, and ice density. The last one can be set to a constant value, which is a proper way to compare various simulations, or using a correlation formula, such as the Macklin formula or Jones formula, both include in the options for ice density. Boundary Conditions and Domain In this section, there are no significant boundary conditions that should be included. In order to have a correct ice accretion simulation, enabling icing in the wall’s boundary conditions is the only thing required. It can be other boundary types in the wall, such as sink conditions, but it is not considered in this project’s scope. 3.3 Simulations in Ansys FENSAP-ICE 33 Figure 3.25 Wall’s boundary condition in ICE3D. Solver Ice accretion is naturally unsteady; thus, equations are integrated in time. Therefore, there are only two parameters that must be specified, which are the time step and total time simulation. Figure 3.26 Solver in ICE3D. Related to the time step, which is the incremental change in time for which ICE3D is solving the governing equations, this selection is critical to maintaining a stable run. Suppose a stable or optimal time step is not found. In that case, ICE3D has an option to calculate the optimal time step for every iteration directly, just checking the automatic time step option in the panel. It can avoid some potential problems related to convergence and save computational costs. On the other hand, related to the total time of ice accretion, it is just the time in which the process of simulation stops. Output In a similar way as the other modules (FENSAP and DROP3D), the main thing in this panel is that various solution files can be obtained, just personalising the parameter Time between solution output However, this panel has an exciting option that has been used in this project’s scope: the generation of a displacement grid. This is a new grid that ICE3D generates automatically based on the accumulated ice on the solid boundaries. This grid has the advantage that it can be used in the following simulations in ANSYS FENSAP. Thus it is perfect for simulations using Multi-FENSAP or Multi-FLUENT, as seen in the following paragraph. 3.3 Simulations in Ansys FENSAP-ICE In the following paragraph, different approaches used during this project will be explained. Each approach is perfectly valid and can be used at all. However, depending on the kind of problem that 34 Chapter 3. METHODOLOGY Figure 3.27 Output files and displaced grid in ICE3D. has to be solved, it is more desirable to use one approach or another, even if the basis is the same for each approach (FENSAP/FLUENT into DROP3D into ICE3D) 3.3.1 First Approach: Sequential Method As we can guess from the name, this approach consists of the resolution of the problem sequentially with each module of FENSAP. Therefore, in this approach, in order to get a flow solution, the FENSAP module is used over fluent. Figure 3.28 Flux diagram of Sequential Method.. It is a simple approach since only a proper set-up and inputs are entered in each module. Thus, 3.3 Simulations in Ansys FENSAP-ICE 35 outputs from the module above should be dragged as input for the next module. Accordingly, outputs from FENSAP should be entered as input for DROP3D (drag the solution as an input), and outputs from DROP3D should be entered as input for ICE3D, as both geometry and grid should be in every module. Figure 3.29 Sequential Method in FENSAP. Figure 3.30 Multi-Iteration process with Sequential Method. An increase of the number of module used can be seen. As this approach is sequentially, and the user can only start the next module once the current 36 Chapter 3. METHODOLOGY module is finished, it has the advantage that solutions of every module can be viewed and revised before the beginning of the next module. Thus it is a proper approach when the problem to solve is simple or ice accretion time needs to be longer to do various iterations. For ice accretion simulation time large enough, in which various iteration is needed, this approach is unsuitable because the number of modules that should be introduced in the project is too much, as three modules per iteration must be introduced. Therefore, in the next paragraph, a new approach has been considered to solve the problem with a more extensive ice accretion time simulation, the Sequence option. 3.3.2 Second Approach: Multi-FENSAP Method Figure 3.31 Flux diagram of Multi-Fensap. Multi-FENSAP sequence is a built-in module in Ansys FENSAP, which has integrated the past three modules (FENSAP, DROP3D, and ICE3D) into one module. Inputs and outputs are the same as before for three sub-modules but with the difference that Multi-FENSAP permits introduce several iterations for our problem-based on a chosen parameter (for example, ice accretion time simulation). This module can remake the geometry and the grid of the problem as it is one of the outputs of Multi-FENSAP. Thus, this new grid and geometry are used as input for the next iteration. The main advantage of this approach is that more extensive ice accretion time simulations can be solved accurately, more efficiently, and faster than the previous approach. Thus the simulation time and the number of iterations are only the new inputs compared to the previous approach. However, it is more difficult to check the solution, as the final solution is only finished once all iterations have been solved. At the end of the simulation, all solutions for every iteration can be viewed and checked as before, but only when the complete simulation has finished. In any case, every module’s residuals and iterations can be checked during the simulation. Finally, another problem with this approach is that it has been proven that resolution times with FENSAP (flow solver) are too high compared to other CFD solvers, such as Ansys Fluent. This 3.3 Simulations in Ansys FENSAP-ICE 37 Figure 3.32 Multi-FENSAP Panel. Figure 3.33 Configuration of Multi-FENSAP. Important to select a number of iterations and a parameters to mark the end of each iteration.. consideration and the knowledge that Ansys Fluent has more capabilities than the FENSAP module has led to another approach: Multi-Fluent. 38 Chapter 3. METHODOLOGY Figure 3.34 Residual check during multisequence process.. 3.3.3 Final Approach: Multi-FLUENT Method Figure 3.35 Ansys Workbench 2022 version. As Multi-Fensap, it is a built-in module in Ansys FENSAP that is very similar, but in this case, the flow solver used is Ansys Fluent. Thus DROP3D and ICE3d modules remain the same. The main difference between this approach and the previous one is that it has new inputs: the fluent file case and data (file.cas and file.dat, respectively). 3.3 Simulations in Ansys FENSAP-ICE 39 Figure 3.36 Multi-FLUENT panel with Ansys Fluent. Those files have to be obtained from Ansys Fluent software. In order to get those files, a fluent simulation has to be made. At first sight, this approach seems to be more difficult or cost more time than the previous, but simulation time in Fluent is much lower than FENSAP. Figure 3.37 Multi-FLUENT panel with Ansys Fluent file.cas and file.dat. Subsequently, the rest of the approach works the same way as before, as the remaining modules (DROP3D and ICE3D) are the same. The main advantage of this approach is that time simulation costs are much lower than in FENSAP. Thus, it is preferable to obtain the flow solution with Ansys Fluent. Moreover, if a parametric study of a problem is made, it is easier and faster to solve in Fluent. Obviously, as with the previous approach, in this approach is easy to make a simulation with lots of iterations and extensive time simulations, as this Mullti-FLUENT works equal to Multi-FENSAP, with the difference of the flow solver used. 40 Chapter 3. METHODOLOGY Figure 3.38 Flux diagram of Multi-Fluent. Therefore, as done in the Fensap module, a proper set-up has to be made and then proceed to the simulation. When the simulation is finished, a proper check and review of the flow solution should be made. If everything is correct, case and data files can be exported from Fluent and imported to Ansys FENSAP. 4 SIMULATIONS Th is section will briefly present the geometry and meshing chosen for the flat plate and the cylinder. In addition, first of all, we will show how we have worked in Ansys FENSAP to carry out quasi-2D simulations in a suitable way, using quasi-2D meshes to work with the 3D ice formation model implemented in FENSAP. 4.1 Set-Up: 2-D Simulation in Ansys FENSAP Ansys FENSAP is purely 3D software, but by applying a series of boundary conditions and choosing a suitable mesh, we can perform quasi-2D simulations. This is a great advantage over performing a 3D simulation, as it significantly reduces simulation times, the computational cost decreases considerably, and it provides similar solutions to a 3D model of the problem. In Ansys FENSAP, there is mainly a way to perform a quasi-2D simulation, depending on the geometry and physics of the problem, which considers only a single element in the direction normal to the plane of the airflow. One side of the element should be considered a symmetry plane, and the other one a periodic plane or symmetry. With this consideration, one can obtain a 2D flow, droplet, and ice solution faster and more efficiently from a regular quasi-2D grid without solving a full 3D problem. This quasi-2D mesh grid, according to FENSAP Ice Manual, can be hexa, prism, or a mix between both, but that hybrid mesh can give us convergence problems. Furthermore, it can be assessed if the problem and grid are symmetric. If not, it can give a solution that is not symmetric and unstable and, therefore, invalid. As our problem of the flat plate and the cylinder can be considered symmetric geometrically speaking, the 3D model has been simplified to simplify the problem, minimize computational cost, and ensure optimal convergence in the simulations. 4.2 Flat Plate The following paragraph will present the geometry and grid selected for the flat plate problem. It has been a 2D simplification, as said before, in order to improve convergence solution, optimize computational cost, and decrease simulation time. 4.2.1 Geometry and Boundary Conditions For the flat plate problem, the geometry consists of a volume with six surfaces, which are as follows: •One front side surface, considered as a inlet of the problem. 41 48 Chapter 4. SIMULATIONS Figure 4.12 Geometry of the airfoil problem (II). 4.4.2 Mesh In Figure 4.13 can be seen the final mesh of the airfoil NACA 0012 with a length of one meter. For other lengths, the mesh grid has been similar. It can be seen that a hexahedron mesh has been selected, with inflation near the airfoil’s leading edge, as in a flat plate problem, to capture ice accretion on the airfoil as optimal as possible. Figure 4.13 Mesh of the airfoil problem. (a) Orthogonal Quality. (b) Skewness. Figure 4.14 Mesh metrics for airfoil problem. 4.4.3 Cases For accurate prediction of ice accretion on the airfoil, it is important to study the formation of ice on its surface by analyzing the ice accretion time per iteration. Adopting the same method used for 4.4 Airfoil: NACA 0012 49 the cylinder problem, we can select a suitable ice accretion time per iteration, guaranteeing precise and reliable predictions. Table 4.2 Operating conditions for ice accretion simulation for the airfoil.. Parameter Value Air Speed 104 Temperature (K) 262.04 MVD (micron) 20 LWC (g/m3) 0.55 Icing Duration (min) 7 Angle of Attack (α)Variable First, a simulation with those conditions and α=4 will be made to evaluate the result with the ones obtained in the literature (Pena et al. l, [ 10 ]). Later, a series of simulations with the same conditions can be seen in Table 4.2 will be made, but varying the airfoil’s angle of attack, in order to see how the ice accretion shape changes with the angle of attack. Those simulations (and cylinder simulation) are expected to note the importance of a correct ice accretion time simulation per iteration. 5 RESULTS In this chapter the results of the simulations conducted to analyze ice growth on various surfaces will be discussed. Specifically, we will focus on the elementary case of the flat plate. This analysis involves exploring the variation of density, Reynolds and Eckert numbers, which are crucial factors that affect fundamental parameters such as plate length, temperature, or incident velocity during ice formation. Additionally, we will briefly discuss the analysis of ice formation on a cylinder, emphasizing the importance of choosing the correct time step during simulation. Lastly, we will explore a more interesting study focused on ice formation on an aerodynamic airfoil. 5.1 Flat Plate 5.1.1 Density variation In these simulations it can be considered two boundary conditions for the density of the incident airflow, in order to see if the variation of this parameter has any influence on ice accretion behaviour. First, it has been simulated an ice accretion over the flat plate with a variable density for the incident airflow, considering that density follows ideal gas law. Figure 5.1 shows the result of the mass caught and the total mass caught (both water and ice) for an imposed range of temperature. (a) Mass Caught. (b) Total Mass Caught. Figure 5.1 Mass Caught on a Flat Plate considering ideal gas (variable density) along the chord for a range of temperature during a minute of simulation. Total mass caught is just the mass caught on the flat plat but for the total time of simulation, which in this case is 60 seconds, thus trends shown on those graphics are the same. It can be seen that the 51 52 Chapter 5. RESULTS mass caught over the flat plate surface is higher if the temperature is lower. It is normal because of the higher density of the incident of the airflow. Even though, trends for every temperature are similar to each other, with a high increase of the mass caught at the leading edge of the flat plate, until a point in which the mass caught remain a little bit “constant” (mass caught growth is lesser). It is also due to the fact that with higher temperature, droplets that impact into the flat plate not freeze instantly with the contact with the surface, thus a runback water over the surface is formed and its leave the flat plate at the trailing edge. For example, Figure 5.2 shows ice growth over the flat plate for the same range of temperature. It can be seen that ice growth is naturally higher with lower temperature in the simulation. For lower temperature, ice formation is increased because the temperature is lower than the freezing temperature of the water. At higher temperature, water does not freeze and runback is present. Figure 5.2 Ice accretion on a Flat Plate considering ideal gas (variable density) along the chord for a range of temperature. Figure 5.3 shows the same pictures as before but zoomed near to the leading edge of the flat plate. It can be seen that ice growth is much higher for lower temperature and it is formed closer to the leading edge than with higher temperature, thus mass caught in the leading edge also is lower for higher temperature. It occurs because in those simulations it is used a glaze ice model. This model is more advanced that rime ice, and it can co-exist both ice and runback water over the surface. With rime ice model, all droplets freeze instantly at contact with the flat plate surface. However, even an advanced glaze ice model is used, we can reach a point in which all droplet freeze with contact with the surface, only by lowering the temperature sufficient for incident air flow. In fact, if a comparison is made between Figure 5.2 and Figure 5.2, it can be seen that total mass caught for a temperature equal to 100 K is the same as the total ice growth over the flat plate surface. This solution is the same as if we would have considered a rime ice model. 5.1 Flat Plate 53 (a) Mass Caught. (b) Ice Growth. Figure 5.3 Mass caught and ice growth over the flat plate for variant density near the leading edge.. Table 5.1 Ratio between ice accretion over the flat plate and total mass caught. Temperature Ratio 267 0.07494942 261 0.142970082 255 0.206442597 252 0.236620035 249 0.265808368 246 0.294073881 200 0.636905254 150 0.883786544 100 1 In Table 5.1 can be found the ratio between total mass caught and ice growth over the surface. There is a point in which the glaze and rime ice models give us the same results, but for higher temperatures; it must be noted that this ratio is not equal to one. It means that the rest of the mass caught over the flat plate is liquid (film water over the flat plate surface or runback water). The same analysis has been done for the same conditions of ice simulations but using a constant density value. In Figure 5.4, the mass and total mass caught over the surface of the flat plate. The first difference that can be seen is that the mass and total mass caught is independent of the temperature, as the density is the same. However, it is different for ice growth on the flat plate. As can be seen in Figure 5.5, the trend showed before for ideal gas density is the same as the trend for constant density, in which ice accretion is naturally higher when the temperature is low. 54 Chapter 5. RESULTS (a) Mass Caught. (b) Total Mass Caught. Figure 5.4 Mass Caught on a Flat Plate considering constant density along the chord for a range of temperature during a minute of simulation. Figure 5.5 Ice accretion on a Flat Plate considering constant density along the chord for a range of temperature. Another interesting point in that analysis is the comparison of the ratios between ice accretion and total mass caught on the flat plate for both constant and variable density. Figure 5.6 shows this, in which the blue line shows the boundary of the ratio, which is that total mass caught is the same as ice growth (case of rime ice), where all fallen droplet freeze on the surface, and red line show the freezing temperature for water. 5.1 Flat Plate 55 Figure 5.6 Ratio between total mass caught and ice accretion on the flat plate with glaze simulations. Black line shows the ratio obtained for variable density in Table 5.1 and the cross the ratio for constant density. As can be seen, the ratio at the same temperature for both cases (constant and variable density) are the same. As a consequence, we can say that this ratio is independent of density variation thus the consideration of a constant density for those analysis can be done, in order to minimize ice accretion variables and to realise a parametric study for more significant variables such as temperature, velocity or chord length through the variation of Reynolds and Eckert numbers. 5.1.2 Variation of the Reynolds number This section presents the outcomes of a parametric analysis of the Reynolds Number. The analysis involved sweeping the Reynolds number over several orders of magnitude by adjusting the incident flow velocity on the plate while keeping the chord length constant. In the subsequent analysis, the chord length will be altered with the variation of the Eckert number. The study was conducted for two temperature values and discovered that the trends for a more comprehensive temperature range are comparable. Therefore, performing the analysis for two temperature values is representative. (a) T = 200 K. (b) T = 150 K. Figure 5.7 Mass Caught on a Flat Plate considering constant density along the chord for a range of Reynolds number during a minute of simulation. Figure 5.7 shows the Mass caught on a flat plate with varying Reynolds numbers for both temperatures (T=150 K and T=200 K). It can be seen that for both temperatures, the tendencies are the 56 Chapter 5. RESULTS Figure 5.8 Mass caught comparison between both temperature. same. In fact, in Figure 5.8, it can be seen that the Mass caught is the same for both values; thus, for a given Reynolds number, the Mass caught over the surface of a flat plate is independent of the temperature. As will be seen next, the difference will be given in terms of the proportion between ice and liquid water over the surface. It is noticeable that at higher Reynolds numbers, the weight on the flat plate’s surface tends to accumulate more towards the back of the plate. The Reynolds number is the ratio of inertial and viscous forces, meaning that if the Reynolds number is more significant, the velocity of the airflow is higher, and the inertial forces become more dominant over the viscous forces. As a result, mass accumulation at the leading edge slows down and moves closer to the trailing edge as the airflow velocity (and Reynolds number) increases. This behavior is evident in ice and water film on the flat plate’s surface, as demonstrated in the examples below. (a) T = 200 K. (b) T = 150 K. Figure 5.9 Ice growth on a Flat Plate considering constant density along the chord for a range of Reynolds number during a minute of simulation. Figure 5.9 shows ice growth over the flat plate’s surface for the same temperature values. For both values, we can observe the sale trends as in Figure 5.7 or Figure 5.8, as if the incident airflow 5.1 Flat Plate 57 Figure 5.10 Ice accretion comparison between both temperature. is higher, ice growth formed onto the flat plate is lower. Because of the higher inertial force of the incident airflow, part of the mass falling on the plate does not remain on the plate. It falls off the trailing edge, reducing the final mass caught over the surface and, therefore, the ice accretion. If we lower the Reynolds number sufficiently, we can see that there comes the point where the maximum mass of ice accumulated on the plate reaches its limit and even remains constant along the length of the plate. In this case, there are differences between ice growth over the surface for both temperature values, as its accretion naturally depends on the temperature. In Figure 5.10, we can see ice accretion on the flat plate for both temperatures for the same Reynolds number, in which the solid line is related to 150 K and the dotted line is related to 200 K. Naturally if the temperature is higher, ice accretion is lower, for every value of Reynolds considered. However, it can be seen that trends followed by both temperatures are very similar. The difference between both values of ice accretion for both temperatures is just the amount of liquid water that has not frozen over the surface. It can be seen in Figure 5.11, which shows the thickness of the water film that has been formed over the flat plate for both temperatures. As to the figures above, trends for both values of temperature are similar to each other. Contrary to ice accretion, the film thickness is higher for higher values of temperature, as there are more incident droplets that do not freeze on contact, forming a film of water over the surface. This can be seen perfectly in Figure 5.12, which compares film thickness for both temperature values, where we can observe those trends clearly. 5.1.3 Variation of the Eckert number In this section, we will discuss the findings from our analysis of the parametric variation of the Eckert number. Our analysis involved testing two different values of the Reynolds number and a range of temperatures. To manipulate the Eckert number, we adjusted the temperature and velocity values of the incident air. We ensured that the Reynolds number remained constant across all cases by varying the length of the plate. Please note that the length of the plate is not unitary for this case, 64 Chapter 5. RESULTS (a) Case 1. Eckert equal to 8.7e-6.. (b) Case 2. Eckert equal to 5.4e-5.. (c) Case 3. Eckert equal to 2.2e-4. (d) Case 4. Eckert equal to 8.7e-4.. (e) Case 5. Eckert equal to 5.4e-3.. Figure 5.19 Ice Growth simulations cases divided by each Eckert number. AA = 200K, CC = 225 K and BB = 250 K. 5.2 Cylinder 65 5.2 Cylinder Within this section, the obtained results of the analysis conducted for the cylinder are available for observation. The purpose of the analysis was to confirm the correlation between surface ice formation and the time step selected in the simulation. The cylinder was chosen for this investigation due to an existing case study in the literature on ice formation on its surface. In the original study, a simulation was performed for 30 minutes of extended exposure, but only in a solitary time step. Our aim was to replicate these conditions as accurately as possible, based on information gathered from the literature. We were able to reduce simulation times between each iteration through multiple simulations, ultimately reaching a point where solutions converge. To achieve this, we decreased the time-step of the iteration. Table 4.1 presents the operating conditions we used for the cylinder simulations, which were selected because they are identical to those in Sokolov et al. ([ 13 ]), who analyzed the aerodynamic force on an iced cylinder. Figure 5.20 Ice accretion on a cylinder with no iteration (t = 30 min).. Based on information from Sokolov ([ 13 ]), we conducted a simulation to replicate operating conditions. Our findings, as depicted in Figure 5.20, revealed that the ice accretion on the cylinder closely matched that of [ 13 ]. Although our peak ice accretion was marginally lower than the literature’s report, this may be attributed to variances in the solver’s internal parameters. In Montminy et al. ([ 9 ]), it was observed that the time interval utilized in ice accretion simulations can significantly impact the resulting ice accretion. They utilized the LEWICE ice accretion solver and found that the ice shapes achieved for different time intervals varied, as depicted in Figure 5.21. It is noteworthy that modifying the simulation’s time interval can yield a distinct solution for the ice accretion. From the figures presented (Figure 5.22 and Figure 5.21), it is clear that the amount of ice that accumulates on a cylinder is influenced by the time step used in simulations. It is crucial to choose an appropriate time step value to ensure accurate results. A full simulation with no 66 Chapter 5. RESULTS Figure 5.21 Ice shapes obtained in [9] for various time interval considered. Figure 5.22 Ice accretion obtained on a cylinder for a various values of time-step interval.. intermediate iteration produces a solution that differs significantly from the final one. When using a shorter time step, the solution is narrow and peaked, while a low time step results in very similar ice accretion, indicating that there is a value at which the solution converges. Therefore, when conducting simulations, it is essential to consider the time step interval to obtain the desired results. 5.3 Airfoil: NACA 0012 67 5.3 Airfoil: NACA 0012 In this part of the report, the findings from the analysis of the airfoil, specifically the NACA 0012, are presented. One of the primary objectives was to verify the correlation between the time step iteration for ice accretion and the final ice accretion shape on the airfoil, similar to the approach used for the cylinder. Additionally, we obtained the ice accretion shape for various angles of attack to understand how the ice accretion varied for different angle of attack values. We utilized the same boundary conditions as illustrated in Table 4.2 and Pena et al. ([10]) for these simulations. Figure 5.23 Ice accretion obtained on an airfoil for two values of ice accretion time step iteration. In Figure 5.23, we can observe the ice accretion on an airfoil obtained in Ansys FENSAP under the conditions above for two types of simulations: the first one, a simulation of the full ice accretion time (420 seconds) in a unique simulation, and the second one, introducing three time-steps, one each t = 140 seconds. Noticeable differences in the final ice accretion shape can be seen for both simulations. For example, for the full simulation, a smoother surface geometry can be observed than in the other one. This is because FENSAP calculates the ice formation’s geometry on the clean airfoil, as the simulation has been made in a single iteration. For the second case, the ice formation has been calculated in every iteration on the grid obtained in the previous iteration, giving a different final ice accretion shape in this case. This can be visualized in Figure 5.24, in which the ice accretion shape can be observed for every iteration made. The comparison of different solutions for analyzing ice formation on an airfoil is shown in Figure 5.25. Pena et al.’s practical solution offers a good starting point for analysis. However, FENSAP’s multi-iteration solution is much closer to both the experimental case and Pena’s result than the full simulation in FENSAP. One interesting feature of the ice formation is the horn at the airfoil’s leading edge, which is more pronounced in the experimental case than in FENSAP’s solution with a single iteration. Although the simulation captures the formation of the horn, it is smaller than the actual one. Nevertheless, the multi-iteration FENSAP solution exhibits a much more prominent horn, similar to the one observed in reality. 68 Chapter 5. RESULTS (a) t = 140s. (b) t = 280s. (c) t = 420s. Figure 5.24 Ice accretion shape on the NACA 0012 for each iteration of the simulation.. Figure 5.25 Comparison of ice accretion obtained on an airfoil between the solution obtained by Pena et all ([10], experimental and Ansys FENSAP.. Upon examining the experimental case, it becomes apparent that the surface affected by ice is quite uneven, similar to the results achieved by Pena and the multi-iteration FENSAP solution. The single iteration solution produces a smoother ice accretion shape, but it is still not practical. A closer look at Figure 5.26 reveals the same solution as before, but with a closer focus on the leading 5.3 Airfoil: NACA 0012 69 edge of the airfoil. This figure provides a clearer view of the leading edge horn and the rough ice accretion shape. It is evident that the amount of time allocated for ice accretion and the number of iterations are critical factors in simulating ice accretion. Figure 5.26 Zoom in the leading edge of the airfoil. Figure 5.27 Ice accretion on NACA0012 depending on Angle of Attack. The provided figure depicts various ice accretion geometries obtained while simulating different angles of attack on an NACA 0012 airfoil. It is evident that a higher angle of attack results in ice accretion starting from a deeper point, while a lower angle of attack leads to ice accretion starting near the leading edge of the airfoil. Furthermore, Figure 5.28 showcases some final ice accretion examples for these cases. 70 Chapter 5. RESULTS (a) AoA = -2o.(b) AoA = 2o. Figure 5.28 Example of ice accretion on NACA 0012 for two values of AoA. 6 CONCLUSIONS As seen during this project, the ice accretion problem is complex. It is a phenomenon that involves numerous crucial parameters and factors, making the study of this problem difficult and extensive. That is why researchers and organizations such as NASA or universities have studied and analyzed it throughout the last century. It started in the early 30s with experimental studies in laboratories with clean and iced airfoils and later more extensively with the introduction of CFD calculations. It permitted the study of the ice accretion problem in a more complex and faster way, and many studies related to the airfoil and wing geometry started to develop. As we have previously stated, our project has provided a comprehensive guide for individuals who are new to ice accretion research. This guide covers both the general physics associated with ice accretion and CFD software. Our team has conducted a series of parametric simulations over various geometries to gain insight into how key parameters affect ice formation over different surfaces. We are confident that this guide and the results from our simulations will be an invaluable resource for those interested in pursuing research in this field. The purpose of this paper is to offer guidance to those who are new to the ice accretion field. Our approach involves utilizing Ansys FENSAP ICE, a CFD software specifically designed for ice accretion problems. We have provided an overview of the primary modules integrated into the software, such as FENSAP, DROP3D, and ICE3D, as well as the key parameters and input/output requirements for conducting an ice accretion simulation. We have also evaluated different methodologies for conducting an ice accretion simulation, including replacing FENSAP with Ansys FLUENT as the primary module for air solution. Our analysis includes testing all of these methodologies. Lastly, we have included a summary of common ice accretion geometries and their impact on airfoil aerodynamics. Within the second main objective, in this paper it has been carried out several simulations varying some crucial parameters such as Reynolds and Eckert number, temperature and chord length. An extensive simulation work has been done within this project, where thousands of simulation hours have been spent, including software testing, refinement, set-up of the program and the simulations performed, both correct and failed simulations. It has been an extensive an prosperous work where its principal parametric study gives us those principal findings: • The temperature significantly influences the ice accretion’s final shape and quantity. For higher temperatures, ice accretion formed over the surface is lower, but the film thickness of liquid water over the surface has a greater thickness. • In our study, we carefully examined the impact of Reynolds number on the shape of ice accretion. By manipulating the airflow velocity while keeping the chord length constant, we 71 72 Chapter 6. CONCLUSIONS observed a significant influence of Reynolds number on the kinetic energy of the flow. As a result, the final shape of the ice accretion was greatly affected. Our findings showed that lower velocities reached the peak of mass and ice caught on the surface before higher Reynolds. Conversely, higher Reynolds led to an increase in kinetic energy, causing the mass caught to fall off the trailing edge and a consequent reduction in the mass caught and ice growth. • The Eckert variation gave us insight into the ice shape for various values of chord length. The trends obtained for every Eckert are very similar between them, founding that at lower Eckert, ice accretion is formed nearer the leading edge than at higher Eckert. • The simulation of the NACA 0012 airfoil and the cylinder and its comparison with the bibliography have shown that a correct selection of time steps and ice accretion time simulation is crucial for a correct estimation of the ice geometry over the surface for extensive time simulations, showing a better relationship with the practical solution for lower ice accretion time for every time step. We conducted a thorough analysis of the project objectives, providing valuable insights into critical parameters for ice formation. Despite facing computational constraints and time limitations, I made calculated assumptions and simplifications to ensure the project’s success. Although the selected mesh had limitations, it was able possible to refine it to the best of my ability for complex geometries. These findings mark an excellent starting point for future research, and I am confident that they will guide the team towards even greater success. Through an in-depth analysis of various simulations, including a wide range of values and parameters, I have achieved my objectives. By conducting parametric simulations and referencing current literature, I have laid a strong foundation for future, more comprehensive analyses that consider additional parameters, geometries, and devices such as de-icing and anti-ice systems. My ultimate aim is to make significant contributions to ice accretion research and publish articles based on this document. Working collaboratively with my two supervisors, I am committed to achieving this goal. We already have a research article in progress. List of Figures 1.1 ATR-72 Aircraft 2 2.1 Rime ice example 5 2.2 Glaze ice example 6 2.3 Increasing aerodynamic effects vs. Icing conditions and geometry type [2] 6 2.4 Main ice geometry shapes 7 2.5 Increasing aerodynamic effects vs. Icing geometry shape 8 2.6 Ice Roughness central regions 8 2.7 Horn-ice shape geometry 9 2.8 Bubble zone created on an airfoil surface due to the ice horn shape 9 2.9 Effects of Ice horn shape on a NACA 0012 9 2.10 Key parameters for ice accretion problem 11 3.1 Procedure followed in CFD Simulations. 14 3.2 Ansys logo [1] 15 3.3 Summary of basic ANSYS-FENSAP scheme 15 3.4 FENSAP modules 16 3.5 Model section in FENSAP Panel 17 3.6 Turbulence model integrated in FENSAP Panel 18 3.7 Roughness Model in FENSAP Panel 19 3.8 Reference and flow conditions in FENSAP Panel 20 3.10 Boundary conditions panel. 21 3.11 Steady Flow panel 22 3.12 Example of Streamline upwind scheme. 22 3.13 Output panel. 23 3.14 Model Section in DROP3D Panel 24 3.15 Body Force in DROP3D 25 3.16 Condition Section in DROP3D Panel 26 3.17 Droplet initial solutions 27 3.18 Inlet boundary conditions in DROP3D 27 3.19 Solver panel of DROP3D 28 3.20 Output panel of DROP3D 28 3.21 Ice accretion physics scheme 29 3.22 Icing model in ICE3D 31 3.23 Body forces in ICE3D 31 3.24 Condition Section in iCE3D Panel 32 73 80 Bibliography [13] Pavlo Sokolov and Muhammad Shakeel Virk, Aerodynamic forces on iced cylinder for dry ice accretion – a numerical study, Journal of Wind Engineering and Industrial Aerodynamics 206 (2020). [14] Spencer J. Stebbins, Eric Loth, Andy P. Broeren, and Mark Potapczuk, Review of computational methods for aerodynamic analysis of iced lifting surfaces, 11 2019. [15] K Szilder and W Yuan, In-flight icing on unmanned aerial vehicle and its aerodynamic penalties. [16] William B Wright, User manual for the nasa glenn ice accretion code lewice version 2.0, 1999. [17] Miao Xin, Guo Zhong, and Yihua Cao, Numerical simulation of an airfoil electrothermaldeicing-system in the framework of a coupled moving-boundary method, Fluid Dynamics and Materials Processing 16 (2020), 1063–1092. Index 81