Pore scale modelling of moisture transfer in building materials with the phase field method

This study explores the applicability of the phase field method for modelling moisture storage and transport in porous materials. Accordingly, the system is treated as a continuum where the phases (liquid and humid air) are separated through a diffuse interface, which evolves in the pores until the equilibrium state is reached. The interface thickness is related to the surface tension, while the contact angle is defined as a boundary condition. The mass transfer in the porous matrix is driven by the Cahn-Hilliard equation and the phase transition is controlled by an equation of state. The above method is tested for a simple geometry (infinitely extended parallel plates), by comparing the numerical outcomes against available measured data and analytical solutions. The challenges arising for a further application to complex pore structures and real building materials are discussed.


Introduction
Storage and transport of moisture in porous building materials are often modelled via empirical material functions such as the moisture retention curve and moisture permeability. This approach is typically inapt to capture microscopic absorption and desorption processes. On the other hand, pore-scale simulation can help to explain the correlation between storage/transport of moisture and physics parameters such as the contact angle and pore geometry. Pore-scale approaches lead hence to a better understanding of these phenomena and in the future might be even applied for the characterization of material functions in alternative to or in combination with experiments.
In prior studies [1] [2], this issue was addressed by using pore-network models, which however are affected by strong approximations, concerning for instance the geometry of the porous structure and the dynamics of liquid-gas interfaces through the pore junctions. Here we shall follow a different approach based on the phase field method. Accordingly, the system is treated as a continuum where the phases (liquid and humid air) are separated through a diffuse interface, which evolves in the porous matrix until the equilibrium state is reached. The interface thickness is related to the surface tension, while the contact angle is defined as a boundary condition at the solid surface. The mass transfer in the porous matrix is driven by the Cahn-Hilliard equation [3], in which the liquid and vapour flows are described through a Darcytype model and Fick´s diffusion, respectively, while the phase transition is controlled by an equation of state. The above method is tested for a simple geometry (infinitely extended parallel plates), by comparing the numerical outcomes against available measured data and analytical solutions. First, it is shown that the Young-Laplace and Kelvin equations are fulfilled at steady-state conditions. In this regard, the model sensitivity on the interface width and capillary dimension is investigated. Hence, time-dependent isothermal wetting and drying processes are addressed, showing that the model captures well the liquid-gas interface dynamics.
In the last part of the paper, the challenges arising for a further application to complex pore structures and real building materials are discussed. The main limitation appears to be the high computational cost deriving from the extremely thin interface. Nevertheless, the method seems promising for increasing the understanding on pore-scale wetting and drying dynamics.

Mathematical model
In this section, we describe the phase field method used for modelling two-phase isothermal moisture transfer in a porous medium.

The Cahn-Hilliard equation
According to the phase field method, the fluid is considered as a continuum, in which the different phases (liquid and humid air) are separated by a diffuse interface. The transport of mass is described via the Cahn-Hilliard equation: Here 0 denotes the bulk chemical potential and the term 2 2 represents the interface contribution. In particular, the coefficient [m 7 /(kg s 2 )] affects the interface width and surface tension as shown below in section 3. Eqs. (1) and (2) form a system of coupled partial differential equations in which the dependent variables are and . Note that, sufficiently far from the interface, Eq. (2) implies = 0 and Eq.(1) turns into the wellknown diffusion equation. The bulk chemical potential, density and pressure are related via the following equations [4]: where denotes the Helmholtz free-energy density [m 2 /s 2 ], and the pressure [Pa]. Accordingly, one can write: where = 1/ denotes the specific volume [m 3 /kg]. If the function ( ) is known, one can determine the Helmholtz free-energy density as follows: The required correlation between pressure and specific volume ( ) is called equation of state.

Equation of state
The equation of state (EoS) governs the liquid-gas phase transition by stating the interfacial energy that arises from the molecular interaction between the two phases inside the transition region. In turns, it affects the interface thickness and the surface tension, and it equally controls the saturation pressure and the densities of saturated liquid and vapour. The choice of a proper EoS is thus a key aspect of the phase field model. An extended review of available formulations is reported in [5]. In prior studies on phase field modelling [4], the well-known van der Waals equation has been used: with: Here, R denotes the gas constant [kg/(mol K)] and the molar mass [kg/mol]. The parameters a and b in Eq. (7) are determined from the critical values of temperature and pressure, as shown, e.g, in [5]. Again in [5] the authors point out that Eq. (7) is inaccurate for polar compounds as water. For this reason, we shall consider also the analytical correlation proposed by Patel [6], suitable also for polar molecules: The parameters characterizing water at 20°C are reported in Table 1 for both Eqs. (7) and (9). In Table 2 the values of the saturation pressure and saturation densities of liquid and vapour are compared with measured values from the literature. It can be observed that Eq. (9) describes well the behaviour of water at 20°C, even if some small deviations are present, especially concerning the density of the saturated liquid. This formulation has hence been used for the test cases reported below. For higher accuracy, the EoS reported in [7] can be used, which however presents much higher complexity as Eq.(9). In Fig. 1 the pressure and bulk chemical potential calculated via Eq.(9) are shown as functions of the density.

Determination of the mobility
In this section, we describe one way to estimate the mobility appearing in Eq. (1). To this aim, we shall first consider the transport of the homogeneous liquid phase. Neglecting the gravity force and inertia and assuming that viscous forces are proportional to the bulk velocity of the fluid, the momentum equation reduces to a Darcy-type flux: the plates dH=2h and fRe=96 [8]. Eq. (10) becomes thus: Considering that in the homogeneous phases ≈ 0 and taking into account Eqs. (3) to (6), the Cahn-Hilliard flux is written as follows: Hence, comparing Eq. (11) with Eq. (12) at liquid saturation conditions we obtain: sw = 2 ℎ 2 12 (13) Where sw and sw denote the density and the mobility of saturated liquid, respectively. From Eq. (13), assuming for instance h=10 -8 [m] and =1.5 10 -3 [Pa s] [9], one obtains sw =5.9 10 -9 [kg s/m 3 ]. Let us now consider the flux of water vapour due to Fick´s diffusion in the gaseous phase (humid air) [1]: Where 0 represents the diffusivity of vapour in air [m 2 /s], the Knudsen number, the gas constant for water vapour [J/(kg K)] and the absolute temperature [K]. According to [1], The Knudsen number is defined as = /(2ℎ), with the free path length of the gas molecules =5 10 -8 [m]. Comparing Eq. (12) with Eq. (14) at saturation conditions one obtains: Where sv and sv denote the density and the mobility of saturated vapour, respectively. Assuming T=20 [°C] and N k =2.5 (corresponding to h=10 -8 [m]), it results sv =9. 5 10 -13 [kg s/m 3 ].
The mobility in the interface region is assumed to be a linear interpolation of the mobility of saturated liquid and vapour as follows: Note that for = and = Eq. (16) leads to = and = , respectively.

Boundary conditions
For the cases treated below in this study, the following boundary condition are considered. 1) Wetted surface: Eq. (17) imposes the zero flux condition through the surface, while Eq. (18) defines the contact angle ; n denotes the unit vector perpendicular to the surface and pointing outward of the domain. 2) Open boundary (Dirichlet boundary conditions): Here, is the imposed density value at the open boundary. One might impose, for instance, = or = , whether the domain is bounded by liquid water or unsaturated vapour, denoting with the relative humidity at the boundary. 3) Symmetry condition: A symmetry plane is defined via Eqs. (21) and (22), by imposing zero-flux conditions at the boundary.

Testing the model
In this section, the model proposed above is applied to two-dimensional test cases. In particular the following conditions, typically occurring in building physics applications, are addressed, i.e., 1) liquid-gas equilibrium inside a pore, 2) capillary absorption from a liquid reservoir and 3) vapour absorption/desorption from/to the environment. For the sake of uniformity, in all considered cases the contact angle = 30° is applied. According to [10], this assumption is accurate for liquid absorption from a liquid reservoir. The approach is however valid for arbitrary contact angles.
The independence of the results from the mesh refinement has been tested. Cartesian and free triangular meshes have been considered, by obtaining similar results and solver performances. The numerical outcomes are compared with analytical solutions, when available.

Equilibrium conditions for a liquid drop entrapped inside a closed cavity
The model has been first tested at equilibrium conditions, by investigating the behaviour of a liquid drop entrapped inside an infinitely extended cavity with rectangular section, as shown in Fig. 2 (a). Just one quarter of the cavity has been modelled, considering the symmetry planes characterizing the problem. The equilibrium is reached as the final state of a transient simulation, starting with plane liquid-gas interface. It has been observed that the value of the mobility M does not influence the equilibrium state, which is defined by the equation of state and parameter k only.
The simulated density distribution at steady-state conditions is shown in Fig. 2 (b) for two different values of the parameter k. According to prior studies [4], the interface width δ becomes larger with increasing value of the parameter k. This fact is also visible by plotting the density and pressure profiles along the symmetry plane of the cavity (Fig. 3).
Note that the interface width might be defined in different ways, as discussed in [11]. Here, we report the original definition by Cahn and Hilliard [3]: where denotes the critical density.
In Table 3, the values of δ and σ obtained for different values of the parameter k are compared with a measured outcome from the literature [13]. It appears that the measured surface tension can be well predicted with k=0.12 10 -16 [m 7 /(kg s 2 )], while k=1.54 10 -16 [m 7 /(kg s 2 )] is required to fit the measured interface width. In order to achieve a good agreement for both surface tension and interface width further improvement of the equation of state is needed. This task overcomes however the scope of the current study.
In the test cases reported in the following sections (liquid absorption from a liquid reservoir and vapour absorption/desorption from/to the environment) a realistic interface thickness is maintained. Let us now define the capillary pressure as follows: In Eq. (25) and denote the pressure in the liquid phase and the partial pressure of water vapour in the gaseous phase, respectively, which can be determined as output of the simulation. As one notes in Fig. 3 (b), depends on the parameter k, and in turn, for the above considerations, on the surface tension. Moreover, the simulated capillary pressure must fulfil the Young-Laplace equation, which for the given geometry, is written as follows: Here, denote the contact angle and ℎ the cavity width. As shown in Fig. 4 (a) the values of the capillary pressure obtained by numerical simulation are in good agreement with Eq. (26), whether the same surface tension is applied. Further analysis has been performed on the relative humidity of the gaseous phase. In Fig. 4 (b) the simulation results are compared with the values of the relative humidity calculated via the Kelvin-equation: The relative humidity has been determined by numerical simulation in two different ways: 1) according the definition of relative humidity ( = ⁄ ); 2) via Eq. (27) by employing simulated -values. Again, a good agreement between numerical and analytical results can be observed.

Capillary absorption from a liquid reservoir
This section concerns modelling of capillary water uptake from a liquid reservoir. The exemplary case of infinitely extended parallel plates is considered. The following parameters have been applied for the simulation: k=1.54 10 -16 [m 7 /(kg s 2 )] and h=10 -8 [m]. In Fig. 5, the problem geometry and the resulting density distributions at different times are shown. Just the half of the cavity is modelled, considering the symmetry plane characterizing the problem.
The density and pressure profiles in the symmetry plane are reported at different times in Fig. 6 (a) and (b), respectively. In accordance with the above made assumption of a Darcy-type flux, linear pressure profiles can be observed in the liquid phase. The evolution of such linear region is shown in Fig. 6 (c). It appears that the liquid pressure at the border of the linear region decreases, approaching the steady-state value obtained in the previous section. Note that the position of the liquid-gas interface y * can also be determined analytically. Considering Eq. (11) in the one-dimensional case, with j w =ρ dy * /dt and dp/dy=p c /y * , one obtains by integration: * = −√ ℎ 2 6 (28)  In Fig. 7, the result by numerical simulation is compared with the curve obtained with Eq. (28) in which the steady-state capillary pressure has been employed. A good agreement can be observed. The delay of the simulated trend is a combined effect, due to the meniscus formation and the lower pressure drop over the interface observed in the transient water-uptake case.

Capillary drying and vapour absorption
In this section, we shall investigate the case of a liquid drop entrapped between parallel plates and subject to evaporation or vapour condensation, as shown in Fig. 8. The two cases differ only in the value of the relative humidity imposed at the open boundaries (top and bottom of the cavity). In fact, evaporation or condensation occurs depending on whether the relative humidity imposed at the open boundaries is respectively lower or higher than the equilibrium value defined by the Kelvin equation (27). In both the investigated cases, a numerical issue arose concerning the steep change of the mobility over the interface. To overcome this challenge, a reduced mobility value of sw =1.25 10 -10 [kg s/m 3 ] has been applied in the liquid phase, while in the gaseous phase the value calculated above in section 2.3. is maintained ( sv =9. 5 10 -13 [kg s/m 3 ]). The fact that a reduced mobility in the liquid phase has been used shall not affect significantly the results, since the process is mainly controlled by vapour diffusion.
In Fig. 9, the evolution in time of a drying drop is shown. In this case, the relative humidity imposed at the upper boundary is = 0.3 as against the equilibrium relative humidity ≈0.7 reported in Fig. 4  During the first drying period (t<7 10 -8 [s]), the liquid phase fills the entire cavity width, with the formation and progressive regression of the meniscus. Afterward, the liquid is confined on the sidewalls, while the central part of the cavity is already empty.
The opposite case of a liquid drop subject to vapour condensation is reported in Fig. 11. Here, the relative humidity imposed at the open boundary is = 0.95. These examples show that the model can capture quite well the liquid-gas interface dynamics during drying and vapour absorption.

Extension to arbitrary pore shapes and complex pore structures
The method described above can be extended, with minor modification, to the modelling of moisture storage and transfer in arbitrary pore shapes and three-dimensional pore structures. A first exploratory study has been performed using the van der Waals equation of state, Eq. (7). In Fig. 11, the density distribution in an exemplary pore structure at different times during drying is shown. The material was initially saturated, with a relative humidity = 0.3 imposed on the upper boundary.
In this concern, an issue on the determination of the mobility arises, since Eq. (13) applies strictly for parallel plates only. This issue can be tackled by defining the mobility as a function of the position, including in this way the effect of the geometry on the flow.
The main challenge is however due to the high computational effort required by the model. In fact, in order to achieve realistic values of the surface tension, interface widths of the order of one nanometre are required. It seems therefore ambitious to perform simulations in the micrometre range, determinant for the pore structure of several building materials. The use of adaptive mesh refinement, ensuring sufficiently fine meshing in the interface region and allowing coarser meshing in the homogeneous regions, might be helpful in this regard to reduce the computational cost.

Conclusions and outlook
In this study, the applicability of the phase field method for simulation of isothermal moisture storage and transfer in porous materials has been investigated. It has been shown that the model gives results in agreement with the Young-Laplace and Kelvin equations in the steady state case. Moreover, the model is able to capture the dynamics of liquid uptake from a reservoir, as well as the vapour absorption/desorption process from/to an unsaturated atmosphere.
The main limitation of the method for a quantitative analysis on complex three-dimensional pore structures is due to its high computational cost. In fact, since interface widths of the order of one nanometre are required to obtain a realistic surface tension, simulations in the micrometre range seem very ambitious. Such numerical challenge is emphasized by the fact that the equation of state is very steep at low densities and the mobility varies of several orders of magnitude over the interface.
In order to tackle this computational challenge, the use of adaptive mesh refinement shall be considered in future research. In addition, a more accurate equation of state has to be employed in order to reach a better agreement with both measured interface width and surface tension.