Modeling evaporation with a meshfree collocation approach

In this paper, a new model for the below-boiling point evaporation process with a meshfree collocation method is developed. In order to capture the phase change process, two different approaches are proposed: multi-phase and single-phase. First, a multi-phase approach is considered, where a novel mass transfer model assumes that the diffusion driven by the vapor concentration gradient in the air phase near the interface is the primary driving force for the mass transfer between phases as both the liquid water and air/vapor phases are simulated. Then, a water-only single-phase approach is also proposed, in which only the liquid water phase is simulated. For this, appropriate free surface boundary conditions are developed based on the convective mass transfer theory to model evaporation and incorporate airflow effects without explicitly simulating the air phase. In order to validate the proposed models, a series of experiments with varying air temperature, relative humidity, and airflow rate is conducted. The numerical results show a good agreement with the evaporation rate measured in the experiments. The multi-phase simulations agree better with the experiments, while the single-phase simulations also produce good results with a much lower computational effort.


Introduction
The below-boiling point evaporation of water can be easily observed in everyday life.However, unlike its occurrence, it is not trivial to accurately model the phenomenon and predict the resulting evaporation rate in a numerical simulation.It is a complex interfacial mass transfer process depending on many different factors [43].Nevertheless, accurate estimation of the evaporation rate and the corresponding total dry-out time under the given conditions are crucial in many applications.One crucial example of this is the prediction of local corrosion effects in automotive applications.Medium-and long-term water accumulation on the metallic surfaces of a vehicle can initiate and accelerate localized corrosion, which eventually harms the life span and durability of the entire vehicle [37,38].Thus, developing a robust and flexible numerical method to resolve evaporation phenomena for the automotive applications can support the development of anti-corrosion measures during vehicle design and, consequently, save time and cost.
In the field of computational fluid dynamics (CFD), most numerical methods require mesh generation over the computational domain to resolve fluid mechanics.While these traditional mesh-based methods such as finite volume method (FVM) thrive in a wide range of applications, many researchers have also developed meshless methods for many other specific purposes over years.One of the first developed and most widely used meshless method is called smoothed particle hydrodynamics (SPH) [30].SPH is a Lagrangian meshfree approach that discretizes the computational domain with mass-carrying particles.It was originally developed for astro-physical applications, where no solid boundaries are present.This implies one of the major drawbacks of SPH.Since it does not have a natural way to impose most of the boundary conditions, it is difficult or needs extra effort to enforce certain boundary conditions [15,44].Another popular meshless methods class is based on radial basis functions (RBFs).RBF-based methods use radial functions, whose values solely depend on the distance from a fixed point.A desired function is approximated by a linear combination of a single radial basis function at different translated points, and the function derivatives are defined by the derivatives of the incorporated radial basis functions [44].Mai-Duy and Tran-Cong [31] and Šarler [40] have deployed this method to resolve Navier-Stokes equations and fluid mechanics in their studies.One weakness of RBF-based methods is that the accuracy and stability of the method highly depend on the choice of the basis functions and their shape parameters [35].Many studies such as [11,20,21] introduce ideas to improve the accuracy and stability/convergence of the method for many engineering and scientific applications.However, they are usually application-specific and ad hoc in nature [44].The meshless approach used for the current study is a meshfree collocation approach referred to as the generalized finite difference methods (GFDM).GFDM generalizes finite difference methods (FDMs) to arbitrary point distributions by using weighted least squares approximations.Here, the desired derivatives to resolve the flow regime are given by linear combinations of the function values on neighboring points [34].It provides a straightforward way to incorporate a wide variety of boundary conditions, which is a considerable advantage over SPH.Moreover, the derivatives in GFDMs are taken as linear combinations of function values.These are not directly dependent on the derivatives of kernel functions unlike in RBF generated finite differences (RBF-FD) or SPH, where the desired function derivatives depend on the chosen radial basis function or kernel function, respectively [44].
In the automotive industry, air-only single-phase CFD simulations have been traditionally used in the fields of aerodynamics and aeroacoustics development for many years.Recently, the role of numerical methods that can efficiently simulate free surface flows to analyze water management processes for vehicles has increased significantly.They enable the optimization of water protection functionalities as well as water drainage systems and potentially also the minimization of corrosion effects.Prominent applications of automotive water management simulations include water wading and fording [32,49], analysis of the impact of rainwater on the windshield, and cleaning of sensors for automated driving.While FVM-based volume-of-fluid (VoF) methods have also recently been used for these simulations [8], meshfree methods are more suitable for these applications due to the ease of handling free surfaces and large domain deformations [23,56].Most of the research related to vehicle water management applications to date focuses on the dynamic behavior of water flowing along complex external vehicle surfaces, e.g., during rainfall or water wading.This merely examines which parts are wet by external water sources, but not how long they stay wet.Thus, under this concrete motivation of extending the scope of CFD simulations in the automotive industries, the present work is proposed as an initial step for the development of numerical methods to expand their water management and corrosion protection capabilities for the vehicle component applications.
Many researchers have developed methods in CFD for multi-phase mass transfer phenomena across a phase interface for many different applications.Wickert and Prokop [55] provide a state-of-the-art overview of already developed water evaporation solvers using several commercial CFD codes.They validated the simulation results with a simple experiment only at one ambient atmospheric condition-35 • C and 32 % relative humidity-without forced convection of air.Furthermore, they analyzed and compared the performance as well as physical assumptions of the different CFD codes.However, detailed analysis of the performance of each commercial code, especially with the existence of forced convection of air flow, is absent in this survey.Vidal-López et al. [52] deployed a CFD simulation technique using Ansys Fluent to analyze the mass transfer coefficients (wind function) in a Class-A evaporation pan, a standardized measurement equipment by the US National Weather Service, and compared the results with empirically derived wind functions.Li and Heiselberg [28] as well as Gallero et al [13] used CFD simulation tools to analyze and estimate water evaporation aspects in swimming pool applications.Li and Heiselberg [28] compared different water evaporation models using Ansys Fluent to determine the mass flow rate of water evaporation and air humidity distributions in a swimming facility for the ventilation system in indoor swimming pools and validated the results with an actual swimming pool facility.Gallero et al. [13] also used Ansys Fluent and proposed an improvement of a formerly proposed evaporation model for swimming pools by [1] in order to enhance the accuracy when forced convection is dominant.Gallero et al [13] Li and Heiselberg [28] and Vidal-López et al [52] all used air-only single-phase simulations by imposing boundary conditions for the water surface instead of defining an air-water interface for the sake of simplicity, which eventually cannot calculate the actual water amount reduced by evaporation.On the other hand, in this current study, the water source for the evaporation is present for both approaches (multi-phase and single-phase approaches) so that the reduction of the actual water amount can be accurately calculated.
Kumar et al [27], Ragab and Wang [36], Saufi et al [39] developed and investigated evaporation models for liquid (water or fuel) droplets.Kumar et al. [27] used the discrete phase model (DPM) in Ansys Fluent to investigate the droplet evaporation behavior for applications such as spray cooling, scramjet combustion, and fuel injection.Saufi et al. [39] developed a comprehensive computational solver in Open-FOAM, called DropletSMOKE++, to resolve the evaporation of fuel droplets using the VoF method.There, the evaporation rate is directly evaluated based on the vapor concentration gradient at the phase interface used in the current study.Ragab and Wang [36] introduced an improved evaporation method to increase the accuracy of the resulting temperature during the evaporation rate and compared the results with the built-in droplet evaporation model in Ansys Fluent.Yang and Kong [57] and Wang et al [54] investigated evaporation models for the liquid droplets using SPH.Yang and Kong [57] proposed a general approach for evaporating multi-phase flows discretized by SPH particles.Here, the mass transfer rate from liquid to gas particle element is calculated based on the vapor concentration gradient as of this study and compared the 2D simulation results with the literature.wang et al. [54] also introduced a 2D SPH framework of the liquid fuel droplet evaporation for the combustion process.They used Fourier heat conduction equation and Fick's law of diffusion for the heat and mass transfer calculation.The single droplet simulation results are validated with the experiment.All these research works for evaporation phenomena ( [27,36,39,54,57]) focus on the analysis of microscopic physical aspects of evaporation for a single small-sized droplet in the range of μm to mm for applications such as spray cooling or fuel injection in combustion engines.The target scale of the evaporation source in these applications is much smaller, and the temperature conditions are much higher than the goal application of the current study.Moreover, even though some researchers such as [39,57] calculated the mass transfer rate based on the vapor concentration gradient at the interface as in the current study, the results were not validated thoroughly in the larger scale also in connection with different ambient conditions.
In this study, we present a new mass transfer model to simulate the below-boiling point evaporation of water implemented in GFDM.Given our overall motivation of proper corrosion protections as well as the ultimate purpose of constructing an automated simulation pipeline for water evaporation on complex vehicle geometry applications without any manual geometry clean-up procedure, which is normally necessary in mesh-based methods, we utilized GFDM to discretize the fluid domain.This approach can easily capture the free surface deformation or disconnection that can occur during the evaporation process.In addition, this makes it possible to keep a sharp phase interface independent of the refinement and geometric location, which are great advantages over mesh-based methods such as finite element method (FEM) or FVM as well.Considering the ultimate automation process for vehicle applications, it also gives more room for potential optimizations with respect to adaptive refinement based on the local humidity and velocity gradients, which is usually more time-consuming in meshbased methods.Another advantage of the meshfree approach is that the free surface of the water can be captured without the need to explicitly simulate the surrounding air-unlike mesh-based methods such as VoF.Thus, evaporation can be simulated either in a multi-phase (water and air) or a singlephase (water-only) scenario, both of which are considered in the present work.Moreover, as explained earlier and in [44], GFDM has several beneficial aspects compared to other possible meshless methods.Besides, the current study extends its proven vehicle water management capability [23] further by allowing seamless migration of previous simulation results from already deployed water wading or rainfall simulations, which can be used as the initial state of the evaporation simulation introduced in this study.In the multi-phase scenario, we assume that the main driving force of the liquid mass transfer into vapor is diffusion due to the vapor concentration gradient in the air close to the free surface of the water [6].Thus, semi-empirical evaporation submodels are not required.We replace this exact treatment of the mass transfer with a suitable approximation in the single-phase scenario, which leads to much faster simulations compared to the multi-phase scenario.
In order to validate both newly implemented approaches, we also conduct a series of experiments with a relatively simple vessel geometry.In the experiments, we measure the evaporation rate at different ambient air temperatures, relative humidities, and airflow rates.While existing experiments [18,19] fixed one variable for each measurement series and measured the evaporation rate by varying the other two variables (e.g., fixed 40 % relative humidity and measured the evaporation rates with different air temperatures and flow rates), all three variables are adjusted and swept in the experiment done in the current study.Furthermore, the existing experimental work typically measures the water temperature at the air-water interface using a thermocouple.In contrast to this, we measure the temperature of the water body in our study, which is then used as initial condition in the simulations for more accurate validation results.
The paper starts with the description of the mathematical models used for the evaporation physics including the underlying conservation equations and mass transfer model (Sects.2, 3).The computational models adopted for the GFDM and the boundary conditions applied for each approach are then explained (Sect.4).In Sect.5, we present the design, build, and conduction of the validation experiments as well as a brief discussion of the results.Following that, in Sect.6, we analyze the simulation results of the two approaches in detail and compare them to the real-world experiments.Finally, the results of this research and specific aspects of each approach are summarized in Sect.7. We also discuss the possible direction future research could take.

Mathematical model: multi-phase flow
We start by introducing the model for the multi-phase approach, where both the water and air phase are simulated.The effects of the bulk fluid motion from the convective mass transfer mechanism at the interface are negligible because of zero relative velocity [4].Thus, in order to accurately model the evaporation physics in multi-phase flows, we assume that the main driving force of the liquid mass transfer into vapor is diffusion only driven by the vapor concentration gradient in interfacial normal direction in the air phase close to the water surface [6].Here, the evaporation rate is evaluated directly without the necessity of any semi-empirical evaporation submodels or assumptions.

Conservation equations
The underlying fundamental governing equations for modeling both the water and vapor phases are the standard conservation equations for mass, momentum, and energy in Lagrangian form given by [23] where ρ is the mass density, p is the pressure, v is the velocity vector, S is the stress tensor, g is the gravitational body force, ρ • E is the total energy, λ is the heat conductivity, T is the temperature, ∇ is the gradient, and ∇ T is the divergence.
From the above general physical model, the stress tensor (S) is split into solid and viscous parts in order to derive a useful numerical scheme as [23] S = S solid + S visc . ( Since the solid part is not under consideration for our specific application, only the viscous stresses are used and defined as where η is the dynamic viscosity (Pa • s) and dε dt is the strain rate tensor, accounting for the local deformation rate of the fluid given by [23] where I is the identity tensor.

Mass transfer rate
The liquid mass is transferred into vapor by the vapor concentration gradient in the vicinity of the interface in the air phase.Thus, the local mass transfer rate or mass flux ṁ water in the unit of kg m 2 •s can be calculated as: where ρ air is the density of air kg m 3 , D vapor is the diffusion coefficient of vapor in air m 2 s , and ∂ ∂n | air,interface represents the vapor concentration gradient at the interface in the direction normal to the interface (m −1 ).This equation is based on Fick's law of diffusion, assuming that the vapor concentration or vapor partial pressure always corresponds to the saturation concentration or saturation partial pressure at the interface for the given (interfacial) temperature [4].
The water vapor concentration (or water vapor mass fraction in air), with respect to the given condition, can be obtained as: where M air , M vapor are the molar masses of air and vapor g mol , ϕ is the relative humidity (%), p vapor,sat is the saturation partial pressure (Pa) of vapor at the given temperature, and p moist_air is the total pressure of air (Pa) or the atmospheric pressure for the simplicity [24].p vapor,sat at the given temperature can be calculated using either the Arden Buck equation or the psychrometric chart [3,41].
The effective diffusion coefficient D vapor of vapor in air is defined as the sum of the laminar and turbulent diffusion coefficient and is given by The laminar diffusion coefficient is given by [17] D vapor,laminar = 2.3056 where p 0 is the atmospheric pressure (Pa), p is the current pressure (Pa), and T is the temperature ( • C).Note that the leading coefficient is adapted from the original version from the literature due to conversion from m 2 h to m 2 s .In the current study, p 0 p is effectively 1; thus, the laminar diffusion coefficient solely depends on the temperature.Using the standard k − ε turbulence model [34], the turbulent diffusion coefficient is given by [51]: where k is the turbulent kinetic energy m 2 s 2 and ε is the turbulent dissipation m 2 s 3 assuming a turbulent airflow.

Energy balance at the interface
Thermal energy at the interface should also be balanced since energy-specifically the latent heat-is consumed during the mass transfer process.Thus, we introduce the following thermal balance equation at the interface by: where h evap is the evaporation enthalpy or latent heat of the water evaporation J kg , λ water is the thermal conductivity of water W m•K , ∂ T ∂n | water,interface is the derivative of water temperature at the interface in inward normal direction K m , and q air to water , q radiation are the heat flux W m 2 from air to water and due to solar radiation, respectively.In this study, we do not consider solar radiation, i.e., q radiation = 0.

Conservation equations
The fundamental conservation equations used for the wateronly single-phase approach are identical to the ones used in Sect.2.1, but are only used here for the liquid water phase.

Mass transfer rate
Since meshfree methods can simulate free surfaces without the explicit presence of the air phase, which is a huge advantage over traditional, mesh-based CFD methods, approximate mass transfer solutions can be obtained with less computational effort using water-only single-phase simulations.
For the single-phase simulation, convective mass transfer theory is used in order to approximate the effect of the bulk air movement as the air phase is not resolved directly.Here, the mass transfer rate (mass flux) is given by: where ρ air is the density of air kg m 3 , α air,mass is the averaged mass transfer coefficient m s , ∞ air is the dimensionless ambient vapor concentration in air, and sat air (T water ) is the dimensionless saturation vapor concentration at the interface for the given temperature T water .This equation is derived from the thermal and mass analogy of Newton's law of cooling [4].
It is not straightforward to obtain an exact, universal mass transfer coefficient α air,mass , and a number of semi-empirical theories and equations to estimate it for various different scenarios exist [4,43].Among them, one approximation can be obtained using Sherwood, Schmidt, and Reynolds number relations (Sh, Sc, and Re, respectively) when replacing the Nusselt and Prandtl numbers in the Nusselt equation as described by [5,22,33] where D vapor is the vapor mass diffusivity in air m 2 s , ν is the kinematic viscosity m 2 s , and L is the length scale of the geometric surface, here the free surface of the water.For the present study, we apply the Sherwood number estimation for forced convection over a flat plate [4], which is a similar scenario to the validation experiments presented in Sect.5, leading to With this, the evaporation rate can be approximated and used as boundary condition at the free surface.

Energy balance at the free surface
In the water-only single-phase simulation, the exact amount of heat flux from air to water and the heat transfer coefficient cannot be obtained accurately.Consequently, the thermal equilibrium during a simulation run cannot easily be achieved.Thus, it is assumed that the consumed latent heat at the free surface is compensated solely by the heat flux from air to the water bulk since equilibrium bulk water and air temperature are imposed as initial conditions and that the equilibrium thermal condition is sustained during the simulation.Furthermore, in the multi-phase simulations and the experiment done by [19], it is observed that the temperature at the free surface of the water is merely 0.2-0.5 • C lower than the bulk water temperature, which would have minimal effects on the results if taken into account.

Computational modeling
In this section, we describe the discretization procedure of the mathematical models introduced in Sects. 2 and 3.The computational framework used in the present work is the simulation suite MESHFREE [12], which is built on a generalized finite difference method (GFDM), and has also been referred to as the finite pointset method (FPM) in our previous work [23].

Domain discretization
The computational domain is discretized with a scattered set of points, referred to as a point cloud, without any underlying mesh topology connecting them.Scattered points are used to discretize the entirety of the domain, including not just the interior, but also domain boundaries and free boundaries.
In the multi-phase approach, a point cloud discretizes each phase and also the interface between the phases.Points of both phases are present at the interface.Examples of point cloud discretizations are shown in Figs. 1 and 2 for the singlephase and multi-phase cases, respectively.All point clouds used in this work are generated using a meshfree advancing front procedure [48].

Derivative computation
Numerical derivatives are computed using a GFDM [10,14,29], a collocation approach, which discretizes the governing equations in their strong formulation.As the name suggests, the GFDM generalizes classical finite differences to scattered point clouds without regular placement.At a numerical point i, located at x i , all derivatives are computed as weighted sums of function values at neighboring points.Thus, the derivatives of a function u are approximated as: where ∂ * is the continuous derivative being approximated such as the x-, y-, z-derivative or Laplacian operator.The numerical derivative at point i is represented by ∂ * i .The coefficients c * i j are determined using a weighted least squares method.We refer to [44] for more details on the derivative computation procedure and to [25] for the use of GFDMs for phase change modeling.

Time integration
The time-integration procedure starts with a discretization of the Lagrangian motion by moving each point with its velocity [45].To update the velocity and pressure, we use an implicit projection-type approach [2,7].All resultant implicit linear systems are solvers with a BiCGSTAB solver [53] without the use of any pre-conditioner.For the multi-phase approach, a monolithic approach is used, in the sense that both phases are solved together in one sparse system.
For an in-depth explanation of the time integration scheme, we refer to our previous work: [23,47] for the overall scheme, [45] for an in-depth discussion on the Lagrangian motion step, and [34] for the discretization of the turbulence model.

Interface
For the multi-phase simulations, an example of the initial domain discretization can be seen in Fig. 2. The figure illustrates the distinct regions of the air and water phases.Points of both phases are present at the interface, as shown in Fig. 2.This enables us to enforce interface the boundary conditions for each phase.For a coupling between the phases across the interface, each interface point is associated with an opposite interface point of the other phase.Consider a point j water of the water phase.At each time step, this point is associated with a point i air of the air phase, which is the interface air point closest to j water .To determine a physical quantity of interest, for example the vapor concentration, of the air phase near j water , the corresponding value at i air will be used.A similar concept of opposite point is also used for the air phase to determine water phase quantities nearby.The actual boundary conditions applied at the interface are elaborated in Sect.4.5.1.

Multi-phase approach
For the multi-phase approach, boundary conditions for the air should be defined to limit the computational domain.This can be thought of as defining an air box, whose inlet boundary conditions are identical to the ambient air conditions.In detail, at the inlet, constant flow rate, and temperature are imposed as conventional Dirichlet boundary conditions.Furthermore, a desired constant ambient humidity or vapor partial pressure is defined in terms of vapor mass fraction, calculated based on the equations explained in Sect.2.2, as the inlet boundary condition.At the outlet, on the other hand, conventional zero-gradient Neumann boundary conditions for these variables are applied.Furthermore, for the pressure, standard Neumann zero-gradient conditions for the inlet while atmospheric fixed pressure value with Dirichlet condition for the outlet boundary conditions are applied.
For other geometric elements, normal no-slip, adiabatic wall, and standard zero-gradient wall pressure conditions are imposed and no vapor concentration interaction is defined.For the water side, since there is no water source defined, simply normal no-slip, adiabatic wall boundary conditions are defined, where the geometric elements are in contact with the water phase.
In the phase interfacial region, which is the most interesting part of this study, the mass and thermal transferring boundary conditions are defined for both phases.First, the water mass is reduced based on the mass transfer rate calculated with the model described in Sect.2.2.The air-side boundary points, on the other hand, hold the saturated vapor concentration for the given interfacial temperature, which would feed the vapor to the air phase just like a source term.Second, based on the mass transfer rate, the consumed latent heat is deducted as described in Sect.2.3 with a negative heat flux source term at the water side boundary condition.Additionally, the heat exchange between air and water based on the temperature difference on each side is also applied at boundary conditions of both sides.

Single-phase approach
For the water-only single-phase approach, since the air phase is not explicitly simulated, the limiting boundary element definition for the air (or air box) is not required.Merely the boundary elements that are in contact with the water phase should be defined as conventional no-slip for the velocity, adiabatic walls for the thermal, and standard zero-gradient Neumann wall for the pressure boundary conditions.Here, the mass transfer rate is simply calculated at the free surface as described in Sect.3.2, and the corresponding mass is reduced as a negative source term at the boundary points.

Free surface correction
Water mass transfer due to evaporation is not directly incorporated into the boundary conditions, but added as a separate step after the standard free surface movement is calculated for both approaches-single-phase and multi-phase.For each free surface point, its additional normal velocity is calculated from the given mass or volume change rate and the free surface area that this point represents.Then the point is moved according to this velocity and the current time step size.During the following point cloud organization step, interior points that may now be too close to the free surface will be merged with the free surface points to maintain a good point cloud quality.Since each free surface point moves independently based on the calculated evaporation rate, the current method can inherently handle the free surface effectively even when the water volume is disconnected during the mass transfer process over time or when discrete water volumes already exist from the initial state.Furthermore, the setup includes a time lapse factor, which can increase the speed of this additional surface movement due to evaporation, such that the evaporation effects can be sped up virtually compared to the airflow and the overall simulation time can be reduced.
Meshfree collocation methods like GFDMs are not inherently volume-or mass-conserving for cases with free surfaces [46].Thus, several options for volume/mass correction are implemented in the software MESHFREE [50].For each of these, a target volume/mass is tracked, accounting, e.g., for changes through inflow or outflow boundaries or points deleted by user-defined criteria.This tracking had to be extended to the changes caused by the additional surface Fig. 3 Schematic view of the evaporation rate measurement experiment.Detailed dimensions of the actual experiment geometry are demonstrated in Sect.5.1 as well as in Fig. 8 movement due to the given volume or mass change rates.In the simplest volume correction case that is used for this study, the free surface points are moved in normal direction again by an average distance determined by the volume error (up to a maximum relative to the total volume) divided by the total free surface area.This works well as long as the total water volume is in one connected domain.For scenarios with multiple disconnected water volumes, such as several distinct droplets, more localized volume or mass correction approaches need to be considered.
A detailed description of conservation aspects and remedies for meshfree methods goes beyond the scope of the present paper and therefore will be the subject of our upcoming work.

Experimental measurements
In order to validate the newly developed evaporation model in terms of the evaporation rate, a series of experiments were performed under different conditions.Here, three variablesambient air temperature, ambient air humidity, and airflow rate-that highly influence the evaporation rate were varied.

Experimental setup
A schematic of the validation experiment can be found in Fig. 3. First, a channel of length 35 cm with a cross section of 14 cm × 14 cm was built in order to generate constant airflow by adjustable exhaust fans placed at the channel outlet.Four EBM Papst 8212 JN exhaust fans were placed at the channel outlet, see Fig. 4. At the channel inlet, a coarse mesh plate (approximately 3 cm × 3 cm) and a fine mesh fabric (less than mm scale) were placed back to back as shown in Fig. 5. Thus, we minimized turbulent effects and obtained constant inlet velocities.Second, a glass petri dish with inner diameter of 80 mm was placed as a water vessel in the middle of the channel, on top of a metal cylinder of height 35 mm and with diameter of 50 mm, see Fig. 6.In this vessel, there were initially 40 g of water.The mass difference between the initial state and the leftover water mass after the experiment was measured to determine the evaporated mass (accuracy of the scale is 10 mg), and together with the experiment duration the evaporation rate mg min .Three temperature/humidity sensors were placed along the channel, see Fig. 6.The sensor located closest to the inlet measured the inlet air temperature (accuracy of 0.5 • C) and humidity (accuracy of 1.8 %).These values were used as inlet boundary conditions for the corresponding simulations.We used the middle and outlet sensors to observe the humidity distribution change caused by the evaporation from the water vessel.Lastly, we placed a thermocouple inside the water vessel to measure the temperature of the water body in equilibrium state.These measurements were also used as initial conditions in the simulations.The entire channel was placed in a 1000 liter climate chamber (Weiss WK1 1000), as shown in Fig. 7. Thus, the ambient air temperature and humidity inside the channel can be controlled as desired.The actual airflow rate was measured with a flowmeter (accuracy of 0.1 m s ) that was placed behind the inlet mesh structure inside the channel.Each flow rate was generated with a fixed voltage and current imposed by a voltage/current regulator.
The target conditions for the different experiment runs are described in the following list.Note that the conditions in actual measurements are not exactly identical to the target values but show slight deviations since the climate chamber allows a difference between the target and actual values.For each experiment run, the actual values are described in Sect.6.
A total of 36 experiment runs were performed, all for an experiment duration of 90 minutes.Note that for chamber stability reasons, 30 % relative humidity was used instead of 20 % for the 10 • C cases.

Experimental results
As expected, the results of the experiments show that the evaporation rate increases when the inlet air temperature is increased as well as when the inlet airflow rate is increased.We observe that the evaporation rate decreases linearly when the relative humidity increases.With respect to the airflow rate, the experiments show a quasi-linear relation.The linear slope that prevails for higher temperatures (at around 32 • C) decreases for lower temperatures.Thus, the effects of the airflow rate are less for lower temperatures.Furthermore, the evaporation rate increases exponentially when the ambient air temperature increases.Since the saturation vapor pressure of water increases exponentially with temperature, this behavior has been observed in the literature as well [5].
Finally, the equilibrium water body temperature is lower when the evaporation rate is higher for all given conditions.

Numerical results and validation
Both air-and-water multi-phase and water-only single-phase simulations are performed using a surface triangulation of a CAD geometry that is constructed to match the experimental setup described in Sect.5, see Fig. 8. Furthermore, we use the ambient air temperature, relative humidity, and airflow rate conditions listed for the experiments.In each simulation case, the required time-step size is automatically calculated with Courant-Friedrichs-Lewy (CFL) conditions with the CFL number of 0.2.The resulting time-step sizes are in the order range from 1.0e −4 to 1.0e −5 s depending on the inlet velocities and the resulting velocity changes over time within a simulation run.

Multi-phase
For the multi-phase simulations, both air and water phases are present to obtain accurate evaporation rates by taking into account the effects of the airflow at the phase interface.
A verification and validation of the basic underlying scheme and implementation to simulate continuous flow can be found in our earlier work [9,42,47].
Furthermore, at the phase interface, surface-surface coupling is used since the two phases interact by a clearly defined common interface, and the most important aspects of the evaporation process take place at this common interface [26].The coupling free surface is defined by a pre-generated closed interfacial geometry for the initialization.At the below the interfacial geometry, water points are filled, whereas air points are organized above it.As an example of a simulation result, the evolution of the humidity profile for 1.0 m s airflow rate, 32.2 • C air temperature, and 78 % relative humidity is illustrated in Fig. 9 as well as in the animation provided in the supplementary material.In order to correctly resolve the evaporation at the water-air interface as well as turbulent effects in the air, we refine the point clouds near the interface as well as the sensors.We observe that the generated vapor is well diffused from higher-concentrated to lowerconcentrated region as well as convected by the airflow from channel inlet to outlet.We determine the numerical evaporation rate by integrating the calculated mass flux kg m 2 •s over the entire interface with respect to the associated area of the interface points.The resulting evaporation rate over time is illustrated in Fig. 10.The evaporation rate reaches a quasi-steady state of 14 mg min after 0.7 s.The small fluctuations (± 5 % around the mean value) are presumably due to the Lagrangian nature of the simulation.Note that finer point clouds reduce the fluctuations and also lead to a quasi-steady state.For details regarding mean evaporation rate convergence with respect to the point cloud resolution, see Fig. 18 in Sect.6.4.

Single-phase
In the water-only single-phase simulations, as a constant mass transfer coefficient (see Sect. 6.2.2) as well as constant ambient air temperature/relative humidity values are imposed in each case, we observe identical evaporation rates for all free surface points, see Fig. 11.Here, the negative value indicates a mass reduction at the free surface.As in the multiphase simulations, we determine the overall evaporation rate by integrating the local evaporation rate at each point over the entire free surface area.We compare the results of both approaches in the following section.

Evaporation rate validation
In this section, we compare the numerical with the experimental evaporation rates.

Multi-phase
A detailed comparison of the results of the experiments with those of the multi-phase simulations can be found in Figs. 12, 13, and 14 for the three airflow rates 1.0 m s , 1.8 m s , 2.7 m s , respectively.In each figure, the x-axis represents the ambient relative humidity, and the y-axis represents the resulting evaporation rate.The blue solid lines and marker points represent the experimental results, and the red dashed lines and marker points represent the simulation results.The three ambient air temperatures are illustrated by different types of marker points-triangular markers for 32.3 • C, square markers for 20.5 • C, and circular markers for 10.6 • C. The temperature of the water body, which is used as initial condition in the simulations, is noted in brackets next to each measurement point.The error bars for the experimental results are generated based on the accuracy of the measurement equipment and the chamber specifications.Note that the ambient air temperatures in the experiments are subject to a measurement accuracy of ± 0.5 • C.
Within the expected measurement errors, the simulation results show a good agreement with the experimental data.For all temperatures, the trend line with respect to relative humidity fits well to the one of the corresponding experiment.A quantified comparison between simulations and experiments is provided in Sect.6.3.

Single-phase
First, we determine the averaged mass transfer coefficient α air,mass using the Sherwood number estimation as described in Sect.3. We use the radius of the water vessel as length scale, i.e., L = 40 mm.Furthermore, we assume a constant value of Sh = 0.7 for the Schmidt number [4,16].For the Reynolds number, we obtain Re = 2222 for 1.0 m s airflow rate, Re = 4000 for 1.8 m s airflow rate, and Re = 6000 for 2.7 m s airflow rate according to Equation ( 13).This results in  where the temperatures specified in the brackets correspond to the ambient air temperature.Note that the resulting mass transfer coefficient depends on the air temperature because the vapor mass diffusivity D vapor depends on the temperature.
The evaporation rate results for the single-phase simulations are illustrated in Figs. 15, 16, and 17.The numerical results show a close correlation with the experimental ones.However, the evaporation rate is over-estimated for higher temperatures with lower relative humidity.Furthermore, for higher temperatures, we observe a steeper decrease as relative humidity increases than in the experiments.This is a current limitation of the single-phase approach.A quantified comparison between simulations and experiments is provided in Sect.6.3.

Error comparison
For both the multi-phase and the single-phase approaches, the relative deviation of the numerical evaporation rates from the experimental ones is summarized in Tables 1, 2, and 3.For a better interpretation of the numerical results, the measurement uncertainties in the experiments are also listed.
Here, the positive sign stands for a higher evaporation rate compared to the corresponding experiment, and the negative sign stands for a lower evaporation rate, noting that highlight means the simulation result is out of the measurement uncertainty range.For high relative humidity as well as low temperature, we observe higher relative errors since the absolute evaporation rates are relatively small and even a small absolute deviation may lead to a big relative error.Furthermore, in these cases, the ambient relative humidity inside the climate chamber has bigger fluctuations in the experiments because it is hard for the chamber regulator to keep stable conditions.This leads to larger errors in the measurements.

Convergence study
For each simulation case, the required point cloud resolution and the number of points are determined by a convergence study.For the sake of brevity, we present only one such convergence study as an example of the performed point cloud refinement study in Fig. 18, namely for the 1.8 m s airflow rate, 20.5 • C air temperature, and 57.3 % relative humidity.The x-axis represents the total number of points in the     The highlighted values indicate that the simulation error is larger than the measurement uncertainty computational domain, and the y-axis represents the resulting evaporation rate.The value used as the minimum point cloud resolution (also known as smoothing length or interaction radius) around the water-air interface is noted in brackets next to each marker point.The evaporation rate is nearly constant for a minimum point cloud resolution of 1.6 mm and shows less than a 2.0 % difference compared to the 1.5 mm case and even less difference between 1.5 mm and 1.4 mm.Thus, we conclude that a converged state with respect to the point cloud resolution is reached.The identical convergence study method was performed for each measurement point in order to obtain a point cloud-independent solution.

Runtime comparison
We present a comparison of the runtimes for both approaches in Table 4.Note that all simulation runs were performed on identical hardware with the same number of cores to primarily evaluate the relative runtime difference.However, the hardware specification is confidential.The runtime corresponds to the total clock time for the entire simulation, which also includes the reading of the bounding geometry, initialization and management of the point clouds, as well as postprocessing to determine mass fluxes and evaporation rates.Note that the comparison is only made with respect to the three airflow rates because the runtime depends highly on the airflow rate, and we observe only negligible effects due to changes in air temperature and humidity.Thus, the simulation runs use the same air temperature and humidity, namely 20.5 • C and approximately 60 %.As shown in the table, for the multi-phase approach, the runtime increases if the airflow rate increases.There are two reasons for this: First, larger airflow velocities result in lower time steps due to a CFL-like restriction in the Lagrangian movement of points [45].Second, we observe that a larger number of points are necessary to correctly resolve the flow physics as well as the concentration gradients at the water-air interface, especially for the air phase.On the other hand, in the single-phase approach the required point cloud resolution is only slightly dependent on the airflow rate.Overall, it is sufficient if the point cloud is dense enough to correctly resolve the free surface.However, for an accurate runtime comparison, we use the same point cloud resolutions for both the multi-and single-phase approaches for the same airflow rates.Note that the total number of points presented in Table 4 includes both air and water phase points for the multi-phase simulations.Thus, even for the same point cloud resolution, the single-phase simulations have a significantly lower number of points.In summary, as expected, the runtime and error rate comparisons show that the multi-phase simulations are much more computationally intensive (factor 20-40), but resolve the evaporation phenomenon more precisely.

Conclusion and outlook
In this paper, we investigated the below-boiling point evaporation phenomenon of water.Using a meshfree collocation approach, or GFDM, we introduced a new evaporation model in two versions.In the multi-phase approach, we developed a new mass and heat transfer model based on the vapor concentration gradient.Furthermore, we presented a simplified single-phase approach using convective mass transfer theory.To validate the two approaches, we designed and built an experiment to measure the evaporation rate for different airflow rates, ambient temperatures, and ambient relative humidities.As expected, in the experiments the evaporation rate increases if the airflow rate or the ambient air temperature are increased, whereas it decreases if the ambient relative humidity is decreased.Within the expected measurement errors, the corresponding multi-phase simulation results show a good correlation to the experiments.We observe somewhat larger but acceptable deviations in the single-phase simulation results.The multi-phase approach is advantageous over the singlephase one because it can accurately determine the local evaporation rate, taking into account the effects of local humidities and airflow characteristics.However, it is much more computationally intensive than the single-phase approach since both phases and the interactions at the water-air interface have to be resolved.The single-phase approach, on the other hand, is much simpler and computationally lighter.It can also be easily applied to complex geometries without the need to maintain phase coupling of air and water phases.However, the mass transfer coefficient is calculated based on simplified model assumptions.
As the newly developed evaporation features and the performed validation results are promising, the initial steps for the ultimate application goal are successfully achieved.In our future work, we aim to apply both models to actual automotive industrial applications.For this, we are currently focusing on the development of customized adaptive refinement criteria to automatically refine the point clouds based on the vapor concentration and velocity gradients in the multi-phase approach.This enables the application of the multi-phase approach to any complex geometry without hassle.Another line of work is the investigation of the simplified single-phase approach for complex geometries.
Here, the focus is on the general determination of the mass transfer coefficient with an increased accuracy.Moreover, an automation pipeline process including a seamless migration of previous raining or wading simulations performed by the identical GFDM is also included in the upcoming task.

Fig. 1
Fig. 1 Domain discretization with a point cloud for the single-phase simulations.The color indicates different types of points: interior points (blue), free surface points (red), and wall boundary points (white).Only a clip of the full point cloud is shown.(Color figure online)

Fig. 2
Fig. 2 Domain discretization with a point cloud for the multi-phase simulations.The color indicates the different phases: air phase (blue) and water phase (red).Points of both phases are present at the interface, marked with lighter shades of blue and red, respectively.(Color figure online)

Fig. 4
Fig. 4 Experimental setup: Four EBM Papst 8212 JN axial fans at the channel outlet to generate constant airflow

Fig. 5
Fig. 5 Experimental setup: Mesh structures at the channel inlet to reduce turbulence

Fig. 6
Fig. 6 Experimental setup: Structure inside the channel

Fig. 9 Fig. 10
Fig. 9 Multi-phase simulation result: Relative humidity profile from 0 s to 5 s in the cross section of the domain for 1.0 m s airflow rate, 32.2 • C air temperature, and 78 % relative humidity.An animation of the simulation result is provided in the supplementary material.(Multimediaview).(Color figure online)

Fig. 11
Fig. 11 Water-only single-phase simulation result: Constant evaporation rate at the free surface points over time indicated by the color scale-identical evaporation rate profile for both 1 s and 5 s.(Color figure online)

Fig. 13
Fig. 13 Multi-phase simulation validation: Numerical vs. experimental evaporation rate for the 1.8 m s airflow rate.The bracketed values next to the marker points indicate the equilibrium temperature of the water body air : 32.2 • C) Multi-phase simulation (T air : 32.2 • C) Experiment (T air : 20.5 • C) Multi-phase simulation (T air : 20.5 • C) Experiment (T air : 10.6 • C) Multi-phase simulation (T air : 10.6 • C)

Fig. 15
Fig. 15 Single-phase simulation validation: Numerical vs. experimental evaporation rate for the 1.0 m s airflow rate.The bracketed values next to the marker points indicate the equilibrium temperature of the water body

Fig. 16
Fig.16 Single-phase simulation validation: Numerical vs. experimental evaporation rate for the 1.8 m s airflow rate.The bracketed values next to the marker points indicate the equilibrium temperature of the water body : 32.2 • C) Single-phase simulation (T air : 32.2 • C) Experiment (T air : 20.5 • C) Single-phase simulation (T air : 20.5 • C) Experiment (T air : 10.6 • C) Single-phase simulation (T air : 10.6 • C) air

Table 1
The highlighted values indicate that the simulation error is larger than the measurement uncertainty 123

Table 2
The highlighted values indicate that the simulation error is larger than the measurement uncertainty • C air temperature, and 57.3 % relative humidity.The bracketed values next to the marker points indicate the minimum point cloud resolution around the water-air interface

Table 4
Runtime comparison for 5 s of simulation time for different airflow rates, 20.5 • C air temperature, and approximately 60 % relative humidity