A study of droplet evaporation coupling model based on Eulerian method

<!-- metadata commented in wiki content


==A study of droplet evaporation coupling model based on Eulerian method==

Yan TAN<sup>1,2</sup>,  Jiupeng SHANG<sup>3</sup>,  Yuyou ZHAN<sup>3</sup>

1. Aeroengine maintenance training center, Civil Aviation Flight University of China, Guanghan Sichuan 618307, China

2. Key Laboratory of Icing and Anti/De-icing, China Aerodynamics Research and Development Center, Mianyang Sichuan 621000, China

3. Aeronautical Engineering Institute, Civil Aviation Flight University of China, Guanghan Sichuan 618307, China
-->

==Abstract==

Research on engine icing is a hot topic among the world. Different from the aircraft wing or airframe icing, the evaporation phenomenon in the internal flow field has a great influence on the engine icing. Moreover, the thermodynamic coupling between droplets and flow field is not available in current particle trajectory calculations, or only for one-dimensional situation. Therefore, a three-dimensional droplet trajectory calculation model based on Eulerian method is used to demonstrate the thermodynamic coupling between droplets and flow field. The model was verified by NRC small engine icing wind tunnel test data and the flow field evolution is obtained which cannot be obtained by the one-dimensional coupling model. In the meanwhile, the effects of different initial LWC, relative humidity and MVD on the internal flow evaporation were studied, and the trends of droplets and flow field affected by evaporation were obtained. The numerical method in this paper can provide guidance for the subsequent research on engine icing.

==Key Words==

<span id='_GoBack'></span>Eulerian method, trajectory calculation, icing, numerical simulation

==Nomenclature==

{| style="width: 100%;border-collapse: collapse;" 
|-
| LWC
| liquid water content
| MVD
| median volume diameter
|-
| RH
| relative humidity
| <math>Sh</math> 
| Sherwood number
|-
| <math>d</math> 
| diameter
| <math>D</math> 
| diffusivity coefficient
|-
| <math>\rho </math> 
| density
| <math>Re</math> 
| Reynolds number
|-
| <math>Sc</math> 
| Schmidt number
| <math>u</math> 
| velocity of flow field
|-
| <math>v</math> 
| particle velocity
| <math>P_{vap}</math> 
| partial vapor pressure
|-
| <math>P_{sat}</math> 
| saturation vapor pressure
| <math>M</math> 
| molar mass
|-
| <math>{\alpha }_p</math> 
| LWC
| <math>V</math> 
| volume fraction
|-
| <math>\phi </math> 
| volume fraction
| <math>{\alpha }_{\mbox{va}p}</math> 
| vapor concentration
|-
| <math>D_{eff}</math> 
| effective mass diffusion coefficient
| <math>e_p</math> 
| particle internal energy
|-
| <math>e</math> 
| air internal energy
| <math>k</math> 
| heat transfer coefficient
|-
| <math>{\tau }_{ij}</math> 
| viscous shear stress tensor
| <math>q</math> 
| heat transfer rate
|-
| <math>\mu </math> 
| dynamic viscosity
| 
| 
|}

==Subscripts==

{| style="width: 100%;border-collapse: collapse;" 
|-
|  style="vertical-align: top;width: 11%;"|<math>p</math> 
|  style="vertical-align: top;"|particle
|  style="vertical-align: top;width: 8%;"|<math>a</math> 
|  style="vertical-align: top;"|air
|-
|  style="vertical-align: top;width: 11%;"|<math>l</math> 
|  style="vertical-align: top;"|laminar
|  style="vertical-align: top;width: 8%;"|<math>vap</math> 
|  style="vertical-align: top;"|vapor
|-
|  style="vertical-align: top;width: 11%;"|<math>conv</math> 
|  style="vertical-align: top;"|convection
|  style="vertical-align: top;width: 8%;"|<math>evap</math> 
|  style="vertical-align: top;"|evaporative
|}


==1. Introduction==

Icing during flight has always been an important threatening to flight safety. An extensive amount of research over the world have carried out and summarized the meteorological conditions of aircraft icing [1]. However, a series of turbofan engine malfunctions characterized by stalling or power loss during high altitude cruise in the past decades [2-4]. As the research moves along, more and more evidences show that the icing at high altitude is the root cause of the above problems. The engine icing is caused by not only supercooled droplets but also ice crystals, which greatly expands the scope of the original aircraft icing meteorological conditions.

At present, engine icing has attracted extensive attention all over the world. The related research has been carried out and some achievements have been made in recent ten years. Under certain assumptions, engine icing process have been simulated by using numerical method and engine icing mechanism have preliminarily explored [5-9]. However, Due to lack of experimental data, the numerical simulation suffered great limitation. Considering experimental research, NRC (National Research Council) of Canada and NASA (National Aeronautics and Space Administration) of America have made great contributions to the experimental study of engine icing mechanism. NRC carried out local flow simulation based on research altitude test facility [10-12]. NASA have established the propulsion system lab (PSL) for full-scale engine icing test [13]. In the experimental study, it is found that the evaporation or sublimation of particles will have a certain impact on the flow field, especially for the internal flow such as engine icing. Therefore, in the study of engine icing, the influence of evaporation on the flow field environment must be considered, and the evaporation phenomenon needs higher requirements for data acquisition.

Although some of the existing particle trajectory calculation models consider the changes of particles which caused by evaporation and sublimation, they are not coupled with the initial flow field [5,14-15], which cannot be used for wind tunnel data verification. For this reason, some scholars have proposed a new coupling model. Davison of NRC proposed an evaporation model for the engine ice wind tunnel. This evaporation model considered the heat and mass transfer and has the capability to run particle distributions for a real time solution [16]. The model proposed by Bartkus takes into account the liquid-solid-gas three-phase thermodynamic coupling and obtains the verification of PSL wind tunnel test results [17]. However, these two models are one-dimensional models, which can quickly calculate the results, but cannot be used in the study of 3D engine icing problem with complex geometry.

A three-dimensional particle trajectory calculation model based on Eulerian method is established to obtain the thermodynamic coupling of particles and flow field without considering the solid-liquid mixing phase. NRC small ice tunnel in Ottawa was numerically studied based on the three-dimensional model. And the effects of different initial LWC, relative humidity and MVD were studied.

==2. Description of the Model  ==

The air flow field is first calculated by the Spalart-Allmaras turbulence model [18], and the calculated results are used as known conditions for particle trajectory calculation.

According to Fick's Law [19-20], the evaporation rate of single particle resulting from evaporation is defined as:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s_p=\frac{Sh\pi d_p{\rho }_aD_l\left({\omega }_{\infty }-{\omega }_{surface}\right)}{\pi d_p^3/6}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Sh=2+0.6Sc^{1/3}{Re}_p^{1/2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{Re}_p=\frac{{\rho }_a\vert u-v\vert d_p}{{\mu }_a}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Sc=\frac{{\mu }_a}{{\rho }_aD_{vap}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}


<math>{\omega }_{\infty }</math> and  <math>{\omega }_{surface}</math> are the mass fraction close to and on the particle surface, respectively (Eqn. (5)).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\omega }_{\infty }=\frac{M_{vap}}{M_a}\frac{P_{vap}}{P_a},{\omega }_{\infty }=</math><math>\frac{M_{vap}}{M_a}\frac{P_{sat}}{P_a}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}


The molar mass is 18 g/mol for water vapor and 29 g/mol for air.

Combining eqn. (5) with eqn. (1), results in eqn. (6) for the evaporation rate of single particle (Eqn. (6)).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s_p=\frac{Sh\pi d_pD_l}{\pi d_p^3/6}\cdot \frac{M_{vap}}{M_a}\cdot \left(P_{vap}-\right. </math><math>\left. P_{sat}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}


It is assumed that in a finite control volume, all particles have the same physical properties, such as the same density and the same volume. Therefore, the LWC in a control volume can be obtained by integrating the control volume (Eqn. (7)).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\alpha }_p\mbox{=}\frac{{\int\!\int\!\int }_V{\rho }_pdV_p}{V}=</math><math>\frac{V_p}{V}\frac{{\int\!\int\!\int }_V{\rho }_pdV_p}{V_p}=</math><math>\frac{V_p}{V}\frac{{\rho }_pV_p}{V_p}={\phi }_p{\rho }_p</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}


Therefore，eqn.(8) is the solution of the evaporation rate in a finite control volume.

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>S_p=\frac{{\int\!\int\!\int }_Vs_pdV_p}{V}=\frac{V_p}{V}\frac{{\int\!\int\!\int }_Vs_pdV_p}{V_p}=</math><math>{\phi }_ps_p=\frac{{\alpha }_ps_p}{{\rho }_p}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}


The particle and vapor mass control equations are given by eqns. (9) and (10) respectively.

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial {\alpha }_p}{\partial t}+\frac{\partial \left({\alpha }_pv_j\right)}{\partial x_j}=</math><math>-S_p</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial {\alpha }_{\mbox{va}p}}{\partial t}+</math><math>\frac{\partial \left({\alpha }_{vap}u_j\right)}{\partial x_j}=</math><math>\frac{\partial }{\partial x_j}\left(D_{eff}\frac{\partial {\alpha }_{\mbox{va}p}}{\partial x_j}\right)+</math><math>S_p</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}


It is assumed that the velocity and temperature of water vapor and air flow field are consistent, so it is not necessary to establish the momentum and energy governing equations of water vapor, but the momentum conservation equation of droplet particles needs to be re-established and shown by eqn. (11).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial \left({\alpha }_pv_i\right)}{\partial t}+</math><math>\frac{\partial \left({\alpha }_pv_iv_j\right)}{\partial x_j}=</math><math>{\alpha }_p\frac{3{\rho }_aC_D}{4{\rho }_pd_p^2}\vert u-</math><math>v\vert \left(u_i-v_i\right)-v_iS_p</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}


Eqn. (12) is the energy governing equation of droplet which is combined with the energy equation in flow field again results in eqn. (13).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial \left({\alpha }_pe_p\right)}{\partial t}+</math><math>\frac{\partial \left({\alpha }_pe_pv_j\right)}{\partial x_j}=</math><math>{\alpha }_p\left(q_{evap}+q_{conv}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial \left(\rho e\right)}{\partial t}+\frac{\partial \left(\rho eu_j\right)}{\partial x_j}=</math><math>\frac{\partial }{\partial x_j}\left(k\frac{\partial T}{\partial x_j}\right)-</math><math>p_a\frac{\partial u_j}{\partial x_j}+{\tau }_{ij}\frac{\partial u_j}{\partial x_j}-</math><math>{\alpha }_pq_{evap}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}


==3. Case Study ==

==3.1 Validation==

Based on the NRC small ice wind tunnel in Ottawa shown in Figure 1, a quarter symmetrical mesh model is built with a total of 20 million tetrahedral meshes.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_TAN_199218120-image44.png|486px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 1. NRC small engine icing wind tunnel geometry (units: m)[16]</div>

One of the test conditions in paper [16] is as follows: air flow velocity at spray mast is 4.8 m/s; initial temperature is 273.15 K; relative humidity is 31%; particle LWC is 0.99 g/m3; particle MVD is 15.1 μm; pressure at engine inlet is 101.35 kPa. From the experimental results, the MVD of the particles at the engine inlet is 22.8 μm. MVD of particles increases in the strong evaporation environment, which seems to violate the law of mass conservation, but the results are in line with the actual engineering. This is because the particles ejected by the spray master are not exactly the same size, with an average diameter of 15.1 μm. In the process of evaporation, the diameter of all particles becomes smaller, but the size of smaller particles becomes smaller after evaporation, and then cannot be detected by the [javascript:; measuring] [javascript:; equipment], thus indirectly increasing the MVD value at the engine inlet [16]. However, there is no problem of data acquisition in the numerical simulation. In order to match the relative humidity at the engine inlet, the boundary conditions are the same as those of the experiment, but particles with same diameter which are 20μm are used as the initial condition, and the initial temperature of the particles is defined as 293.15 K.

It is shown that the calculated results are quite consistent with the experimental results from the comparison of parameters in Table 1. In addition, from the comparison of temperature fields in Fig. 2, it can be found that the temperature of the background flow field at the engine inlet decreases significantly due to evaporation.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Table 1. Comparison of parameters at engine inlet </div>

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Experimental results
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Numerical results 
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Relative error((Numerical results - Experimental results)/ Experimental results)
|-
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|Speed (m/s)
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|35.8
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|35.5
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|-0.84%
|-
|  style="text-align: center;vertical-align: top;"|Relative humidity (%)
|  style="text-align: center;vertical-align: top;"|57
|  style="text-align: center;vertical-align: top;"|57.78
|  style="text-align: center;vertical-align: top;"|1.36%
|-
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Temperature (K)
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|271.55
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|271.21
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|-0.13%
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_TAN_199218120-image45.png|600px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 2. Comparison of temperature (up: Pre- evaporation model coupling; down: Post -evaporation model coupling)</div>

==3.2 Effect of LWC==

In order to study the effect of LWC on the variables at engine inlet, the simulations take the above test conditions as reference, and only change the LWC value at spray mast and keep the remaining conditions unchanged. The calculation cases are shown in Table 2.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Table 2. Cases of different LWC</div>

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Case 11
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Reference Case
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Case 12
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Case 13
|-
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|LWC (g/m<sup>3</sup>)
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|0.5
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|0.99
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|2
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|4
|-
|  style="text-align: center;vertical-align: top;"|RH (%)
|  style="text-align: center;vertical-align: top;"|31
|  style="text-align: center;vertical-align: top;"|31
|  style="text-align: center;vertical-align: top;"|31
|  style="text-align: center;vertical-align: top;"|31
|-
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|<span id='OLE_LINK1'></span>MVD (μm)
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|20
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|20
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|20
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|20
|}


Figure 3 shows the comparison of results for different LWC, in which the comparison data are at the central axis of the engine icing wind tunnel.

The initial LWC of case 13 is 8 times that of case 11, but the loss of LWC at engine inlet of case 13 is only 4.4 times of that of case 11 (Fig. 3e). This shows that for the same particle diameter, the total evaporation is directly proportional to the initial LWC, but the evaporation of a single particle is proportional to the LWC inversely, which results in the MVD of case 13 being larger than that of the other three cases (Fig. 3d), and the final MVD of case 13 is almost 3.3 times that of case 11. Moreover, the changes of temperature and relative humidity in the background flow field are mainly influenced by the total evaporation. The larger total evaporation, the lower temperature of the background flow field (Fig. 3a), and the greater relative humidity at engine inlet (Fig. 3f).

In addition, the evaporation capacity of a single particle has a great influence on the particle MVD, but less than 5 times of variation does not have a significant impact on the drag coefficient (eqn. 11) in the particle motion. Therefore, the particle velocity in the cases of different LWC are almost the same (Fig. 3c). Similarly, the evaporation of a single particle has little effect on the particle temperature (Fig. 3b). The particle temperature of case 13 at engine inlet is only about 0.5K higher than that of case 11. Therefore, the change of LWC at spray mast has little effect on particle temperature and particle velocity.

{| style="width: 100%;border-collapse: collapse;" 
|-
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image46-c.png|294px]] 
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image47-c.png|288px]] 
|-
|  style="text-align: center;"|(a)
|  style="text-align: center;"|(b)
|-
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image48-c.png|282px]] 
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image49-c.png|306px]] 
|-
|  style="text-align: center;"|(c)
|  style="text-align: center;"|(d)
|-
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image50-c.png|300px]] 
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image51-c.png|288px]] 
|-
|  style="text-align: center;"|(e)
|  style="text-align: center;"|(f)
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 3. Comparison of results for different LWC (a. Air temperature; b. Particle temperature; c. Particle velocity; d. MVD; e. LWC; f. Relative humidity)</div>

==3.3 Effect of Relative Humidity==

Similarly, to understand the effect of initial RH on the variables at engine inlet, only change the RH value at spray mast and keep the remaining conditions unchanged. The specific cases of different RH are shown in Table 3.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Table 3. Cases of different RH</div>

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Case 21
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Reference Case
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Case 22
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Case 23
|-
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|LWC (g/m<sup>3</sup>)
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|0.99
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|0.99
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|0.99
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|0.99
|-
|  style="text-align: center;vertical-align: top;"|RH (%)
|  style="text-align: center;vertical-align: top;"|10
|  style="text-align: center;vertical-align: top;"|31
|  style="text-align: center;vertical-align: top;"|50
|  style="text-align: center;vertical-align: top;"|80
|-
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|MVD (μm)
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|20
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|20
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|20
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|20
|}


The comparison of results for different RH are shown in Figure 4, in which the comparison data are also at the central axis of the engine icing wind tunnel.

As the initial relative humidity at spray mast increases gradually on the premise of the same initial LWC and particle MVD, it is obvious that the evaporation phenomenon in the engine icing wind tunnel is gradually weakened. In general, the initial relative humidity is inversely proportional to the total evaporation and single particle evaporation. The relative humidity of case 23 is 80% at the spray mast, and the total evaporation loss is only 34.6% at the engine inlet, while case 21 is almost evaporated completely at the engine inlet (Fig. 4e). Similarly, comparison of particle MVD reflects the evaporation of single particle. The MVD of case 23 is only reduced by 13.5%, while the particles of case 21 is evaporated completely when they moved to half of the wind tunnel (Fig. 4d). This also explains the MVD of NRC particles increases during the test. Although all the particles have some degree of evaporation in the strong evaporation environment (31%), some of the smaller particles evaporate completely which improves the statistical value.

From the comparison of relative humidity changes, it can also be found that the curve of relative humidity changes more smoothly with the increase of initial value, which also reflects the influence of initial relative humidity on evaporation phenomenon (Fig. 4f). Similarly, the change of temperature of the background flow field (Fig. 4a) and the particle temperature (Fig. 4b) are inversely proportional to the initial relative humidity. It should be noted that the initial temperature of particles (293.15K in this paper) has little effect on the subsequent temperature changes (Fig. 3b and Fig. 4b). In addition, as the results of different LWC, the initial relative humidity has little effect on the change of particle velocity (Fig. 4c).

{| style="width: 100%;border-collapse: collapse;" 
|-
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image52-c.png|294px]] 
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image53-c.png|288px]] 
|-
|  style="text-align: center;"|(a)
|  style="text-align: center;"|(b)
|-
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image54-c.png|282px]] 
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image55-c.png|306px]] 
|-
|  style="text-align: center;"|(c)
|  style="text-align: center;"|(d)
|-
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image56-c.png|300px]] 
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image57-c.png|282px]] 
|-
|  style="text-align: center;"|(e)
|  style="text-align: center;"|(f)
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 4. Comparison of results for different RH (a. Air temperature; b. Particle temperature; c. Particle velocity; d. MVD; e. LWC; f. Relative humidity)</div>

==3.4 Effect of Particle MVD==

Finally, the effect of initial particle MVD on variables are obtained. Table 4 shows the cases of different MVD.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Table 4. Cases of different MVD</div>

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|Case 31
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|Reference Case
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|Case 32
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|Case 33
|-
|  style="border-top: 1pt solid black;text-align: center;"|LWC (g/m<sup>3</sup>)
|  style="border-top: 1pt solid black;text-align: center;"|0.99
|  style="border-top: 1pt solid black;text-align: center;"|0.99
|  style="border-top: 1pt solid black;text-align: center;"|0.99
|  style="border-top: 1pt solid black;text-align: center;"|0.99
|-
|  style="text-align: center;"|RH (%)
|  style="text-align: center;"|31
|  style="text-align: center;"|31
|  style="text-align: center;"|31
|  style="text-align: center;"|31
|-
|  style="border-bottom: 1pt solid black;text-align: center;"|MVD (μm)
|  style="border-bottom: 1pt solid black;text-align: center;"|10
|  style="border-bottom: 1pt solid black;text-align: center;"|20
|  style="border-bottom: 1pt solid black;text-align: center;"|50
|  style="border-bottom: 1pt solid black;text-align: center;"|100
|}


The comparison of results for different MVD are shown in Figure 5, in which the comparison data are also at the central axis of the engine icing wind tunnel.

With the same initial LWC, increasing the MVD of particles reduces the amount of particles and the total surface area of contact with air which reduces the total amount of evaporation.

Therefore, the larger the initial MVD of particles, the less LWC loss at the engine inlet (Fig. 5e), and the smaller the change of MVD due to the influence of evaporation on particle diameter (Fig. 5d). However, the LWC variation [javascript:; trend] of case 33 in Fig. 5e is particularly noteworthy. Which decreases slowly and then increases again. This is caused by the inertia of particles.

Figure 6 shows the LWC comparison of case 33 and reference case. the particles of case 33 are difficult to change the direction of motion in the straight section of the wind tunnel due to its large particle diameter and large inertia. Therefore, a shadow zone is formed on the wall of the straight section, and the particles concentrate towards the center of the flow field, resulting in an increase of LWC at the central axis.

In fact, case 33 is affected by evaporation, and the average LWC at engine inlet is only 82.5 g/m<sup>3</sup>. The larger the initial MVD of particles, the smaller the total evaporation and single particle evaporation, and the smaller the final relative humidity change. Relative humidity of case33 only increased by 4%, while case 31 increased by nearly 28 % (Fig. 5f), which also had less effect on the temperature of background flow field (Fig. 5a). In addition, different from the first two groups of comparison cases, the initial MVD of a particle has a more obvious change on the particle temperature. The larger the MVD of a single particle, the greater the internal energy of the particle, and the slower the particle temperature drops (Fig. 5b). The MVD of case 33 at spray mast is 10 times of that of case 31. Due to different evaporation degrees, the proportion of MVD at engine inlet increases to about 43 times, which is enough to influence the drag coefficient of particles. The larger the MVD, the slower the particle velocity. However, in the straight section of the wind tunnel, the particle velocity difference becomes smaller again with the influence of the particle inertia, (Fig. 5c).

{| style="width: 100%;border-collapse: collapse;" 
|-
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image58-c.png|288px]] 
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image59-c.png|288px]] 
|-
|  style="text-align: center;"|(a)
|  style="text-align: center;"|(b)
|-
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image60-c.png|282px]] 
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image61-c.png|300px]] 
|-
|  style="text-align: center;"|(c)
|  style="text-align: center;"|(d)
|-
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image62-c.png|300px]] 
|  style="text-align: center;width: 0%;"|[[Image:Draft_TAN_199218120-image63-c.png|288px]] 
|-
|  style="text-align: center;"|(e)
|  style="text-align: center;"|(f)
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 5. Comparison of results for different MVD(a. Air temperature; b. Particle temperature; c. Particle velocity; d. MVD; e. LWC; f. Relative humidity)</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_TAN_199218120-image64.png|600px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 6. Comparison of LWC (up: Case 33; down: Reference Case)</div>

==4. Conclusions ==

A three-dimensional particle trajectory calculation model based on Eulerian description was developed. The NRC small wind tunnel is taken as the research object and the calculation results are verified by the experimental results. In addition, the effects of LWC, relative humidity and MVD at spray mast on the variables at engine inlet were analyzed.

(1) The results show that increasing the initial LWC will greatly increase the total evaporation, but will reduce the evaporation of single particle, while simply increasing the initial relative humidity or initial MVD will reduce the total evaporation and single particle evaporation.

(2) The particle with larger MVD has higher resistance and is affected by inertia. The flow field in the straight section of the wind tunnel concentrates to the center of the wind tunnel and forms a shielding area on the wall of the straight section, which may cause data acquisition error.

(3) The temperature change of particles is greatly affected by relative humidity. The particles with larger initial MVD have a slower cooling rate because of higher initial internal energy. Otherwise, the change of particle temperature is not affected by the initial particle temperature.

'''Funding: '''This research was funded by the Open Fund of Key Laboratory of Icing and Anti/De-icing (Grant No. 1901IADL20190306).

'''Conflicts of Interest: '''The authors declare no conﬂicts of interest.

==References ==

:1. Yuan, Q. H.; Fan, J.; Bai, G. C. Review of Ice Crystal Icing in Aero-Engines. J. Propul. Technol. 39, 2641-2650, 2018.

:2. Bravin, M. and Strapp, J.W., “A Continuing Investigation of Diurnal and Location Trends in an Ice Crystal Icing Engine Event Database,” SAE Int. J. Advances & Curr. Prac. in Mobility 2(1):90-105, 2020

:3. Mason, J. G; Strapp, J. W; Chow, P. The Ice Particle Threat to Engines in Flight. AIAA rep. AIAA 2006-206.

:4. Pasztor, A. Airline Regulators Grapple with Engine-Shutdown Peril. Wall Street Journal. 251, 2008.

:5. Villedieu, P.; Trontin, P.; Chauvin, R. Glaciated and mixed-phase ice accretion modeling using ONERA 2D icing suite. AIAA rep. AIAA 2014-2199.

:6. Zhang, L. F.; Liu, Z. X.; Zhang, M. H. Numerical simulation of ice accretion under mixed-phase conditions. P I Mech. Eng. G. J. Aer. 230, 2473-2483, 2016.

:7. Norde, E.; Weide, E. T. A.;Hoeijmakers, H. W. M. Eulerian Method for Ice Crystal Icing. AIAA J. 1-13, 2017.

:8. Nilamdeen, S.; Habashi, W. G. FENSAP-ICE: Modeling of Water Droplets and Ice Crystals. AIAA Aerospace Sciences Meeting and Exhibit, 2009-4128.

:9. BU X Q, LI H, HUANG P, et al. Numerical simulation of mixed phase icing on two-dimensional airfoil. Acta Aeronautica et Astronautica Sinica, 2021, 42(3): 124085(in Chinese).

:10. Currie, T. C.; Struk, P. M.; Tsao, J. C.; Fuleki, D.; Knezevici, D. C. Fundamental Study of Mixed-Phase Icing with Application to Ice Crystal Accretion in Aircraft Jet Engines. AIAA Aerospace Sciences Meeting and Exhibit, 2012-3035.

:11. Fuleki, D. M.; Macleod, J. D. Ice Crystal Accretion Test Rig Development for a Compressor Transition Duct. AIAA Aerospace Sciences Meeting and Exhibit, 2010-7529.

:12. Mason, J. G.; Chow, P.; Dan, M. F. Understanding Ice Crystal Accretion and Shedding Phenomenon in Jet Engines Using a Rig Test. ASME rep. ASME GT 2010-22550.

:13. Struk, P. M.; Tsao, J. C.; Bartkus, T. P. Plans and Preliminary Results of Fundamental Studies of Ice Crystal Icing Physics in the NASA Propulsion Systems Laboratory. AIAA Aerospace Sciences Meeting and Exhibit, 2016-3738.

:14. Wright, W., Jorgenson, P., and Veres, J., Mixed Phase Modeling in GlennICE with Application to Engine Icing, in: AIAA Atmospheric and Space Environments Conference 2010, Canada, August 2-5, 2010.

:15. Hauk, T., Roisman, I., and Tropea, C., Investigation of the Melting Behavior of Ice Particles in an Acoustic Levitator, in: 11th AIAA/ASME Joint thermophysics and heat transfer conference, USA, June 16-20, 2014.

:16. Davison, Craig R., James D. MaCleod, and Jennifer L. Chalmers. Droplet Evaporation Model for Determining Liquid Water Content in Engine Icing Tunnels and Examination of the Factors Affecting Liquid Water Content, in: 9th AIAA Atmospheric and Space Environments Conference, 2017.

:17. Bartkus, T. P., Struk, P., & Tsao, J. C.  Development of a coupled air and particle thermal model for engine icing test facilities. SAE International Journal of Aerospace, 8(2015-01-2155), 15-32.

:18. Spalart, P. R.; Allmaras, S. R. A one-equation turbulence model for aerodynamic Flows. AIAA Aerospace Sciences Meeting and Exhibit, 1992-0439.

:19. Crowe, C., Sommerfeld, M., and Tsuji, Y., Multiphase Flows with Droplets and Particles. CRC Press, 101-103, 1998.

:20. Zhang, Y., Ozcer, I., Nilamdeen, S., Baruzzi, G.S. et al., Numerical Demonstration of the Humidity Effect in Engine Icing. SAE Technical Paper, 2019-01-2015, 2019.