Assessment of liquid film evaporation modelling in a turbulent channel flow

In Gasoline Direct Injection (GDI) engines, the interaction between the premixed flame front and the liquid wall film is considered as the primary source of soot. Modeling soot formation in this configuration requires the local gas mixture to be accurately predicted above the wall film which, in turn, requires an accurate wall film evaporation model. Models derived from semianalytical solution suitable for wall-modeled LES have been proposed in the literature, but their validation stay limited. This work aims to investigate the dependence of the film evaporation model to wall modeling and mesh refinement at the wall in LES. A turbulent channel flow is considered for this purpose. Wall-refined LES cases, using an evaporation model based on the local species gradient, serve as reference results. For the center of the wall cell located in the viscous sub-layer of the boundary layer (respectively within the buffer layer), the evaporation is found underestimated by about 45% (respectively by about 60%). According to these results, wall heat-flux models are not the main source of errors of film evaporation and the usage of classic wall laws seem not suited for wall-modeled LES of evaporating wall film.


Introduction
In Gasoline Direct Injection (GDI) engines, fuel is directly injected in the combustion chamber, leading to the formation of a thin wall film on the piston and cylinder walls. Depending on the engine condition, the formed film may not evaporate completely till the ignition phase. As observed experimentally [1], flame film interaction is considered as the primary source of soot formation in GDI engines. High fidelity simulation tools are essential to improve the understanding of soot formation process in this configuration, and thus, reduce particulate matter emissions. The film evaporation modeling is a key step in this process since it is the source of the locally rich regions where soot are formed. Consequently, this model accuracy is of utmost importance for further investigating soot emissions. The present study focuses on the film evaporation model itself. Thus, wall roughness and coke coating effects will be neglected. Moreover, wall temperatures will be constant. Cold start conditions are targeted because they are known to contribute largely to GDI engines soot emissions [1]. Accordingly, the imposed temperatures will be below the saturation temperature of the fuel so that the film evaporates in the complete wetting regime, i.e. without boiling. Finally, the wall-film is modeled using a Lagrangian approach as in actual engine simulations. Several models, with different physical and numerical complexity, can be found in the literature to express the evaporation rate. In the following and throughout the paper, all quantities with a ( + ) superscript, a ( s ) subscript, or a ( 1 ) subscript correspond respectively to wall units, wall or film surface values, or the values at the center of the wall mesh cell. The y direction corresponds to the wall normal direction. For a single-component fuel (denoted below with subscript F ), the evaporation rate can be derived from the mass balance equation at the liquid/gas interface as follows: where ρ s , D F,s and Y F,s are respectively the density of the gaseous phase, the fuel diffusion coefficient and the fuel mass fraction at the interface. Y F,s can be computed using the fuel partial vapor pressure P v,F given by the Raoult law: where W a is the molar weight of the mixture without fuel. Finally, the fuel mass fraction gradient reads: where ∆y is wall mesh cell size in the wall normal direction. The use of the expression (1) to express the evaporation rate requires a fine discretization of the gas phase near the wall film to well capture the species gradient. This approach will be referred here after as local-gradient approach. A less computationally demanding approach consists in modeling the evaporation rate as follows [2]: where c p is the mass specific heat capacity at constant pressure, Le is the molecular Lewis number, W is the molar weight of the gas mixture, B M is the Spalding mass transfer number defined as The exponents over the correction terms are set to n = 1/3 and m = 0 [3]. The heat transfer coefficient h is given by where φ s is the heat flux across the surface and T is the temperature. This approach directly depends on the wall heat flux model through the term φ s . It is referred hereafter as far-field approach. Although this approach is widely used for film evaporation in practical cases, its validation remains limited. In particular, to the authors' knowledge, its dependence on mesh refinement and heat transfer model are not yet investigated in the literature. Such validations are even more needed in the context of industrial simulations in which multiple models are used together and some, such as combustion models during flame-wall interaction may require very small cells and specific heat-transfer wall laws. Moreover, mesh size may also vary locally during the simulation with Adaptive Mesh Refinement (AMR). The stationary turbulent channel flow is an academic case suited for the study of near wall flows. It was used by Desoutter et al. [4] to conduct a Direct Numerical Simulation (DNS) of turbulent boundary layer over an evaporating liquid film. A similar configuration is used in the present work in order to investigate the dependence of the film evaporation model to wall modeling in Large Eddy Simulation (LES). A version of the CONVERGE [5] solver has been specifically developed for this work. Its performance on such academic cases is evaluated in the following against literature references starting with an isothermal case without film [12]. Given the used numerical methods, it is not expected that the results of this step could fully match the reference results. The case being also simulated with wall model allowing to choose which momentum wall law to use for the rest of the study. Then, the accuracy of the local-gradient approach and its validity as a reference for comparison against wall-modeled results are assessed by comparison against Desoutter et al. [4] results.  Finally, the channel flow cases specific to this study will be presented in order to evaluate the farfield approach results with four heat-transfer wall laws against local-gradient approach results. First, two cases with thermal stratification and without film are performed to evaluate the heattransfer wall laws on dry wall conditions. Then, two cases with wall-film are also presented, allowing the final model comparison. The dry-wall and wall-film pairs of cases are each done at two different Reynolds numbers Re τ (defined in the next section) to assess the dependence of mesh refinement at the wall.

Configuration
As illustrated in Figure 1, the standard sizing of the minimal channel flow configuration is applied. The top and bottom boundaries of the computational box are isothermal no-slip walls, while the streamwise and spanwise boundaries are periodic. The gas flow consists of a mixture of air and fuel. The fuel is iso-octane in all the performed simulations, except for the validation case of the local-gradient approach where the fuel is n-heptane to be consistent with the work of [4]. In all computations, the Reynolds number Re τ based on channel half-height H, the wall shear velocity u τ = τ s /ρ s with τ s being the wall shear stress, and the kinematic viscosity ν s , is targeted. CONVERGE uses Cartesian mesh. For the wall-modeled simulations, a uniform mesh grid with cubic cells is used. The mesh cell size in this case is set such that ∆x = ∆y = ∆z = H/36. In the case of wall-refined simulations, the base grid is more refined in the wall-normal direction such that ∆x = ∆z ≈ 2.3∆y. The refinement close to the walls is case dependent and ∆y + 1 is given later for each case.
Numerical method CONVERGE uses a finite-volume method to numerically solve the conservation equations. Spatial discretization is second-order accurate and time stepping is with the semi-implicit Crank-Nicolson scheme. Implicit LES filtering is used for turbulence modeling. Temperature dependent properties of air are used for the molecular dynamic viscosity µ mol and molecular thermal conductivity λ mol of the mixture. The molecular diffusion coefficient of each species k is based on its given Schmidt number Sc k and its molecular kinematic viscosity ν mol following D k = ν mol /Sc k . The σ-model [6] is used to express the turbulent Sub-Grid Scale (SGS) viscosity with a default value of C σ = 1.5. This SGS model is wall adapting and ensures the expected cubic behavior near solid boundaries [6]. The turbulent thermal conductivity and turbulent diffusion coefficients are computed respectively using turbulent Prandtl and Schmidt numbers P r t = Sc t = 0.7. Two momentum wall-laws are compared for the wall-modeled LES of the first isothermal case without film. The standard law-of-the-wall model, referred hereafter as lin-log, and the Werner-Wengle wall-law [7]. The influence of wall heat flux modeling is evaluated through the comparison of four heat flux models: O'Rourke-Amsden [8], Han-Reitz [9], Angelberger [10] and GruMo-UnoMORE [11]. The film initialization procedure was modified in order to ensure the film thickness on the wall is uniform. The film is considered sufficiently thick so that its thickness and surface temperature could be set as constant during the simulation, allowing a statistically steady evaporation configuration.
Finally, source terms are added to the momentum, energy and fuel mass fraction equations in order to compensate for the effects of the wall shear stress, heat flux and mass evaporation. The streamwise momentum source term is adjusted for each case to reach the target Re τ . Similar to the procedure followed by [4], expressions depending on the target values in the center of the channel have been implemented for energy and mass source terms.

Isothermal turbulent channel flow
The quality of the predicted aerodynamic fields is assessed by comparing the results against the reference DNS by Lee & Moser [12]. Simulations were performed at two different turbulent conditions, Re τ = 550 and Re τ = 1000, as summarized in Table 1. As mentioned previously, a wall-refined case and two wall-modeled cases, using the standard lin-log and the Werner-Wengle wall-laws, are evaluated. Figure 3a and Figure 3b show the normalized streamwise velocity u + = u/u τ mean value and its fluctuations u + rms = √ u 2 /u τ at Re τ = 550 and Re τ = 1000 respectively. The wall-refined results are in a good agreement with the results of [12] in both cases. The shape of the normalized fluctuations are well reproduced although they are under-estimated compared to the reference results at Re τ = 1000. This under-estimation can be accepted considering the difference in accuracy between the numerical methods (Lee & Moser DNS code based on spectral method [12]). In the following, the wall-refined configuration will serve as reference. Using the two wall-models, the normalized fluctuations of the three velocity components are similar in shape and scale as the DNS results (only u + rms is shown) at Re τ = 550. However, both fail to capture the pic of fluctuations at Re τ = 1000 as the mesh is not refined enough at the wall. While very close at Re τ = 550, the two momentum wall-laws predictions for u + are significantly different at Re τ = 1000, Werner-Wengle wall-law being much closer to the DNS. From these results, the Werner-Wengle wall-law is retained for the rest of this study.

Validation of the local-gradient evaporation model
In order to validate our implementation of the local-gradient evaporation model in CONVERGE, we use Desoutter et al. [4] DNS of an evaporating film in a turbulent channel flow as a reference. They used the AVBP [13] code with a third-order time and space Taylor-Galerkin scheme [14], along with an evaporation model based on the local-gradient approach (Eq. (1)). In this case the mixture is air/n-heptane. The film surface temperature is set to 309.4 K, the target temperature and fuel mass fraction at the center of the channel are T target = 400 K and Y F,target = 0.1. For a given case conditions, the fuel mass fraction at the film surface Y F,s (ex-  pression 2) is constant in time and space and it is highly dependent on the case pressure and on the used P v,F . In order to avoid any differences related to Y F,s evaluation, we made the choice to impose Y F,s = 0.3 as obtained in the computation of [4]. Desoutter et al. [4] and this work results are summarized in Tab. 2 and Fig. 4. The normalized velocity and its fluctuation profiles, the normalized temperature and its fluctuation profiles (T + = (T s − T )ρ s c p u τ /φ s and T + rms = √ T 2 ρ s c p u τ /φ s ), and the fuel mass fraction and its fluctuation profiles agree relatively well with the reference results. The temperature gradient at the film surface is underestimated by about about 17% leading to an error on the heat flux of about 8%. In our case, the fuel mass fraction in the wall cell center is slightly underestimated compared to the reference result at the same distance from the wall, leading to an over-estimation of the fuel gradient at the film surface (about 23% higher than reference result), and thus, an over-estimation of the evaporation rate by about 25%. Given the differences of mesh refinement at the wall, and given the numerical differences between the two codes, the present results are acceptable and the local-gradient approach is validated. In the following, the wall-refined configuration along with the local-gradient will be used to get reference results for far-field approaches assessment.

Assessment of wall heat flux models: Wall heat-flux
First, the four considered heat wall models are evaluated in a turbulent channel flow with thermal stratification and dry walls. In this case the mixture is air/iso-octane at an equivalence ratio of 0.8. The imposed wall temperature is 373K and source terms are set to keep a mean temperature of 600K at the center of the channel. The different cases first cell distance to the wall ∆y + 1 , reached Re τ and T c , predicted heat flux φ s with its deviation to the wall-refined reference φs are summarized in Tab. 3. Figure 5 shows the normalized profiles of the temperature mean value and fluctuations. At Re τ = 550, the wall-modeled profiles are almost the same and agree very well with the wall-refined profiles. At Re τ = 1100, the wall-modeled normalized mean temperature profiles differ slightly with globally the same slope an magnitude as the wall-refined case. At the first node corresponding to y + ≈ 15, all the wall-modeled cases recover the wallrefined temperature with an error smaller than 10 K. Table 3 presents the obtained heat flux at the wall for each case along with the relative error using the wall-refined case as reference.At Re τ = 550, all the heat models predict very well the wall heat flux. At Re τ = 1100, the errors are more important (about 10%). Han-Reitz and Angelberger predict almost the same wall heat flux as these two models have very close formulations. Han-Reitz heat model will not be considered for the assessment of the film evaporation rate in the following.

Assessment of wall heat flux models: Film evaporation rate
The last step of this study is the assessment of the wall-modeled film evaporation predictions with each heat-transfer wall law. The mixture used is again air/iso-octane. The imposed wall temperature is 370K which is below the saturation of iso-octane at atmospheric pressure (T sat = 372.3K), ensuring to stay in the wetting regime. The source terms are set to keep a mean temperature T target = 600K and a mean fuel mass fraction Y F,target = 0.05053 at the center of the channel. Table 4 summarizes the actual Re τ , mean temperatures and fuel mass fractions at the center of the channel for each case, as well as the heat-fluxes φ s and evaporation ratesṁ and their deviation from the wall-refined reference. All the wall-modeled simulations at Re τ = 550 predict a similar evaporation rate, about 45% smaller than the wall-refined result. At Re τ = 1100, there are slight differences between the modeled cases, around 60% smaller than the wall-refined result. The errors on the evaporation rate estimation are much larger than those of the wall heat-flux estimation, suggesting that the estimation of the heat transfer coefficient is not the main source of errors on the evaporation rate. A poor modeling of the turbulent mixing near the film surface by the far-field approach maybe the main source of errors on the evaporation rate. Figure 6a shows a significant difference on the velocity profiles between the wall-refined case and the wall-modeled cases, leading to 30% and 40% smaller mass flow rates for the Re τ = 550 and Re τ = 1100 cases respectively, while it didn't exceeded 10% for dry wall cases.

Conclusions
In the present study, an academic case is used to evaluate the accuracy of film evaporation modeling based on a far-field approach in conditions that may be encountered in applied multiphysics cases with CFD code using AMR. The dependence to the heat-transfer wall-modeling strategy and to wall cell size was thus focused on. When using such models, the error was found to be around 45% with respect the wall-refined case when the center of the wall cell is located in the viscous sub-layer, and around 60% when the center of the wall cell is located in the buffer layer. In all wall-modeled cases, it was found also the evaporation rate prediction does not depend on the heat-transfer wall model accuracy. It was noticed also that the presence of an evaporating wall-film increases the errors on the mass flow-rate in the channel for the same Re τ . These results suggest that film evaporation has a significant influence on near wall physics and the usage of the classical wall-laws with film evaporation modeling based on a far-field approach is not suited for wall-modeled LES cases. Finally, the results of this study were obtained for fixed values of the empirical exponents n and m (expression (4)), which need therefore to be investigated. In addition, the film thickness and surface temperature were set constant in present simulations. Transient response of these quantities along with the influence of AMR usage will be investigated next.