CFD Model for Aircraft Ground Deicing: Veriﬁcation and Validation of an Extended Enthalpy-Porosity Technique in Particulate Two Phase Flows

: Researchers have focused in the last ﬁve years on modelling the aircraft ground deicing process using CFD (computational ﬂuid dynamics) in order to reduce its costs and pollution. As preliminary efforts, those studies did not model the ice melting nor the diffusion between deicing ﬂuids and water resulting from the melting process. This paper proposes a CFD method to simulate this process ﬁlling these gaps. A particulate two-phase ﬂow approach is used to model the spray impact on ice near the contaminated surface. Ice melting is modelled using an extended version of the enthalpy-porosity technique. The water resulting from the melting process is diffused into the deicing ﬂuid forming a single-phase ﬁlm. This paper presents a new model of the process. The model is veriﬁed and validated through three steps. (i) veriﬁcation of the species transport. (ii) validation of the transient temperature ﬁeld of a mixture. (iii) validation of the convective heat transfer of an impinging spray. The permeability coefﬁcient of the enthalpy-porosity technique is then calibrated. The proposed model proved to be a suitable candidate for a parametric study of the aircraft ground deicing process. On the validation test cases, the precision of heat transfer prediction exceeds 88%. The model has the ability of predicting the deicing time and the deicing ﬂuid quantities needed to decontaminate a surface.


Introduction
Several researches have focused on the development of CFD models for deicing aircraft in flight to help in the design of ice protection system [1][2][3]. On the ground, aircrafts also undergo icing phenomenon. This phenomenon caused several accidents [4,5], because of which the aircraft ground deicing was designed and imposed. The aircraft ground deicing process used at most airports is carried out using an impinging spray. Process guidelines are suggested in [6][7][8][9] based on technical and manufacturer reports, guidance materials and standards and technology patents and procedures. Considering the criticality of the process, the safety is prioritized over economic aspects. Therefore, the guidelines focus more on the results than on how to achieve them. SAE AS 6285 [10] provides industry standards for the methods and procedures used in performing the treatments necessary for the proper deicing and anti-icing of aircraft on the ground. AS6285 forms with AS6286 and AS6332 the "global aircraft deicing standards". The revised version AS6285D includes references to the Federal Aviation Administration (FAA) and Transport Canada (TC) deicing program guidance. The suggested procedures give general guidance such as to start the deicing at the top for the vertical surfaces or the direction of the spray shall be from the wings leading edge to the trailing edge. The inclination angle between the spray direction and the surface normal is not mentioned. However, numerous experimental studies [11][12][13][14] have shown that variation in spray parameters, such as the inclination angle and the distance,

•
Convection fusion (natural or forced): It occurs if a solid is placed in a fluid domain with certain pressure and temperature fields. In natural convection, no flow is imposed, the pressure or/and temperature gradients generated at the solid-fluid interface induces the movement. This movement prevents the solid from having a cold layer of fluid around it. This movement persists until a point of equilibrium is reached by melting the solid [17]. In forced convection, the flow lowers the pressure and/or increases the temperature at the solid-fluid interface, causing it to melt [18]. • Close Contact Fusion: It occurs if a heat source and a solid are brought into contact with each other during the solid fusion. The physical situation involves the movement of the heat source or the solid, which prevents the accumulation of the melt between the source and the solid [19].
The ice melting in AGD process occurs because of the spreading of an ADFs film. The film is induced by the spray flow; therefore, this process involves forced convection. From a numerical point of view, numerical methods modelling convection fusion can be divided into two families: (i) The multi-region methods [20] consists in defining two regions in which two systems of different equations are solved. A dynamic mesh is used to track the melting front. (ii) The single region methods [21] consist of using a single region and defining a scalar field, modelling the liquid/solid volume fraction in each cell. Source terms are added to the conservation equations to model the ice presence effect on the flow. The liquid/solid volume fraction is updated at each time-step, according to the energy transfer and a threshold temperature. In this method, a single region is sufficient, hence a static mesh is used, which makes it more interesting in the AGD process.
Indeed, the AGD process is characterized by a large dispersion of the geometric length scales. [16] opted for a particulate two-phase flow to lower the computations needs in term of mesh refinement. They validated the convective heat transfer of an impinging spray on a flat plate. This paper suggests coupling a solidification/melting model to a particulate two phase flow. The single region method is selected for this task due to its simplicity of implementation. It consists of adding (i) a solid fraction field that is updated at each iteration based on the temperature field and (ii) the adequate source terms to the governing equations. The enthalpy-porosity technique [21] is used in this paper in a Eulerian particulate two-phase flow based on [22] already implemented in OpenFOAM, more precisely in its native solver twoPhaseEulerFoam. The Eulerian approach is favored over the Lagrangian for the fact that the liquid/solid phase fraction used by the melting method is a Eulerian field. Using the Lagrangian method will lead again to the use of a dynamic mesh. Instead, a multi-region approach would be adopted. Figure 1 presents a typical computational domain of the problem. A first region is defined for the liquid spray, where the multiphase flow is treated with a Eulerian-Lagrangian model. The second region relatively thin, is defined for the spray impingement on the ice. The multiphase flow in this region is carried out with a Eulerian-Eulerian model based on [22] works coupled to the enthalpy-porosity technique to model ice melting. The Lagrangian-Eulerian region is fully described in [16]. The present paper focuses on the verification and validation of the Eulerian-Eulerian model where an extended version of the enthalpy-porosity technique is proposed and implemented. The original version of the enthalpy-porosity technique is designed for a liquid-solid system. It is widely used in internal flow simulations, particularly in the context of Phase-Change Material (PCM) for heat storage. In this paper, it is extended to a gas-liquid-solid system. It is also calibrated for a computational domain characteristic of the AGD process. Previous works [15,16] presented an impinging spray model and analyze the convective heat transfer on a clean surface. This paper is a complementary effort to simulate the ice decontamination in the context of a numerical test bench. To the knowledge of the authors, this is the first effort in literature that implements the enthalpy-porosity technique to a gas-liquid-solid system. The extended enthalpy-porosity technique is a new candidate for simulating a variety of problems involving solidification and melting phenomena in external flow.
The following section describes the mathematical model and the numerical method of the Eulerian-Eulerian region. The third section presents a sequence of three test cases for the V&V assessment which followed by a calibration of the solidification/melting model.

Conceptual Model
The problem could be seen as the interaction of two phases (gas and liquid). The solid phase is implicitly treated as a liquid phase with a velocity field reduced to zero. The Eulerian-Eulerian region of Figure 1 presents a typical computational domain of the problem. Initially, the water is at solid state at an ambient temperature below the melting point. The hotly injected ADF heats the ice and melts it. The liquid water coming out from the melting process dissolves in ADF, forming a mixture. First, the conservation equations are presented. It is followed by the numerical method describing the solver algorithm and the numerical schemes.

Mathematical Model
The model consists of continuity, species, momentum and energy equations. Those equations are designed for a gaseous phase and a liquid phase. To model ADF and water mixing, the liquid phase is considered as a mixture of two liquid species. A specie transport equation is solved just after the continuity equation. The enthalpy-porosity technique is used to model the ice presence in the momentum and energy equations. In the following equations, indexes "g" and "l" refer, respectively, to the gas and liquid phases; "s1" and "s2" refer, respectively, to the liquid 1 and liquid 2 species. Below, liquid 1 is set as ADF and liquid 2 is set as water.

Mass Conservation
The mass conservation (continuity) equation is analogue to the VOF method [22]: with α i , ρ i , u i are, respectively, the volume fraction, density and velocity of entity i. Then, the specie transport [23] equation is solved to update concentrations and the liquid thermophysical properties.
with D l is the mutual diffusion coefficients of specie 1 and specie 2 defined as: D l = ν l /Sc l with ν l and Sc l the kinematic viscosity and the Schmidt number of the liquid phase. Right after solving Equation (2) the liquid phase thermophysical properties are updated using the mixture rules presented in Table 1.

Thermophysical Properties Expression
Density

. Momentum Conservation
After resolving the mass conservation equations, the phases are considered (in each cell) as continuous and dispersed following two threshold values: (i) A F the minimum volume fraction of a phase to be considered fully continuous and by default the other phase will become a dispersed phase. The dispersed phase is modelled as a set of droplets characterized by a diameter and a shape. For sake of simplicity, the shape of the droplets is considered as a perfect sphere. (ii) A P the minimum volume fraction of a phase, for which it can be treated as partly continuous phase, i.e., the cell is part of an interface (A F > A p ).
Two momentum equations of the same form are solved for each phase (gaseous and liquid) [22]: where µ I the viscosity of phase i and the rate of strain (deformation) tensor D(u i ) is defined as D(u i ) = 1 2 ∇u i + (∇u i ) T . S sms,u i is an added source term modelling the melting/solidification of phase i, it is discussed in Section 2.1.4 S u i is the inter-phase momentum transfer with the convention of S u g + S u l = 0 . The two phases exchange momentum in form of drag and virtual mass S u,i = S drag,i + S vm,i . The drag and virtual mass forces are due, respectively, to the velocity difference and the acceleration difference between the two phases. In the following equations, quantities difference presented by the symbol ∆ are defined as: ∆χ = χ other phase − χ phasei .
For cells with dispersed phase [24]: For interface [25]: with C drag and K drag the drag coefficients and A disp the dispersed phase surface. The drag coefficient models implemented in OpenFOAM are presented on Table A1 in the Appendix A.
Since droplets are supposed to have a spherical shape, the virtual mass coefficient is fixed to 0.5 [24].

Energy Conservation
Two energy equations are solved for gaseous and liquid phase [26], as shown by Equation 6.
with h, k and α e f f , respectively, the enthalpy, the kinetic energy and the effective thermal diffusivity. S sms,h is a sink term modelling phase change. K ht is the heat transfer coefficient between the two phases. Two models are already implemented in OpenFOAM. The first is an analytical model for a perfect sphere and the second is a correlation for turbulent heat transfer also for a sphere. The heat transfer coefficients are presented in Table A2 in the Appendix A.

Solidification/Melting Source Terms
The source terms S sms,u i and S sms,h i the solid phase presence effect. They are calculated following the enthalpy-porosity technique. In this technique, the solid-liquid free boundary is not tracked explicitly. Instead, a phase fraction indicator field δ (0 for solid and 1 for liquid) is updated at each iteration based on an enthalpy balance. Originally, the enthalpyporosity technique was developed and used for single phase flows, in this paper a modified version of the technique is proposed.
The momentum sink S sms,u i , introduced in Equation 3, classically models the buoyancy effect S buoy due to the thermal expansion of ice and a drag force S sd exerted by ice on the liquid. However, for the gas phase the buoyancy term should be removed since the implicit solid phase belongs to the liquid phase. Additionally, in the proposed problem, the buoyancy effect for the liquid phase can be neglected since the flow is under forced convection (impinging spray).
In the original version of the enthalpy-porosity technique [21], the solid phase is not designed to move. The present work focusses on the decontamination of a surface, and since no ice breakup is modelled, this characteristic will be preserved.
The Darcy law for flow through porous media is used. For isothermal solidification/melting, permeability has no physical significance; however, it is used classically as a numerical technique to estimate the velocity at the mushy region [21].
where q and Cu are model coefficients, which are discussed in the numerical method paragraph.
where h l = h l,s + h l,L , with h l,s = c p,l ·T l is the sensible enthalpy and h l,L = δL is the enthalpy of fusion (or latent heat of fusion).
Moving the D Dt (α l ρ l h l,L ) on the right side of the energy equation and developing it, gives: which by identification with Equation (6), gives: where S sms,h l is the energy sink term due to phase change of the liquid phase (solid/liquid). It represents also the latent heat released during solidification. The gaseous phase does not undergo a phase change (the deposition/sublimation phenomenon is not modelled), thus the term S sms,h g is null.

Numerical Methods
This section describes the numerical methods used in all the computations presented in this paper. It is divided into numerical treatment and solution control. The numerical treatment of the particulate two-phase flow solver is fully described in [26]. The next paragraph discusses the numerical treatment of the solidification/melting source. It is followed by the numerical schemes and the algorithm controls details for reproducibility purpose.

Numerical Treatment
Concerning the Darcy law coefficient presented at the Equation (9), q is a small number to prevent division by zero (default 10 −3 ). Cu is the mushy region sink coefficient which measures the amplitude of the damping. It is also known in the literature as the permeability coefficient [27]. The higher this value, the steeper the transition of the velocity of the material to zero as it solidifies. Very large values may cause the solution to oscillate (OpenFOAM default value 10 5 ). In the liquid region (δ = 1), S sd is equal to zero and does not affect the momentum equation. In the solid region, (δ = 0) S sd is equal to −Cu/q which is high enough to dominate all other terms in the momentum equation. Therefore, the momentum equation will be reduced to S d ·u = 0, and, consequently, u = 0.
Ebrahimi, Kleijn and Richardson [27] investigated the sensitivity of numerical predictions to the permeability coefficient Cu. The study revealed that reducing the cell size at the mushy zone (free boundary at the ice vicinity) diminishes the influence of Cu on the results. According to [27], values between 10 4 and 10 8 are often applied in the modelling of energy storage systems. The Cu value is often tuned for every set of boundary conditions and material properties.
In the case of AGD process, there is a remarkable lack of experiment data. Tuning the permeability coefficient to match numerical results with experimental ones is impossible for the moment. Therefore, the default value of 10 5 is used at first. Once the computation parameters are fixed after a grid convergence study, the sensibility to the permeability coefficient will be studied. If dependencies are revealed, the mesh should be refined. At that state, the permeability coefficient could be tuned for the first mesh selected through the mesh convergence study.
Finally, the solid phase fraction is updated at each time step i in every cell n in two steps [21]: with T sol and Cp l the temperature of the phase change and the heat capacity of the liquid phase and r an under-relaxation factor.

Numerical Schemes and Algorithm Controls
Conservation equations are solved using the PIMPLE algorithm of OpenFOAM [28]; Within one time-step, the algorithm search a steady-state solution with under-relaxation and then go in time. Three outer correction loops are used to ensure that the explicit parts of the equations are converged. After reaching a defined tolerance criterion of 1e −6 within the steady-state calculation, the algorithm leaves the outer correction loop and move on in time. Conservation equations are solved using a Gauss Siedel smooth solver, and the pressure is then solved using the geometric algebraic multi-grid (GAMG) solver [29] with a diagonal incomplete-Cholesky (DIC) smoother/preconditioner.
OpenFOAM's standard numerical schemes are used to solve the mathematical model. Time derivative terms are computed using an implicit, bounded first-order transient discretization scheme [30]. Gradient terms are computed using a face-based third-order Gauss scheme, which is based on the Gaussian divergence theorem for the volume integral of a gradient field [31]. The third order is used to improve the spatial accuracy by reducing the numerical diffusion especially for volume fractions α i , Y i and δ. Second-order tensors divergence terms such as divergence terms of the momentum equations are computed using a limited second-order Gauss scheme with a single limiter. The limiter is calculated based on the direction of most rapidly changing gradient. Divergence terms including energy, pressure and specie terms are computed using a limited second-order scheme that limits towards upwind in regions of rapidly changing gradient. Divergence terms including only phase fractions are discretized using the Van Leer divergence scheme. Finally, Laplacian terms are discretized using an uncorrected second-order Gauss scheme.

Verification, Validation and Calibration
In this section, three test cases are set for the V&V assessment of the mathematical model and the numerical method. The first test case is a simple diffusion test between two liquids. The results are compared to the ones computed with the native solver laplacianFoam. The second test case is for the validation of the computed temperature of a mixture. The results are compared to the experimental results published in [32]. The third test case is a validation of the convective heat transfer of an impinging hot spray on a flat plate [33]. At the end, the permeability coefficient is calibrated for an AGD test case.

Diffusion Verification Test
The goal of this test is to verify the code consistency. This is done by reducing the problem to a simple diffusion case. The results can then be validated against results from the laplacianFoam solver. laplacianFoam resolves a simple heat equation and was validated in several works [34][35][36].
The problem is reduced to a simple diffusion case by setting adequate initial and boundary conditions. Initially, the computational domain is filled completely by a liquid phase with uniform temperature and pressure, and a zero velocity. The thermophysical properties of the two species composing the liquid are identical to avoid buoyancy effect due to density gradients or variation of the diffusion coefficient due to viscosity gradients. The solidification melting source is deactivated for this test. The control volume, presented in Figure 2, is a 2D square of 1 m edge with wall boundaries on its four sides. Initially, the species compose a disc shape at the center of the control volume with an inner diameter of 0.3 m and an outer diameter of 0.6 m. This form is characteristic of what could be resulting from an impinging spray. Finally, the boundary field value for the ADF specie is set as a uniform constant equal to 1 to let the water species dissipate across the walls.
For the laplacianFoam case setup, the same mesh is used with similar initial and boundary conditions for the temperature field (conform to the water specie). All computations are done with a mutual diffusion coefficient D l = 2·10 −3 m 2 ·s −1 giving a Schmidt number Sc l = 5·10 −4 for a kinematic viscosity of ν l = 1·10 −6 m 2 ·s −1 . This setup permits to have a maximum volume fraction inferior to 0.1 after 50 s of time. It allows to evaluate the code precision for high and low volume fractions of the species. Twelve tests were performed involving three uniform meshes with different degree of refinement (50 × 50, 100 × 100 and 200 × 200 cells) and four different time steps (10 0 , 10 −1 , 10 −2 and 10 −3 s) . The computational domain is discretized into uniform square cells with constant edge length in both X and Y directions. For better understanding of the test case, 5 spatial distributions of Y water along the line probe are presented in Figure 3 with a dimensionless distance X = x/ √ 2m. The finest mesh and the smaller time step were used in this figure. Computations with the other meshes are used to demonstrate the grid convergence. The distribution keeps its symmetry. The integral encompassed by one curve and the x axis decreases over time. At the final time (50 s), the field's gradient ∂Y water /∂X is attenuated, meaning the mixture is more homogeneous. A maximum value less that 0.1 is present at the center of the computational domain. The verifications are made in two steps. The first step is to generate the reference data. It consists of a time-convergence study followed by a grid-convergence study of the laplacianFoam case. For each mesh, the time step induced error is evaluated for each 1 s of time following Equation (14).
where dt designates the time step and t the time. The time step induced error distribution function through time is presented on the Figure 4 for the three meshes. For all the computations, the error has the same time evolution. It decreases for the first 10 s then stabilize at a constant value. The time convergences are also quite similar for the three meshes. The constant value of the error decreases from 0.18% for a time step of 1s to 0.002% for a time step of 0.01 s. For the grid convergence study, the results for the smallest time steps are used. The grid induced error is evaluated following Equations (15) and (16). Figure 5 presents the grid induced error evolution through time. The two errors are similar during the first 8 s. They start at 2% and attain 0.75% after 8 s of time. After 8 s, the relative error between the Coarse mesh and the Medium mesh keeps decreasing progressively until the end of the simulation. The relative error between the Medium mesh and the Fine mesh decreases faster to reaches 0.1% at 13 s and stay constant until the end of the simulation. The mean values of the relative errors are 0.61% and 0.34%, respectively, for Coarse-Medium and Medium-Fine meshes.
The same calculations are done with the AGDEulerFoam. The results precision is assessed using a relative and an absolute error presented, respectively, by Equations (17) and (18).
with Y A and Y laplacian , the results of the AGDEulerFoam and the laplacianFoam solvers. Y abs stands for the results of the laplacianFoam solver with the finest mesh and smallest time step. The maximum and mean values of the errors are presented in Table 2.

Time
Step For all tests, the mean relative error is of the order of 10 −4 %, which shows that the configured solver gives practically the same results as the laplacianFoam solver. The relative error is independent from the time step. Refining the mesh reduces the relative error as seen on the maximum value of the relative error. The analysis of the absolute error leads to the same conclusion as the time step and grid convergence studies of the laplacianFoam solver.
This test shows the validity of the diffusion module for the solver responsible for the miscibility water in the ADF without heat exchange. It also reassures about the code consistency. Indeed, by defining a computational domain containing a single phase, no additional phase was created through the iterations.

Miscibility -Energy Validation Test
The second test case is for the validation of the computed temperature of a mixture. Results are compared to the experimental results published in [32]. The test case is a dambreak designed for the mixture of two species at different temperature. The temperature is monitored during 44.5 s with 10 thermistors placed at unique positions. Figure 6 presents the computational domain at the initial time and thermistors positions. The two species are presented by "L 1 " and "L 2 ", thermistors positions are noted by "S 01 S 02 . . . S 10 ". Initially L 1 is placed in the compartment C 1 at 291.95 K and L 2 is placed in the compartment C 2 at 325.35 K. Four different uniform meshes, which characteristics are presented in Table 3, with four different CFL numbers, were used for a grid and time-step convergence study. The results were treated with the same methodology used in [32]. The probes results were compared to the thermistors results presented in [32]. Table 4 present the maximum and means errors for the different computations. The presented error is the average error over the 10 probes. The best result is reported for the Fine mesh with the larger CFL number (CFL = 1). However, for the same CFL number, the Medium and Extra-Fine meshes presents quite comparable errors. It is also to be mentioned that for the Fine mesh, the CFL dependency of the error is low. In the current test case, the temperature is monitored on discrete points and the dam break creates eddies that may or may not pass by the probe locations. Initial perturbations because of numerical errors may also affect the eddies positions and sizes which reflects on results. Figure 7 shows the eddy presence in the flow. It presents the liquid temperature field at 0.5 s, 4 s and 10 s for the four meshes and a CFL = 1. The black dots reflect the probe locations. This figure reveals the formation of two eddies; (i) A hot eddy is created at 0.5 s caused by the propagation of L 2 in the compartment C 1 and (ii) A cold eddy is present at 4 s caused by the propagation of L 1 in the compartment C 2 . For the Fine and extra Fine meshes, the cold eddy is divided in two eddies. For comparison, Figure 8 present the best result obtained with the AGDEulerFoam solver versus the result of I MTFoam validated in [32]. Initially, the two solvers present the same accuracy level. After approximately 7 s, the AGDEulerFoam shows better results. The error does not exceed 5%. The maximum error for I MTFoam is of the order of 7%. Overall, the I MTFoam solver presents an average error of 4.85%, greater than the 3.30% error of the AGDEulerFoam. Again, the space-discrete data monitoring of an eddy characterized phenomenon prevents conclusions about which solver is more accurate.
However, the AGDEulerFoam provides results comparable to the ones from the validated I MTFoam. Therefore, results are considered satisfying. It is also to be mentioned that [32] is, to the best of the authors' knowledge, the only study presenting experimental results about energy transfer between miscible liquids in a two-phase flow.

Impinging Jet Validation
The AGD is carried out by a hot dispersed spray. For validation, the solver is tested on a less dispersed non-isothermal free-surface flow of a liquid jet impinging on a heated surface [33]. This test case serves to assess the accuracy of the predicted pressure distribution exerted by the jet on a flat plate. It also assesses the accuracy of the convective heat transfer between the liquid film resulting from the impingement on the heated surface. The pressure distribution is compared to Tong [37] numerical results. The convective heat transfer is compared to Edin and Šefko numerical results [38] and Stevens experimental results [33].
The computational domain of the test case is presented in Figure 9. It consists of an axis-symmetric (2D) domain. Initially, the domain contains only air at ambient temperature (20 • C). The water is injected at the same temperature from the inlet with a velocity U 0 . A constant heat flux q wall is applied at the plate. The liquid jet impinges the surface, and the liquid is heated due to the constant heat flux applied to the wall. The results are monitored once the steady state is reached. The relevant test case parameters are presented in Table 5.  The Reynolds number is defined as Re = 4Q/πdν with Q the inlet volume flow rate. The turbulent flow effects are neglected at the nozzle exit since the radial velocity of the thin liquid film spreading over the heated surface is relatively small. Therefore, it is assumed, as in [38], that the flow relaminarizes and stays laminar after the impingement. A uniform velocity is set at the inlet with a slip condition at the nozzle wall to get a uniform velocity at the nozzle-tip (at 3.7d from the plate). The domain is discretized using uniform square cells with constant edge length in both X and Y directions. Four cell sizes, presented in Table 6, were tested for a grid dependency study. The four meshes predicted the same pressure distribution (less than 1% of maximal variation). Thus, the effect of the cell size is presented, in this paper, only for the Nusselt number defined in Equation 19 [38]. Figure 10 presents the Nusselt number distributions and variations for the four tested meshes. The horizontal axis presents the dimensionless distance (r/D) to the stagnation point. The four profiles have the same tendency, the heat transfer increases slightly outwards the stagnation point (r/D < 0.7) then decreases hyperbolically (r/D > 0.7). Refining the mesh increases globally the convective heat transfer. The variation between the meshes increases outwards the stagnation points presenting a local extremum at r/D = 1.4 for the Coarse and Medium meshes. The variation between the Fine and Extra-Fine meshes increases slightly outwards from the stagnation point and remains under 3% and the local extremum is not present. Therefore, a cell size of D/40 suffices to reach a grid-independency for the prediction of the convective heat flux.  In what follows the Extra-Fine mesh results are compared to the literature results. Figure 11 presents a comparison between the pressure distribution predicted by the present simulation versus Tong [37] numerical results. As seen, the pressure distributions show a very good agreement. A maximum difference of 2.28% is monitored at 0.54D from the stagnation point. The mean difference between the two data set is 0.77%. Finally, Figure 12 compares the Nusselt number distribution of the present simulation versus Edin and Šefko numerical results [38] and Stevens experimental results [33]. Edin and Šefko OpenFOAM numerical results agreed well with the experimental results with a maximal error of 20% at approximately r/d = 0.6. The present simulation presents the same error profile with a maximal error of 12%. The solver is considered valid for the heat transfer prediction of an impinging jet. In the next paragraph, the same computational domain is used for a deicing simulation.

Permeability Coefficient Calibration
In this paragraph, an ice block is introduced to the computational domain of the impinging jet test case. Figure 13 presents the main modifications made to the initial conditions. The meshes are the ones used for impinging test case. The liquid is injected with the same Reynolds number as in Section 3.3 at 60 • C and it is composed by ADF instead of water. The internal field contains air at −5 • C. An ice block (9D × D/2) composed by water is set on the bottom-wall boundary at −5 • C. Finally, the constant heat flux boundary condition on the bottom-wall is replaced by a zero-gradient boundary condition.
The thermophysical properties of the 3 species (air, water and ADF) are presented in Table 7. At each iteration, the density field is updated following the phase equation of state. The gaseous phase is treated as an ideal gas and the liquid phase is treated as a perfect fluid. ADFs are mainly composed of ethylene glycol. They also contain corrosion inhibitor and other fluids to reduce their surface tension and increase their viscosity. Their composition varies also depending on the manufacturer. As the aim of this paragraph is to calibrate the permeability coefficient, then the ADF will be considered as pure ethylene glycol.  First, the Fine mesh is used to run a first set of computations with four permeability coefficient values (10 4 , 10 5 , 10 6 and 10 7 ). Then, if a sensibility is detected, the mesh will be refined until the results become unsensitive to the permeability coefficients. Figure 14 presents the ice volume evolution through time for the different permeability coefficients using the Fine mesh. The permeability coefficients 10 5 , 10 6 and 10 7 gave approximately the same results. The permeability coefficient of 10 4 provides, however, a different melting rate. A maximal difference of 13% was monitored between this computation and other results. Figure 15 presents the same results using the Extra-Fine mesh. The same inference could be made for this mesh as for the Fine mesh with the exception that maximal difference was reduced to 4%. It can be concluded that the unsensitivity to the permeability coefficient is reached using the Extra-Fine mesh.
With the Fine mesh, the results are unsensitive to permeability coefficient between 10 5 and 10 7 . It is interesting for an engineering purpose to re-run the test-case using a Coarse mesh with a permeability coefficient of 10 6 . Figure 16 presents on the same graph the melting rates for the three tested meshes using a permeability coefficient of 10 6 .   The Coarse mesh provides a result similar to the Extra-Fine mesh with a maximal error of 4.19%. It is an interesting result giving that the Extra-Fine mesh contains 64 times cells number than the Coarse mesh. Ice melts slower using the Coarse mesh, which is consistent with the fact that the Nusselt number was also under-predicted for the impinging test-case using the Coarse mesh.

Discussion and Conclusions
A new OpenFOAM-V6 based solver (AGDEulerFoam) was developed for simulating a dispersed two-phase flow enabling the definition of the phases as a mixture of multiple species with different thermophysical properties. The enthalpy-porosity technique is extended to handle phase change in multiphase flows. It is coupled to the new solver to model the solidification (and melting) of a liquid phase (of a solid phase). The AGDEulerFoam was developed by programming a species transport equation in the pre-existing OpenFOAM twoPhaseEulerFoam native solver. The thermophysical properties of the multi-species phase are updated following the rule of mixtures at each iteration. The species miscibility module precision is assessed with a comparison versus the laplacianFoam solver validated in [34][35][36]. A maximum difference of 0.01% was recorded between the two solvers. Developments achieved in this study are part of a design methodology for a numerical bench test of the aircraft ground deicing process in order to improve its efficiency.
It is to be noted that no CFD method is defined in the literature to simulate the ice shape evolution in the aircraft ground deicing process. This paper presents a model first of its kind developed on an open-source CFD framework for the aircraft ground deicing simulation. This paper is a complementary effort to simulate the ice decontamination in the context of a numerical test bench. The main findings of the present paper are summarized into three points: 1.
The developed solver predicts unsteady temperature evolution of a two-phase flow in which one phase is a mixture of two species with a precision of 95%. This precision is observed in a 2D dam-break simulation involving mixture.

2.
The proposed solver predicts the convective heat transfer between liquid formed by an impinging jet and a heated wall with a maximal error of 12%. This result is comparable to the existing models.

3.
The permeability coefficient of the extended enthalpy-porosity technique is calibrated through a sensitivity study proposed in [24]. The results are insensitive to the permeability coefficient with an Extra-Fine mesh. The ice shape evolution can be well predicted with a Coarse mesh with a permeability coefficient of 10 6 , which is within the interval stated by Ebrahimi, Kleijn and Richardson [24].
In other words, the proposed CFD method is able to simulate a dispersed two-phase flow with the solidification/melting of one phase. Those characteristics make it able to simulate an AGD process. For better computations efficiency, the present model will be coupled to a Lagrangian model. The Lagrangian model will be used for simulating the majority of the spray far from the contaminated surface. The present model is used for the simulation of a relatively thin layer of the computational domain where the ice melting is modelled. The final solver will be used to simulate a 3D AGD test case and a parametric study of the process to determine the AGD process sensitivity to the ADF spray parameters.  Acknowledgments: CFD computations were made on supercomputers managed by Calcul Québec and Compute Canada.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.