CFD modelling of the effect of capillary pressure on retention behaviour of water menisci at inter-particle contacts

This paper presents a Computational Fluid Dynamics (CFD) model on the effect of capillary pressure on the retention behaviour of a granular material. The model proposes an unprecedented CFD insight into the onset of liquid menisci at the inter-particles contact under varying hydraulic conditions. The present work models the material grains as smooth spherical particles that define a porous network filled by two interstitial fluids: air and silicon oil. The numerical model has been subsequently validated against experimental measurements of the degree of saturation at different capillary pressures taken by Dullien et al. [F.A. Dullien, C. Zarcone, I.F. MacDonald, A. Collins, R.D. Bochard. J. Colloid Interface Sci. 127, 2 (1989)] in a system of smooth glass beads flooded with silicon oil. Results from the numerical simulations confirm the good capability of the model to reproduce the experimental retention behaviour of the granular material. Finally, the present paper laid the basis for future CFD studies on the effect of various factors (e.g. hydraulic hysteresis, surface roughness and/or grain shape) on the capillary pressure acting at the interparticle contact.


Introduction
The capillary actions that raise at the interfacial equilibrium (i.e. meniscus) between vapour and liquid, play a fundamental role in several physical phenomena: flow of granular materials, friction between surfaces, adhesion between particles or particles and surfaces [1][2]. These phenomena are crucial in a very broad range of industrial and engineering applications, such as food industry, pharmaceutics, construction industry and agriculture to name a few.
The importance of capillary forces on the abovementioned physical phenomena has been first recognised starting in the '20s with the seminal models of Fisher [3] and Haines [4]. On the other hand, several experimental techniques have been developed to measure capillary pressure in porous materials, such as pressure plates and high-capacity tensiometers [5][6][7]. More recent techniques enabled instead an unprecedented visualisation of liquid menisci forming in unsaturated granular materials. An example of these experimental techniques is reported in Figure 1, which shows Scanning Electron Microscopy (SEM) of water menisci between contacting glass beads (scanned at LEMAS centre at the University of Leeds).
Increasingly more sophisticated models have been devised in parallel with the continuous development of experimental techniques to study capillary phenomena. Early retention models postulate a relationship between the degree of saturation (or equivalently the gravimetric * email: alejandro.lopez@deusto.es or volumetric water content) and the capillary pressure [8,9]. More recent models have incorporated the influence of hydraulic hysteresis and material deformation on the variation of the degree of saturation of granular and porous media [10][11][12]. More recently, some advanced capillary models have also considered the effect of the pore size distribution and microstructure on the retention behaviour of granular and porous materials [13,14]. The present work consists of the advanced modelling of capillary actions forming in an unsaturated granular material. The modelling approach is based on the Computational Fluid Dynamics (CFD) technique.
Results show that the proposed model captures well the experimental retention behaviour of a pack of smooth glass beads flooded with silicon oil at different pressures. Outcomes from the present study also lay the foundation for future work on the effect of particles surface roughness or grain shape on the intensity of capillary actions. This work will couple both experimental characterisation via micro Computed Tomography (CT) scans and CFD numerical modelling. Furthermore, the present paper can also trigger future research into the application of the CFD modelling approach to other types of granular materials and fluids (e.g. partially saturated soils) and/or to tackle other features of the retention behaviour of granular materials (e.g. material deformation and hydraulic hysteresis).

Numerical model
The fluid flow was solved using ANSYS Fluent CFD software tool. CFD and the Finite Volume Method [15][16][17] have been successfully used in simulations of microfluidic devices [18][19][20]. Given the quasi-static nature of the process, the values of the dynamic viscosity and density shown in Table 1 [1] and the characteristic length (inlet diameter of 10 µm), the fluid flow was determined to be in the laminar regime.
Turbulence is in fact rarely modelled due to the low Reynolds numbers obtained for these type of flows [19] and so is gravity [21,22], thus significantly simplifying the modelling. Simulating droplets and capillarity implies being able to computationally represent free surfaces between the oil and the air phases. In order to do this, the Volume of Fluid method (VOF) is used [23]. The VOF method is a common approach when simulating multiphase flow in microfluidic devices [24][25][26][27].
The VOF method calculates the interface between two fluids by tracking a scalar field variable F. This variable represents the volume fraction distribution of the fluids in the domain. When the value of F in a cell is between 0 and 1, the interface will pass through that cell. The equation defining the distribution of the volume fraction F is the passive transport equation, which is represented in Equation 1. This continuity equation assumes an initial distribution of the volume fraction and a velocity vector field.
This equation should be satisfied for both fluids and it defines the interface between different phases for values of 0 < F < 1. The domain considers the flow in the interstice between 16 particles enclosed in a twodimensional square-shaped domain. Four types of boundary conditions have been used in these simulations. An inlet pressure at the top boundaries and an outlet pressure at the bottom boundaries of the domain have been considered. The sides of the geometry have been set to behave like symmetry planes and all the circular boundaries are set to walls with no-slip condition and a contact angle of zero degrees, as experimentally observed by Dullien et al. [1] by means of a horizontal microscope. For the phase interaction between air and the silicon oil, a surface tension of 0.025 N/m has been considered. Initially, the domain is filled with oil and a different pressure is specified at the inlet for each case. The area-weighted average of silicon oil volume fraction is monitored, and the simulation is considered to be a steady state once the average has reached a constant value.

Experimental validation
Numerical simulations have been validated against experimental data by Dullien et al. [1], who determined retention curves in packs of smooth and rough etched glass beads by means of the porous diaphragm tensiometry technique. According to this technique, Dullien et al. [1] tested two distinct packs composed either by smooth glass beads with an average size of 545 µm or by rough glass beads with a larger average diameter of 631 µm. These packs of glass beads were placed inside a glass container and on top of a fritted disk in hydraulic connection with both an oil reservoir, which provided the wetting fluid during testing, and a levelling bulb partly filled with mercury, which controlled the level of oil injected within the glass beads packs. For this experimental campaign, Dullien et al. [1] used three different types of low viscosity oils (i.e. Soltrol 170, silicon oil and styrene monomer). The main properties of these three oils are summarised in Table 1. The numerical model proposed in the present paper focuses on the retention behaviour of smooth glass beads flooded with silicon oil, as it will be described in the following sections. Future work will instead focus on the CFD modelling of rough particles.
The area-weighted average of oil volume fraction (ratio of face area occupied by oil and total area of the domain) against computational iteration (Fig. 2) was plotted to check if the steady state is achieved. Both the degree of saturation and corresponding capillary pressure were computed at the final iteration.
In order to compare the 2D computational estimation of the degree of saturation with 3D experimental results, an image-based methodology is proposed to calculate the number of pixels representing void, liquid and solid phases as well as the average width of the liquid menisci. Figure 3 shows the process of converting Ansys output to a grayscale image and a binarized image. ImageJ software was used for this purpose. The 2D image-based estimation was converted to 3D by considering a diameter of 545 μm for the spherical particles and a distance between particles of 10 μm, thus resulting in a void ratio of 1.017. The distance of 10 μm was assumed based on experimental observation from Figure 1. The degree of saturation was estimated as the ratio between the computed volume of oil and the physical volume of voids. In particular, the volume of oil was calculated by assuming the length and height of 2D menisci is equivalent to the diameter and height of a disk like menisci in 3D. In our calculation, we also considered that each particle forming 4 menisci with neighbouring particles in 2D, while it is 6 particles in 3D regular packing. So, we increased the calculated volume by a factor of 1.5. Possible sources of errors include image resolutions and higher coordination number in experiments, and this may explain the discrepancy between computed and experimental retention behaviour. The capillary pressure was instead determined from the ANSYS outputs on both air and oil pressures. Figure 4 shows the computational formation of liquid menisci at particle contacts under a capillary pressure of 3.45 kPa. Table 2 instead summarises the computational values of the degree of saturation for the four simulations with varying capillary pressure. Numerical values of the ratio between capillary pressure and imposed pressure on the boundaries of the domain are also reported. in Table 2. The computational results are then plotted in Figure 5 using red asterisk markers to compare with experimental data. Inspection of Figure 5 indicates a good agreement between the experimental and the computed retention behaviour at both pendular and residual saturated states.

Concluding remarks
The paper has presented a CFD model that predicts the retention behaviour of a granular material subjected to various hydraulic loadings. In particular, the proposed model has been validated against the experimental data by Dullien et al. [1] obtained on a pack of smooth glass beads flooded with silicon oil at different levels of capillary pressures. The numerical 2D model consists in a system of 16 smooth and circular particles with a diameter of 545 µm and spaced with a minimum distance of 10 µm. During the simulations, the 16 particles pack was flooded with silicon oil at four distinct values of pressures applied at the domain boundaries, which resulted in different levels of capillary actions. The degree of saturation for each simulation was retrieved via image processing, i.e. by converting Ansys images in both grayscale and binarized images. The numerical values of the capillary pressure were instead directly obtained from ANSYS outputs in terms of both air and oil pressures. Results from the numerical simulations showed a good capability of the proposed CFD model to reproduce the retention behaviour of a granular material. Based on the present results, future work will be addressed to a) the CFD study of capillary actions rising in different granular materials under partially saturated conditions, b) the analysis of the effect of surface roughness and grain shape on the formation and evolution of liquid menisci at the inter-particle contact of granular materials and c) the study of other features of the retention behaviour of granular materials, such as the effect of hydraulic hysteresis and/or material deformation on the variation of the degree of saturation.