Droplet Formation, Growth and Detachment at the Interface of a Coupled Free-FLow–Porous Medium System: A New Model Development and Comparison

Abstract

Coupled free-flow–porous medium systems are of great importance in various natural and industrial applications. Modeling of such systems is always challenging, especially when droplets form at the interface between the two domains. We propose a new concept to take droplet formation, growth and detachment at the interface into account. In this concept, we use pore-network modeling to describe the porous medium and the Navier–Stokes equations for the free-flow domain. New coupling conditions are developed which include droplet interactions with the free flow and the porous medium. Impacts of using different descriptions of the forces acting on the triple contact line and contact angle hysteresis on the predicted onset of the droplet detachment are examined. In addition, we compare the new approach with another model built using ANSYS Fluent based on the volume of fluid method. The results show that the new model is able to describe the droplet formation, growth and then detachment by the free flow. The proposed model provides a base for further developments to handle formation of multiple droplets at the interface between a free flow and a porous medium as well as to include the evaporation in future works.

1 Introduction

Formation of droplets at the interface of a coupled free-flow–porous medium system occurs in many engineering applications such as fuel cells, cooling systems, membrane emulsification and filtration, thermal insulation and air conditioning of buildings (Zhu et al 2007; Arai and Suidzu 2013; Glass et al 2001; Charcosset 2009; Rashidi et al 2018). Droplet formation tends to increase the complexity of the system through turning a simple interface, which only handles the exchange between the free flow and the porous medium, to a complex interface, which is also able to store mass and energy. Such a complex interface must be able to capture the droplet dynamics and their impacts on the interaction between the free flow and the porous medium. Depending on the surrounding flow conditions, droplets might spread and merge, or be detached by the free flow (Baber et al 2016; Ackermann et al 2021).

In addition, having a clear description of the porous medium is a crucial requirement in modeling of such coupled systems. In macro-scale modeling of the fluid flow in porous media, the fluid and porous medium properties are defined by averaging the microscopic properties over a representative elementary volume (REV) and Darcy’s law describes the fluid flow (Zhang et al 2000). Pore-scale models, on the other hand, provide more details of the phenomena occurring in porous media. These models, as their name implies, give an understanding of the fluid movement in a pore using a detailed description of the fluid configurations as well as fluid–fluid and fluid–solid interactions in the pore space (Blunt 2017). Direct numerical simulations, Lattice Boltzmann method and pore-network models are examples of methods which can be used in pore-scale modeling. In comparison with other methods of modeling of pore-scale phenomena, pore-network models offer a low computational cost by eliminating the need for tracking the interfaces between the fluids and solid, while preserve a high level of accuracy (Joekar-Niasar and Hassanizadeh 2011).

Fig. 1

A schematic of a coupled free-flow–pore-network system (Weishaupt et al 2019a)

Mosthaf et al (2011) proposed a coupling concept for a compositional system of a single-phase free flow and a two-phase porous medium on REV-scale. They used a two-domain approach, such that the Navier–Stokes equations and Darcy’s law are employed to respectively describe fluid flow in the free-flow domain and the porous medium. The domains are coupled through a simple interface which transfers mass, energy and momentum, but the interface cannot store them. The numerical implementation and performance of such a coupling concept for an evaporation process from a porous medium is presented in Baber et al (2012). Due to the importance of the pore-scale effects in porous media for interface-driven coupled systems, Weishaupt et al (2019a) used a pore-network model to describe the porous medium in a single-phase system and employed a fully monolithic coupling approach at the interface between the free flow and the porous medium. Figure 1 shows a schematic of a coupled free-flow–pore-network systems used in Weishaupt et al (2019a). Weishaupt et al (2019a) also showed the high efficiency of the coupled model in terms of computational cost and accuracy in comparison with a reference case, which used the Navier–Stokes equations to describe the two sub-domains in a direct numerical simulation. In Weishaupt et al (2019b), the coupling conditions between the pore network and the free flow are improved and the model reduction in the free-flow domain is implemented. Weishaupt (2020) also analyzed the impact of pore geometry and pore-network heterogeneity on the numerical behavior of a non-isothermal compositional coupled system consisting of a single-phase free flow–two-phase pore network for modeling of evaporation from porous media.

To benefit from local accuracy of a pore-network model in modeling of large porous domains, a multi-scale approach could be employed. In such an approach, a porous medium could be divided into sub-domains according to the transport process intensity. Then, a pore-network model could be used to describe the parts of the porous medium, which have major impacts on the system (e.g., near-interface regions), coupled with the remaining parts of the porous medium modeled using an REV-scale approach (Weishaupt 2020).

To capture the impacts of droplet dynamics at the interface between a free flow and a porous medium on REV-scale, Baber et al (2016) extended the coupling concept for a simple interface to include the mass and momentum of the droplet by developing a REV-drop concept based on the concept presented by Mosthaf et al (2011) and Baber et al (2012). Baber et al (2016) used an area-weighted averaging of the coupling conditions of the uncovered part of the interface and the part covered with the droplet. The free flow and the porous medium interact directly at the uncovered part, while the interaction between the two domains occurs through the droplet at the covered part of the interface. Thus, they described the droplet dynamics in an averaged manner and do not resolve the droplet. The droplet volume is computed by introducing an additional degree of freedom at the interface using the mortar method. However, their concept requires prior knowledge about the number and location of the droplets forming at the interface and is not able to capture the impact of the droplets on the free flow. Qin et al (2012) investigated water flooding in gas channels of polymer electrolyte fuel cells by applying Darcy’s law in the porous medium (gas diffusion layer) and in the free-flow domain (gas channel). In their approach, it is assumed that the gas channel behaves as a porous medium with porosity of one; hence, Darcy’s law is applicable to describe the two-phase flow therein. Ackermann et al (2021) developed a three-domain approach which describes the interface as a lower-dimensional domain to include the droplet influence on a coupled free-flow–porous medium system. In this approach, the interface is physically a thin layer of the free-flow domain, which is in contact with the upper layer of the porous medium. The impacts of the droplet behavior are included in the interface description on the REV-scale. To account for film flow on the surface of the porous medium, they employed a similar approach as in Qin et al (2012), which used Darcy’s law in the free-flow channel. Accordingly, Ackermann et al (2021) developed a liquid saturation–relative permeability relationship for the free-flow channel based on the amount of liquid in the channel. Although their approach can describe the droplet formation and growth in the interface domain, it uses averaging to develop the coupling conditions between the domains and also cannot capture the impact of the droplets on the free-flow field. Using the volume of fluid method, Niblett et al (2020) analyzed the impact of pore morphology in porous media on free-flow–porous media interactions and the performance of fuel cells. Their results stressed the importance of the interface conditions in optimization of fuel cell efficiency. They also found that the location of droplet formation on the surface of the porous medium plays a crucial role in water management process.

These previous studies show that we need to deepen our understanding of the impact of multiple droplets emerging at the interface in a coupled free-flow–porous medium system. Considering the complexity of such systems, development of a model which provides sufficient details in predicting the system behavior with a reasonable computational time is desirable.

The aim of the present study is to develop a model being able to handle the formation, growth and detachment of multiple droplets at the interface, and integrate them in the coupling conditions between the free flow and the porous medium. Since in this coupled system, interface processes dominate the behavior, we focus on understanding the coupling conditions applying to mass and momentum transport between the two flow domains through the interface and in the presence of droplets. Considering the huge impact of the phenomena occurring in the porous medium on the interface processes, a pore-network model is used to describe the fluid flow in the porous medium and to take the pore-scale processes into account. The Navier-stokes equations are used to describe the free flow. It should be noted that the models developed in this study are applicable to hydrophobic interfaces and porous media and can also be used for hydrophilic media with minor adjustments.

The paper is structured as follows: In Sect. 2, we address the key parameters determining droplet behavior on a solid surface (Sect. 2.1), the models used in each domain and for the droplet (Sects. 2.2–2.4). We explain interactions of a droplet with a free flow and a porous medium and elaborate on the development of the coupling conditions including the droplet impact as well as the modeling of the droplet behaviors in Sect. 2.5. For comparison purposes, we also use a model based on the volume of fluid method implemented in ANSYS Fluent (refer to Appendix A for more detail). The results of droplet dynamics simulation as well as comparison with the ANSYS Fluent model is discussed in Sect. 3. The focus of Sect. 4 is on the features of the developed model and the next steps to improve it. Finally, we give conclusions and outlook in Sect. 5.

2 Model Description

In this section, firstly, we discuss how the wettability and interfacial tension between phases impact a droplet on a solid surface. In the second part, we describe the model and assumptions that we use to describe a droplet at an interface. After that, we explain the fundamentals of the pore-network model used in the porous medium. Then, the model used to describe the free-flow domain is discussed. The last part focuses on the droplet dynamics occurring at the interface and the coupling conditions.

All models which are introduced in this section are implemented in DuMu$^\text {x}$ , an open-source simulation toolbox for transport in porous media (Koch et al 2021).

2.1 Droplet on a Solid Surface: Interfacial Tension and Wettability

In a system composed of two or more phases, an interface separates two immiscible fluids (e.g., liquid and gas) or a fluid and a solid. Formation of an interface between two phases is a result of the difference between the inter-molecular forces in the phases. As an example, if we have a water droplet surrounded by air, the net inter-molecular force acting on a water molecule due to its neighboring molecules is zero. Surrounded partially by other water molecules, a water molecule at the interface between the phases, however, experiences inter-molecular forces different from a water molecule inside the droplet. At the interface, the net inter-molecular force on the water molecule points inwards of the droplet. The combined effect of the radial components of inter-molecular forces across the entire interface surface is to make the surface contract, thereby increasing the pressure on the concave side of the surface. The difference between the inter-molecular forces of the phases can be reflected by surface or inter-facial tension, which is a universal property of the interface (Brutin 2015). The surface tension acts along the interface to minimize the free energy by decreasing the interface area between the phases. For instance, a similar strength of inter-molecular forces in adjacent phases results in a smaller interfacial tension between them. This might cause a weaker interface and ultimately partial or complete miscibility of the fluids.

In a system consisting of a sessile droplet of liquid, l, on a solid surface, s, surrounded by a gas, g, three phases exist. The surface tension works between each pair of phases. In such a system, a three-phase contact line forms on which, in equilibrium state, a balance between the three surface tensions, $\gamma$, needs to be fulfilled. This is described by Young’s equation (Young 1832),

$$\begin{aligned} \gamma _{lg}\cos \theta = \gamma _{gs} - \gamma _{ls}. \end{aligned}$$

(1)

In the above equation, $\theta$ denotes the equilibrium (static) contact angle of the liquid phase with the solid surface. The contact angle indicates the degree of wettability between a fluid and a solid. Considering liquid water, a surface is called hydrophilic when $\theta < 90^{\circ }$ and hydrophobic when $\theta > 90^{\circ }$. Figure 2a illustrates a droplet on a hydrophobic surface in static condition and the surface tensions acting on the triple contact line.

Fig. 2

a A static sessile droplet: the surface tensions acting on the three-phase contact line on a hydrophobic surface in equilibrium (static) condition and b the advancing and receding contact angles of a sessile droplet surrounded by a fluid flowing from left to right

The contact angle used in Eq. (1) is the equilibrium (static) contact angle. The equilibrium (static) contact angle is observed all along the triple contact line if the solid surface is homogeneous and no external forces act on the droplet. In real applications, however, various elements perturb the equilibrium assumption in the system, which leads to variation of the contact angle along the contact line. This variation is called hysteresis. The main origins of contact angle hysteresis are surface roughness, heterogeneity of the solid surface and chemical interactions between fluids and the solid surface. In a case where external forces are present, when contact angle hysteresis is very little (e.g., an ultra hydrophobic surface), the droplet starts moving easily (Li et al 2007). The contact angle hysteresis allows a deformation of the droplet before it starts moving due to the external forces.

Figure 2b shows the advancing contact angle, $\theta _\text {adv}$, and the receding contact angle, $\theta _\text {rec}$, which determine the contact angle hysteresis in a system. The maximum advancing contact angle is reached once the contact line is about to start moving in favor of the advancing of the contact area. On the other hand, the minimum receding contact angle is the contact angle right before the contact line moves toward the decreasing of the contact area. The ultimate contact angle hysteresis is the difference between the maximum advancing and minimum receding contact angles. For instance, a sessile droplet on a solid surface in a field of gas flow can experience various contact angles along its contact line, with values between the advancing contact angle and the receding contact angle.

2.2 Droplet Description

Assuming a static condition, a homogeneous surface and a small droplet such that the effect of the gravity field on it is negligible, the droplet deformation can be neglected and a droplet forms as a part of a sphere having a circular contact area with the solid surface. Equations (2) and (3) describe a spherical droplet shown in Fig. 3, using the droplet radius of curvature, R, the contact radius of the droplet with the solid surface, $r_\text {CA}$, the height of the droplet, h, and the droplet contact angle, $\theta$. The capillary pressure, $p_c$, due to the curvature of the droplet surface and the surface tension between the phases, is described using Eq. (4) (Baber 2014).

Fig. 3

A sessile droplet and its descriptive parameters

$$\begin{aligned} R&= \frac{h}{1-\cos {\theta }} = \frac{r_\text {CA}}{\sin \theta }\, , \end{aligned}$$

(2)

$$\begin{aligned} V_\text {drop}&= \frac{\pi }{3}h^2(3R-h) = \frac{\pi }{3}r_\text {CA}^3\frac{(1-\cos (\theta ))^2(2+\cos (\theta ))}{(\sin (\theta ))^3} \,, \end{aligned}$$

(3)

$$\begin{aligned} p_c&= \frac{2{\gamma }_{lg}}{R} \frac{2(1-\cos \theta )+\cos \theta {(\sin \theta ) }^2}{(1-\cos \theta )^2(2+\cos \theta )}\,. \end{aligned}$$

(4)

2.3 Porous Medium: Pore-Network Model

In this section, we briefly discuss the pore-network model which is used to model the porous medium flow. For more details, we refer to Weishaupt (2020), which thoroughly discusses the different aspects of the pore-network model.

We model the porous medium using a pore-network model, which describes the porous medium as a network of pore bodies scattered in a solid bulk interconnected through tubes known as pore throats. Figure 4 shows two pore bodies, i and j, connected to each other by pore throat ij. In this model, the pore throats determine the conductivity behavior of the system and the pore bodies are responsible for the storage capacity.

The pore-network model uses the well-known Hagen-Poiseuille equation as the basic principle to describe the one-dimensional flow in each pore throat of the network. The Hagen-Poiseuille equation can be derived from the Navier–Stokes equation by assuming a one-dimensional, fully developed, stationary laminar flow in a single pore throat (Weishaupt 2020). Neglecting the impact of gravity, the general equation to calculate the flow rate of phase $\alpha$, $Q_{\alpha ,ij}$, in pore throat ij connecting pore bodies i and j is:

$$\begin{aligned} Q_{\alpha ,ij} = \frac{A_{\alpha ,ij} r_{\alpha ,ij}^2}{8\mu _\alpha l_{ij}}(p_{\alpha ,i} - p_{\alpha ,j}), \end{aligned}$$

(5)

where $p_\alpha$ is the pressure of phase $\alpha$ in the pore body, $\mu$ is the fluid viscosity, l is the pore throat length. $A_{\alpha }$ and $r_{\alpha }$ are the effective pore throat cross-sectional area and the effective pore throat radius for phase $\alpha$, respectively, which can be calculated as functions of the pore throat shape and the fluid phase presence in the pore throat.

Fig. 4

Pore bodies i and j are connected by pore throat ij

In a non-compositional flow, Eq. (6) describes the mass balance of phase $\alpha$ in each pore body, i.

$$\begin{aligned} V_i\frac{\partial (\rho _\alpha S_\alpha )_i}{\partial t} + \sum \limits _{j}(\rho _\alpha Q_\alpha )_{ij} - (Vq_\alpha )_i = 0. \end{aligned}$$

(6)

In the above equation, the storage term is the first term on the left-hand side in which V is the volume of the pore body and $S_\alpha$ indicates the saturation of phase $\alpha$, which in a single-phase system is equal to one. The second term, the advection term, is the sum of the fluid flow into/out of the pore body in interaction with the connected neighboring pore bodies. The last term on the left-hand side refers to the possible source/sink in the pore body.

In the pore-network model, the primary variables are assigned to the center of the pore body, where the balance equations need to be fulfilled.

2.4 Free-Flow Domain

The continuity equation, Eq. (7), is used to fulfill mass balance in a single phase free-flow domain,

$$\begin{aligned} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho {\mathbf{v}}) - q = 0, \end{aligned}$$

(7)

while the Navier–Stokes equations describe the fluid flow in the free-flow domain,

$$\begin{aligned} \frac{\partial (\rho {{\mathbf{v}}})}{\partial t} + \nabla \cdot \left( \rho {{\mathbf{v}}}\otimes {{\mathbf{v}}}\right) + \nabla \cdot \left( p {\varvec{l}} - \mu (\nabla {{\mathbf{v}}} + \nabla {{\mathbf{v}}}^T)\right) - \rho {\varvec{g}} = 0. \end{aligned}$$

(8)

In the above equations, $\rho$ is the fluid density, $\mathbf{v}$ is the velocity vector, q is the possible source/sink, p is the fluid pressure, $\mu$ is the fluid viscosity and $\varvec{g}$ is the gravity vector.

We employ a staggered-grid finite volume approach to discretize Eqs. (7) and (8) in space. In such an approach, the velocity degrees of freedom are located at the face of the primary cell, which is the center of the relevant secondary grid cell. The rest of the degrees of freedoms such as pressure, temperature (in case of non-isothermal flow), and mole/mass fraction (in case of compositional flow) are located at the center of the primary grid cell. Schneider et al (2020) provides a detailed description on the staggered-grid finite volume approach.

2.5 Interface Region

In a coupled system consisting of a free flow and a porous medium, the coupling conditions between the two domains at the interface determine the behavior of the system. In the absence of droplets on the surface of the porous medium, the free flow and the porous medium directly interact with each other. Formation of droplets causes complexity in the coupling conditions through forming two additional interfaces: one droplet–free-flow and another droplet–porous medium interface. Having a droplet at the interface, the free flow and the porous medium are not only directly connected with each other but also through the droplets. In the following, we discuss the droplet formation and growth on the surface of a porous medium, as well as the impact of such a process on the coupling conditions.

2.5.1 Droplet–Porous Medium Interactions

Mass balance A droplet on the surface of a porous medium stores the mass coming from the porous medium, which here is considered as a pore network. Equation (9) shows the volume change of the droplet due to the interaction with a single pore. In this equation, $\rho _l$ is the density of the droplet which is assumed to be the same as the porous medium liquid phase, $V_\text {drop}$ is the droplet volume, ${{\mathbf{v}}}_l$ is the velocity vector of the fluid in the pore, ${\varvec{n}}_\text {pore}$ is the unit normal vector and $A_\text {drop}^\text {pore}$ is the area of the droplet–pore interface.

$$\begin{aligned} \left( \frac{\text{ d }(\rho _l V_\text {drop})}{\text{ d } t}\right) _\text {drop}^\text {pore} = [\rho _l {{\mathbf{v}}}_l\cdot {\varvec{n}}]_\text {pore}A_\text {drop}^\text {pore}. \end{aligned}$$

(9)

Momentum balance At the interface between the droplet and the pore, as shown in Fig. 5a, continuity of forces in the direction normal to the drop surface is described as:

$$\begin{aligned} \varvec{F}_\text {pore} + \varvec{F}_\text {drop}=0, \end{aligned}$$

(10)

where $\varvec{F}_\text {pore}$ and $\varvec{F}_\text {drop}$ are the forces exerted on the droplet–pore interface due to the pore and the droplet pressure, respectively. Substitution of the forces by the multiplication of the pressure and an infinitesimal area of the droplet–pore interface, $\text {d}a_\text {drop}^\text {pore}$, results in the mechanical coupling condition based on the droplet and pore pressures.

$$\begin{aligned}{}[(p\varvec{n})_\text {pore} + (p\varvec{n})_\text {drop}]\cdot \varvec{n}_\text {drop}\ \text {d} a_\text {drop}^\text {pore} = 0. \end{aligned}$$

(11)

Rearrangement of Eq. (11) gives:

$$\begin{aligned} p_\text {drop} = p_\text {pore}, \end{aligned}$$

(12)

which shows that at every point on the droplet–pore interface, the droplet pressure is equal to the pore pressure.

Fig. 5

a Forces acting on the drop–pore interface and b Forces acting on the drop–free-flow interface

2.5.2 Free-Flow-Droplet Interactions

Mass balance A droplet and the surrounding free flow exchange mass through evaporation from the surface of the droplet.

$$\begin{aligned} \left( \frac{\text{ d }(\rho _lV_\text {drop})}{\text{ d }t}\right) _\text {drop}^\text {ff} = -f_\text {evp}A_\text {drop}^\text {ff}, \end{aligned}$$

(13)

where $f_\text {evp}$ is the evaporative flux from the droplet surface to the free flow and $A_\text {drop}^\text {ff}$ is the area of the droplet–free-flow interface.

Momentum balance on the droplet surface Continuity of forces in the normal direction at the drop surface, Fig. 5b, at each infinitesimal area, $\text {d}a$, needs to be fulfilled:

$$\begin{aligned} {\varvec{F}}_\text {ff} + {\varvec{F}}_{p_c} + {\varvec{F}}_\text {drop}=0, \end{aligned}$$

(14)

where $\varvec{F}_\text {ff}$ is the total free-flow force acting on each point of the droplet–free-flow interface. Considering an infinitesimal area of the droplet surface, $\text {d}a_\text {drop}^\text {ff}$, the forces composing the free-flow force are as follows:

• Inertial force: ${\varvec{F}}_{\rho \mathbf{v}^2} = ((\rho {{\mathbf{v}}} \otimes {{\mathbf{v}}} ) {\varvec{n}})_\text {ff}\cdot {\varvec{n}}_\text {drop}\ \text{ d }a_\text {drop}^\text {ff}$,

• Pressure force: ${\varvec{F}}_p = (p {\varvec{n}})_\text {ff}\cdot {\varvec{n}}_\text {drop}\ \text{ d }a_\text {drop}^\text {ff}$,

• Shear force: ${\varvec{F}}_{\tau _\text {ff}} = ((-\varvec{\tau} ) {\varvec{n}})_\text {ff} \cdot {\varvec{n}}_\text {drop}\ \text{ d }a_\text {drop}^\text {ff}$, where $\varvec{\tau}$ is the shear stress tensor.

The capillary force, $\varvec{F}_{p_c}$, is defined as:

• Capillary force: ${\varvec{F}}_{p_c} = (p_c {\varvec{n}})_c\cdot {\varvec{n}}_\text {drop}\ \text{ d }a_\text {drop}^\text {ff}$, where the subscript c indicates that the capillary force stems from the curvature of the droplet surface.

The droplet force exerted by the fluid inside the droplet on the droplet surface is as follows:

• Droplet force: ${\varvec{F}}_\text {drop} = (p {\varvec n})_\text {drop}\cdot \varvec{n}_\text {drop}\ \text{ d }a_\text {drop}^\text {ff}$.

Thus, Eq. (14) can be rewritten as:

$$\begin{aligned} {\varvec{F}}_{\rho \mathbf{v}^2} + {\varvec{F}}_p + {\varvec{F}}_{\tau _\text {ff}} + {\varvec{F}}_{p_c} + {\varvec{F}}_\text {drop}=0. \end{aligned}$$

(15)

Substituting the definition of each force, we derive Eq. (16), which shows the mechanical coupling condition at the interface between the droplet and the free-flow domain.

$$\begin{aligned}{}[((\rho {{\mathbf{v}}} \otimes {{\mathbf{v}}} ) {\varvec{n}})_\text {ff} + (p {\varvec{n}})_\text {ff} + ((-\varvec{\tau} ) {\varvec{n}})_\text {ff} + (p_c {\varvec{n}})_c + (p {\varvec{n}})_\text {drop}]\cdot {\varvec{n}}_\text {drop}\ \text{ d }a_\text {drop}^\text {ff} =0. \end{aligned}$$

(16)

Droplet impact on the free-flow field To calculate the free-flow forces acting on the droplet surface, we need an estimation of the droplet influence on the free-flow field. Such an estimation should be applicable in a free-flow region discretized in space using the staggered grid scheme as well as being able to include the impacts of multiple droplets in the system. To do so, we take the following approach:

• 1. The droplet properties are calculated using Eqs. (2) and (3) according to the droplet–porous medium interaction (described by Eq. (9)). In case of evaporation of the droplet in the free-flow region, its impact on the droplet volume should be included according to Eq. (13).

• 2. Having the droplet physical properties and position, we identify the grid faces in the free-flow region invaded by the droplet. To determine how much of a grid face is occupied by the droplet, the dimensionless variable $\beta$ is defined as:

  $$\begin{aligned} \beta = \frac{\displaystyle \text {The area of the face occupied by the droplet}}{\displaystyle \text {The total area of the face}}. \end{aligned}$$

  (17)

  If $\beta = 1$, the grid face is completely inside the droplet, otherwise if $0<\beta <1$, the face is partially occupied by the droplet and it is recognized as an interface grid face.

• 3. The interface grid cells are the cells which contain at least one interface grid face. Figure 6 illustrates the interface grid cells in red and the cells inside the droplet in blue.

• 4. If a grid face is completely inside the droplet, we set the velocity on the face to zero, ${\mathbf{v}} = 0$.

• 5. The coupled system is solved, including the droplet impact on the free-flow field explained in the previous step.

• 6. We calculate the free-flow shear and inertial forces on the interface grid faces and the pressure force on the interface grid cells.

Fig. 6

The blue grid cells are completely inside the droplet and the red ones are on the interface

2.5.3 Coupling Concept: Free Flow–Droplet–Porous Medium

To get the final picture, we need to combine what we discussed before in Sects. 2.5.1 and 2.5.2 for the coupling conditions.

Mass balance The total mass balance of the coupled system can be derived through adding up the droplet mass exchange with the porous medium and the free flow.

$$\begin{aligned} \frac{d(\rho _lV_\text {drop})}{dt} = [\rho _l {{\mathbf{v}}}_l\cdot {\varvec{n}}]_\text {pore}A_\text {drop}^\text {pore} - f_\text {evp}A_\text {drop}^\text {ff}. \end{aligned}$$

(18)

Momentum balance To derive the momentum coupling condition, we substitute the drop pressure in Eq. (15) with the pore pressure according to Eq. (11), which leads to:

$$\begin{aligned}{}[((\rho {{\mathbf{v}}} \otimes {{\mathbf{v}}} ) {\varvec{n}})_\text {ff}+(p {\varvec{n}})_\text {ff}+((- \varvec{\tau} ) {\varvec{n}})_\text {ff} + (p_c {\varvec{n}})_c + (p_\text {pore} {\varvec{n}})_\text {drop}]\cdot {\varvec{n}}_\text {drop}\ \text{ d }a_\text {drop}^\text {ff} =0. \end{aligned}$$

(19)

2.5.4 Droplet Formation and Growth

Both the free flow and the porous medium flow can affect the formation of droplets on the surface of the porous medium. In the case of droplet formation as a result of the porous medium liquid breakthrough, it is the porous medium that plays the role of a supplier for the droplet and acts in favor of the droplet growth. In fact, the droplet stores the mass coming from the porous medium until it is detached as a result of the forces that the free flow exerts on it. Evaporation from the surface of the droplet in the free-flow region, on the contrary, causes shrinkage of the droplet, as expressed in Eq. (18). In this study, we ignore evaporation from the surface of the droplet to the free flow for the sake of simplicity and describe the droplet volume as a function of the porous medium out-flow, as shown in Eq. (9).

To model how the droplet grows, we divide the droplet growth process into two periods:

• 1. Initially, when the droplet starts growing, the droplet contact angle is less than the advancing contact angle of the hydrophobic surface. In this stage, the droplet contact line sticks to the pore and the droplet contact angle grows until it reaches the surface advancing contact angle. In other words, a constant contact radius and increasing contact angle mode describes the droplet growth.

• 2. The second stage begins once the droplet contact angle reaches the advancing contact angle of the surface. At this time, we switch to a constant contact angle approach which means that the droplet grows while its contact angle is constant, i.e., equal to the surface contact angle, and the droplet contact radius increases. Thus, the droplet contact area expands beyond the connected pore.

Fig. 7

Droplet formation and growth on a pore, where $\theta$ is the droplet contact angle before it reaches the surface advancing contact angle, $\theta _\text {adv}$, while the droplet contact radius is equal to the pore radius, $r_\text {pore}$. After the triple contact line expands and leaves the pore, the droplet contact radius is $r_\text {CA}$

It should be noted that in both stages, we assume that the droplet is part of a sphere (i.e., non-deformed droplet), whose radius can be calculated based on the droplet contact angle and contact radius, as described in Sect. 2.2. Figure 7 shows the stages of a droplet growth on the surface of a pore.

2.5.5 Droplet Detachment

The detachment of the droplet from the porous medium surface has a significant impact on the coupled system of free flow–porous medium. To predict the droplet detachment due to the free-flow influence, it is necessary to recognize and calculate the forces acting in favor and against the droplet detachment.

As discussed in Sect. 2.5.2, free flow around the droplet exerts forces on the droplet surface on every point of the droplet–free-flow interface. Integration of the forces over the droplet surface gives the total free-flow force which acts in favor of the droplet detachment. This force can be decomposed to the force in the flow direction (i), called drag force, and perpendicular to the flow direction (j) called lift force:

$$\begin{aligned} {\varvec{F}}_\text {total}^\text {ff} = {\varvec{F}}_\text {drag} + {\varvec{F}}_\text {lift}. \end{aligned}$$

(20)

The drag force is the main force that acts to detach the droplet from the solid surface, which can be calculated as:

$$\begin{aligned} {\varvec{F}}_\text {drag} = {\varvec{F}}_{\tau i} + {\varvec{F}}_{pi} + {\varvec{F}}_{pv^2i}, \end{aligned}$$

(21)

where

$$\begin{aligned}& {\varvec{F}}_{\tau i} = \int _{A_\text {drop}^\text {ff}}((- \varvec\tau ) {\varvec{n}})_\text {ff}\cdot {\varvec{n}}_i~\text{ d } a_\text {drop}^\text {ff}\,, \end{aligned}$$

(22)

$$\begin{aligned}& {\varvec{F}}_{pi} = \int _{A_\text {drop}^\text {ff}}(p {\varvec{n}})_{\text {ff}}\cdot {\varvec{n}}_i~\text{ d } a_\text {drop}^\text {ff}\,, \end{aligned}$$

(23)

$$\begin{aligned}& {\varvec{F}}_{\rho \mathbf{v}^2i} = \int _{A_\text {drop}^\text {ff}}((\rho {{\mathbf{v}}}\otimes {{\mathbf{v}}}) {\varvec{n}})_\text {ff}\cdot {\varvec{n}}_i~\text{ d } a_\text {drop}^\text {ff}\,, \end{aligned}$$

(24)

where $\varvec{F}_{\tau i}, \varvec{F}_{pi}$ and $\varvec{F}_{\rho {\mathbf{v}}^2i}$ are the total shear, pressure and inertial force components in the flow direction. $\varvec{n}_i$ is the unit vector in the flow direction (i).

The impacts of the drag force on the droplet before detachment are the deformation in the shape of the droplet as well as the contact angle hysteresis along the triple contact line between the free-flow fluid, the droplet fluid and the solid. The droplet contact angle hysteresis creates a force; we call it “hysteresis force”, which is acting on the contact line in the opposite direction to the drag force and tries to prevent contact line movement. The hysteresis force is influenced by the droplet deformation, the droplet liquid–free-flow gas interfacial tension and the wettability state of the solid surface of the porous medium. From the force balance, it follows that this force is equal to the drag force when the droplet is still attached to the surface, i.e., no movement of the triple contact line. The maximum hysteresis force occurs when the droplet experiences the ultimate difference between the advancing and receding contact angles, occurring just before the droplet detachment. Figure 8 shows a sessile drop under influence of the free-flow drag force acting in the free-flow direction and the hysteresis force acting on the triple contact line.

Fig. 8

Free-flow drag force and hysteresis force acting on a droplet

Several expressions to calculate the hysteresis force have been proposed (Chen et al 2005; Kumbur et al 2006; Cho et al 2012). Kumbur et al (2006) assumed a linear variation of the contact angle along the triple contact line and derived Eq. (25) by integrating the forces on the contact line to describe the maximum hysteresis force, ${\varvec{F}}_\text {hyst}^\text {max}$:

$$\begin{aligned} \begin{aligned} {\varvec{F}}_\text {hyst}^\text {max} = \gamma _{lg}\pi R\sin \theta \left[ \frac{\sin (\Delta \theta ^\text {max} - \theta _\text {adv}) - \sin (\theta _\text {adv})}{\Delta \theta ^\text {max} - \pi } +\frac{\sin (\Delta \theta ^\text {max} - \theta _\text {adv}) - \sin (\theta _\text {adv})}{\Delta \theta ^\text {max} + \pi }\right] . \end{aligned} \end{aligned}$$

(25)

In Eq. (25), $\Delta \theta ^\text {max}$ is the maximum difference between the advancing and receding contact angles. Assuming two semi-circles with contact angles of receding and advancing contact angles, Chen et al (2005) proposed the following equation to calculate the maximum hysteresis force on the triple contact line.

$$\begin{aligned} {\varvec{F}}_\text {hyst}^\text {max} = 2\gamma _{lg}\pi R\sin \theta \sin \left(\frac{\theta _\text {adv} + \theta _\text {rec}}{2}\right) \sin \left( \frac{\Delta \theta ^\text {max}}{2}\right) . \end{aligned}$$

(26)

Cho et al (2012) simplified Eq. (26) by assuming $(\theta _\text {adv} + \theta _\text {rec})/{2} = \theta$ and proposed the following equation:

$$\begin{aligned} {\varvec{F}}_\text {hyst}^\text {max} = 2\gamma _{lg}\pi R\sin ^2\theta \sin \left( \frac{\Delta \theta ^\text {max}}{2}\right) . \end{aligned}$$

(27)

When the free-flow drag force exceeds the maximum hysteresis force, the triple contact line becomes unstable and starts to move in the direction of the drag force. Whether a droplet first slides for some time at the interface and then detaches or detaches rather immediately once the triple contact line becomes unstable, depends on different factors, e.g., surface wettability, fluids surface tension, droplet contact angle hysteresis, droplet growth rate, pore geometry and free-flow velocity (Li et al 2012; Wang et al 2011). Nevertheless, in the developed model here, we assume that the movement of the triple contact line due to the free-flow drag force is the detachment criterion, which determines when the droplet detaches from the solid surface. Such a criterion has shown to be a proper indicator of the droplet detachment (Kumbur et al 2006; Cho et al 2012).

3 Results

3.1 Simulation Setup

Figure 9 shows the setup in which air flows through a channel with dimension of $1.35\ {\text{mm}}\times 4.15\ {\text{mm}} \times 1\ {\text{mm}}$ (width $\times$ length $\times$ height). A vertical pore throat with circular cross section of $0.15\ {\text{mm}}$ radius and $0.5\ {\text{mm}}$ length is connected to the middle of the bottom wall of the channel. Air enters the channel at the inlet with a fully developed laminar velocity profile given by (Hartnett and Kostic 1989):

$$\begin{aligned} {\mathbf{v}}(y,z) = {\mathbf{v}}_\text {max} \cdot \left( 1-\left( \frac{y-a}{a}\right) ^{2}\right) \cdot \left( 1-\left( \frac{z-b}{b}\right) ^{2.3}\right) , \end{aligned}$$

(28)

where ${\mathbf{v}}_\text {max}$ is the maximum flow velocity at the center line, and a and b are half of the width and height of the channel, respectively. Atmospheric pressure is assigned as the Dirichlet boundary condition for pressure at the channel outlet. Initially, there is no flow in the channel and the air pressure is equal to the atmospheric pressure.

Water is injected into the inlet of the pore throat at a constant mass flow rate of $10^{-5} \ {kg}/{s}$. A droplet forms and grows in the channel at the outlet of the pore throat until it is detached due to the air flow in the channel.

For comparison purposes, we use the above-described setup for simulations using a model based on the volume of fluid method, which is implemented in ANSYS Fluent. ANSYS Fluent R2019 is chosen for this study due to its established and stable implementation of a multiphase solver and the applicability on high performance computing clusters for fast calculations (Michalkowski et al 2022). Appendix A explains the basics of volume of fluid method used in ANSYS Fluent model in more details.

Fig. 9

The simulation setup

3.2 Grid Resolution in the Free-Flow Domain

To find the suitable grid resolution in the free-flow domain such that it preserves accuracy but at the same time prevents high computational cost, we apply three different grid resolutions in a case with maximum free-flow velocity of $20 \ {m/s}$, mean velocity of $9.3 \ {m/s}$, at the inlet of the free-flow channel. For the fine grid resolution, we use 110400 grid cells to discretize the domain. We use 16875 cells for the medium grid resolution and 6000 cells to generate the course grids in the free-flow domain. However, in all cases, a refinement algorithm is employed to have smaller grids around the droplet and larger grids in areas far from the droplet. For instance, in the fine grid resolution, we have 48000 cells in the area where droplet presence is probable, whereas there are 2700 grid cells in the coarse grid resolution. Figure 10 illustrates the aforementioned grid resolutions. In this figure, the area in the middle of the channel with smaller grid cells is where the droplet is more likely to invade, and therefore, a finer grid is used. It should be noted that the fluid flows in x-direction in the free-flow channel.

Fig. 10

The three grid configurations used to analyze the impact of grid resolution: a the fine grid, b the medium grid and c the coarse grid resolution

Table 1 compares the predicted droplet properties once the detachment occurs for the three grid resolutions in the free-flow channel. It can be seen that using the coarse grid in the channel results in more deviation from the predicted values by the simulation using the fine grid. The simulation results using the medium and fine grid resolutions show less than one percent difference in prediction of the detachment height and a bit more than two percent difference for the detachment volume.

Table 1 Predicted droplet height and volume at the detachment time using three different grid resolutions in the free-flow channel with $v_\text {max} = 20 {m/s}$. $\Delta$ is the difference between the predicted value using the fine grid resolution and the other resolutions

Since the grid resolution in the free-flow domain mainly affects the forces acting on the droplet, we compare the variation of the free-flow drag force during the droplet growth predicted by the three aforementioned grid resolutions in Fig. 11. As can be seen, the drag force increases almost steadily and free of fluctuations when the fine grid is used. Whereas by coarsening the grid, the fluctuation increases, although the trend remains the same. This happens because having a fine grid resolution provides a better estimation of the free-flow grid faces invaded by the droplet in each time step. In fact, since in our approach the invasion of the grid face happens once the droplet occupies the whole grid face, the impact of this sudden invasion on the free-flow field is larger when the grid is coarser. It is also noticeable that the level of fluctuation becomes stronger over time due to the fact that the droplet invades and blocks the grid cells, which have a stronger influence on the free-flow velocity field.

Fig. 11

The impact of the grid resolution on the free-flow drag force in the free-flow channel with a maximum flow velocity of 20 m/s

Taking all these into account, it is reasonable to use the medium grid resolution for the remaining simulations, as it provides less computational cost while gives almost the same results as the fine grid. It is worth to note that the ANSYS Fluent model uses more than $2 \times 10^6$ grid cells.

3.3 Droplet Dynamics

Figure 12 shows how droplet properties and forces vary over time during a one formation-growth-detachment period of one droplet at the interface. We use the contact angles reported by Theodorakakos et al (2006) for carbon cloth. The equilibrium contact angle of water with the solid surface is $145^\circ$ in both channel and throat. The maximum advancing contact is $150^\circ$ and the minimum receding contact angle is $90^\circ$. The results shown in Fig. 12 belong to the case that the maximum velocity of free flow at the inlet of the channel is $15 \ {m/s}$, i.e., mean velocity of $6.98 \ {m/s}$. Figure 12a shows that the volume of the droplet increases with a constant rate over time due to the constant injection rate into the throat. Figure 12a also shows that the variation rate of the radius of curvature and the height of the droplet decreases by time, i.e., by the growth of the droplet. We can see in the Fig. 12b that the contact radius of the droplet stays constant and the contact angle increases over a time frame at the beginning of the droplet growth, followed by a constant contact angle and growing contact radius period. These two stages of the droplet growth are discussed in Sect. 2.5.4. If we look at Fig. 12c, the drag force increases almost monotonically with the droplet volume. The variation of the drag force becomes more extreme over time, although the rate at which the droplet height changes is decreasing. That means that when the droplet is larger, a change in the droplet height affects the free-flow field stronger than when the droplet is small. The hysteresis force peaks locally when the droplet contact angle reaches $90 ^\circ$. After that, it dips and reaches a local minimum when the droplet contact angle reaches the ultimate contact angle of the surface. By expansion of the droplet contact area, the hysteresis force recovers and starts an upward trend.

Fig. 12

Droplet dynamics: a variation of the droplet height, radius of curvature and volume over time, b variation of the droplet contact angle and contact radius during the droplet growth and c the drag and hysteresis forces variation until the droplet detaches

3.4 Droplet Detachment

In this section, first, we compare the results predicted by the developed model for the droplet detachment with the ANSYS Fluent model. Then, the separation line obtained from the developed model is compared with the line predicted by analytical/empirical approaches. As mentioned before, the model introduced in this study is implemented in DuMu$^\text {x}$ (Koch et al 2021) and is referred to as the “DuMu$^\text {x}$ model” in the current section.

3.4.1 Comparison with the ANSYS Fluent Model

Figure 13 compares the predicted droplet detachment volume and height versus the mean free-flow velocity for the DuMu$^\text {x}$ and ANSYS Fluent model. In the DuMu$^\text {x}$ model, we apply different descriptions of the hysteresis force, provided by Eqs. (25)–(27), to predict the detachment of the droplet and compare the results. It should be noted that in this comparison, the highest value used as the maximum free-flow velocity is $30 \ {m/s}$, the mean velocity is $13.95\ {m/s}$, while the lowest value is determined as the minimum velocity which detaches the droplet before reaching the upper wall of the channel. The static contact angle of the surface is $145 ^\circ$ and the advancing and receding contact angles used in the DuMu$^\text {x}$ model are $150 ^\circ$ and $90 ^\circ$, respectively. It can be seen in Fig. 13a that the predicted detachment volumes by the ANSYS Fluent model and the DuMu$^\text {x}$ model are closer at higher free-flow velocities. By decreasing the velocity, however, the difference between the results of the two models for the droplet detachment volume increases. The DuMu$^\text {x}$ model and the ANSYS Fluent model both predict almost the same minimum free-flow velocity by which the droplet can be detached by the free-flow. However, the detachment volume predicted by the ANSYS Fluent model is considerably larger than the value predicted by the DuMu$^\text {x}$ model.

To get a better understanding, Fig. 13b compares the height of the droplet when it detaches, which are predicted by the DuMu$^\text {x}$ and ANSYS Fluent model. It shows that although the different versions of the DuMu$^\text {x}$ model and the ANSYS Fluent model predict a similar trend, the difference of the predicted values using Eq. (25) with the other results is significant, especially at high free-flow velocities. Although such an observation may seem opposite to what we see in Fig. 13a, where the graphs converge by increasing velocity, it can be explained by looking back to Fig. 12a, where the variation of the droplet volume and height is compared. As explained in Sect. 3.3, a small change in the droplet volume, when the droplet is small, can drastically affect the droplet height. At a high free-flow velocity, the droplet detaches when it is small. Thus, although the predicted values for the droplet detachment volumes differ slightly, the corresponding difference in the predicted detachment heights is noticeable. On the other hand, the predicted values by all models for the least free-flow velocity at which the droplet can be detached are close. Such behavior might be related to the deformation of the droplet which is included in the ANSYS Fluent model, whereas is not taken into account by the DuMu$^\text {x}$ model. In fact, by lowering the free-flow velocity, the droplet can grow larger before the detachment. Therefore, the droplet experiences more deformation during its life at the interface, which enables it to gain more volume before it touches the upper wall of the channel.

Fig. 13

Comparison between the simulation results of the DuMu$^\text {x}$ model using various formulation for the hysteresis force and the ANSYS Fluent for a the droplet detachment volume and b the droplet detachment height

To have a deeper look at the impact of the deformation, we compare the detachment values predicted by the ANSYS Fluent model and the corresponding value computed by assuming no deformation in the droplet (see Fig. 14). For example, in Fig. 14a, the detachment volumes reported for the undeformed droplet are for a droplet with the same height as the simulated droplet by the ANSYS Fluent model, if there were no deformation. That means that the detachment height predicted by ANSYS Fluent model is used in Eqs. (2) and (3) to calculate the droplet volume. Similarly, for the detachment height of the undeformed droplet, it would be the height of a droplet with the same detachment volume as predicted by the ANSYS Fluent model, if the deformation was not included. In other words, the detachment volume predicted by ANSYS Fluent model is used to back calculate the detachment height by using Eqs. (2) and (3). According to Fig. 14a, the difference between the detachment volume for the deformed droplet predicted by ANSYS Fluent model, and for the undeformed droplet increases by decreasing the free-flow velocity. The comparison of detachment height by Fig. 14b also shows a similar behavior. That is to say, the impact of the droplet deformation on the detachment becomes more significant when the free-flow velocity is lower. As mentioned before, when the droplet becomes larger, the change in its height has more impact on the drag force that it experiences. Therefor, the change in the height of a large droplet due to the deformation affects the detachment process more significantly than a small droplet.

One of the main differences between the DuMu$^\text {x}$ model and ANSYS Fluent model is that the ANSYS Fluent model does not take the contact angle hysteresis into account. That means that the two models apply different criteria for droplet detachment. In the DuMu$^\text {x}$ model, the droplet detaches when the free-flow drag force becomes larger than the hysteresis force on the contact line and consequently the droplet contact line starts to move in the drag force direction, whereas in the ANSYS Fluent model, since the contact angle hysteresis is not included, the triple contact line can move freely over the surface while the droplet is still attached to the pore. That is to say, in the ANSYS Fluent model, the surface tension force tries to hold the liquid phase together while maintaining the minimum surface area rather than keeping the droplet on the surface. That means that the droplet detachment occurs once the free-flow drag force overcomes the surface tension force and phase rapture occurs. To take the impact of contact angle hysteresis in volume of fluid method into account, for instance, Theodorakakos et al (2006) took an approach to recalculate the contact angle by considering the droplet deformation due to a free flow. In another study, Fang et al (2008) used a transient modeling approach, in which the contact angle is updated locally in each time step for each cell at the triple contact line during the simulation. Since including the contact angle hysteresis in the volume of fluid method is not in the scope of the present work, we refer to Theodorakakos et al (2006) and Fang et al (2008) for more details.

To show how the maximum contact angle hysteresis impacts the detachment prediction by the DuMu$^\text {x}$ model, three different values for the maximum contact angle hysteresis are used in the simulations and the results are compared with the ANSYS Fluent simulation in Fig. 15. As expected, a droplet on the surface of the porous medium with higher ultimate contact angle hysteresis, can remain longer attached to the surface and withstand higher free-flow velocities. It can also be seen that a smaller free-flow velocity is needed to detach a droplet with lower contact angle hysteresis before it touches the upper wall of the channel. The results confirm that the contact angle hysteresis as a reflection of the solid surface properties plays a crucial role in the droplet detachment.

Fig. 14

Comparison between the predicted value by ANSYS Fluent model and the undeformed estimation of: a the droplet detachment volume and b the droplet detachment height

Fig. 15

Impact of the maximum contact angle hysteresis on the detachment prediction by the DuMu$^\text {x}$ model using Eq. (26) to describe the contact angle hysteresis: a drop detachment volume and b drop detachment height

3.4.2 Comparison with Analytical and Semi-analytical Derivations for the Droplet Detachment

There are several analytical and semi-analytical approaches describing the detachment behavior of liquid water droplets from hydrophobic surfaces in the literature (Chen et al 2005; Kumbur et al 2006; Cho et al 2012; Hao and Cheng 2010). Chen et al (2005) developed an analytical relation to predict the detachment of the droplet using a two-dimensional (2D) description of the force balance around a droplet on a solid surface, influenced by the surrounding gas flow. They used Eq. (26) to compute the hysteresis force. The 2D approximation represents an infinitely wide channel (parallel plates) but also an infinitely wide droplet (deformed cylindrical shape) such that shear forces at the sides of the droplet are not included. Except using Eq. (25), Kumbur et al (2006) used a similar approach to Chen et al (2005) to predict the detachment of the droplet. Cho et al (2012) extended the concept presented by Chen et al (2005) by computing the drag force employing a drag coefficient derived from numerical simulations. Figure 16 compares the separation line generated by the analytical/empirical approaches, the ANSYS Fluent model and the DuMu$^\text {x}$ model. In this figure, the dimensionless drop height is the ratio of the drop height to the channel height. The static contact angle and contact angle hysteresis used in this section are $145^\circ$ and $60^\circ$, respectively.

Fig. 16

Comparison of separation lines predicted by analytical/empirical approaches, the ANSYS Fluent and the DuMu$^\text {x}$ model using various description of the hysteresis force

The separation line generated by Cho et al (2012) derivation shows the best agreement with the DuMu$^\text {x}$ and ANSYS Fluent model. However, it does not include the possible contact of the droplet with the top wall of the channel, such that it predicts a dimensionless drop height larger than one for low free-flow velocities. The separation line generated by Chen et al (2005) derivation has the same issue at low velocities. Although the separation line predicted by Kumbur et al (2006) approach tends to a dimensionless drop height of one at low velocities, it suffers from lack of accuracy in prediction of the droplet detachment. The disability of the analytical approaches to generate a reliable separation line might stem from their two-dimensional derivation, which Cho et al (2012) tried to improve by using numerical simulations in a three-dimensional channel. Another point is that the DuMu$^\text {x}$ and ANSYS Fluent models show results for a channel width of $W_\text {channel} = 1.35 \ {\text{mm}}$ and height $H_\text {channel} = 1 \ {\text{mm}}$ (aspect ratio $\approx 0.74$), which is common in modern proton-exchange membrane fuel cells. Other studies, however, used channels with higher aspect ratios for fitting and validating the data.

3.5 Free Flow–Droplet–Pore Network

To show the application of the model in modeling of formation, growth and detachment of multiple droplets at the interface, we use a setup shown in Fig. 17a. In this setup, a two-dimensional channel with dimension of $10 \ {\text{mm}} \times 1 \ {\text{mm}}$ (length $\times$ height) is used as the free-flow domain. A pore network as the porous medium is connected to the bottom wall of the channel. The pore network consists of 142 pore bodies connected by 190 pore throats to each other. The radii of the pore bodies vary between 0.063 to 0.019 mm. The pore throats have radii from 0.04 to 0.014 mm. The free-flow channel is initially filled with air and the pore network is initially fully saturated with water. Air flows to the channel with maximum velocity of $15\ {\text{m}/\text{s}}$ in a fully developed laminar profile and the outlet of the channel is exposed to the atmospheric pressure. Water is injected to the inlet pores of the network with the rate of $5 \times 10^{-7}\ {\text{kg}/\text{s}}$. The droplets form and grow onto pore bodies with various sizes at the interface. Figure 17b shows a snapshot of such a process, where the multiple droplets formed at the interface influence the free-flow velocity field.

Fig. 17

Simulation of formation and growth of multiple droplets at the interface of a coupled free-flow–pore-network system: a the simulation setup, b formation of multiple droplets at the interface and their impacts on the velocity field of the free flow

4 Discussion

In this section, we discuss the features which distinguish the model introduced in this study, from the other developed models up to this time.

Employing a pore network to model the porous medium enables us to capture pore-scale phenomena. In addition, it makes it possible to develop the coupling conditions between the free flow and the porous medium without using average quantities in the porous medium. Furthermore, in a system consists of a free flow and a pore network, the exact locations of the pores in the porous medium and at the interface are known. Thus, the coupling conditions can be applied to the pores located at the interface, which interact directly with the free-flow domain. Accordingly, the exact position of the droplets emerging at the interface can be located, which is beneficial in modifying the coupling conditions by taking the droplet impact into account. In fact, the pores at the interface containing the same fluid phase directly interact with each other, while in the presence of the droplet they interact via the droplet. Thus, averaging the coupling conditions to include the droplet, as seen in Baber et al (2016), is not necessary.

No need of tracking the phase interface in the porous medium and using a new simplified approach to model the droplet growth in the free-flow domain reduce the computational costs of the new developed model. This model enables us to pursue our aims of including the formation, growth and detachment of multiple droplets at the interface between a free flow and a porous medium and developing new coupling conditions by taking the droplet impact into account. However, it focuses on the droplet as long as it is attached to the surface of the porous medium (i.e., by taking this approach, the life of the droplet after its detachment cannot be monitored). The presented model takes also the impact of a growing droplet on the free-flow field into account, which is of great importance in computing the free-flow forces acting on the droplet to predict the droplet detachment.

As discussed before, in the introduced model, the droplet detaches as soon as the free-flow drag force defeats the hysteresis force acting on the triple contact line. In other words, the onset of moving of the triple contact line of the droplet due to the drag force is considered as the time when the droplet completely detaches from the surface. However, a droplet may also slide or roll at the interface. A positive point of using the hysteresis force is to include the interface physical properties such as the wetting state, roughness and chemical composition in the detachment prediction by reflecting them in the static contact angle and the maximum contact angle hysteresis specific to the surface. These quantities need to be measured via experiments.

The developed model is able to handle supplying liquid to a growing droplet by more than one pore at the interface. However, how to describe a growing droplet which spreads over another pore filled with gas is still an open topic.

5 Summery and Outlook

We have developed a new model to deal with the formation, growth and detachment of droplets at the interface between a coupled free-flow–porous medium system. Pore-network modeling is used as a tool to enable us to capture pore-scale phenomena occurring in porous media. The formation and growth of a droplet is described and a new approach is taken to include the impact of the growing droplet on the free-flow field. New coupling concepts between the free flow and the porous medium are developed, which includes storing mass and momentum in the droplet. Description of the forces acting in the system is given and accordingly the droplet detachment is predicted. We have also analyzed the impact of the contact angle hysteresis on the detachment predictions. In addition, we have used and compared several formulations for computing the hysteresis force, which are derived assuming different descriptions of the contact angle hysteresis on the triple contact line. The ANSYS Fluent has been employed to build a model based on the Volume of Fluid method and the simulation results were compared with the new developed model introduced in this study. Despite differences in the two models such as not including the impact of droplet deformation on the free-flow field in the developed model and the different detachment criteria as well as the required high grid resolution by the ANSYS Fluent model, the results show a good agreement for the prediction of droplet detachment. The relatively low number of grid cells required by the new approach to describe the droplet growth in the free-flow domain is one of its advantages, which provides a low computational cost. However, the developed model enables us to describe the droplet behavior as long as it is at the interface and the behavior of the droplet in the free-flow domain after the detachment is not our focus in this study.

To improve the new concept, we need to include the influence of the droplet deformation on the droplet shape and consequently on the free-flow field. We also want to extend the model further to be able to simulate a non-isothermal compositional coupled system which can describe the droplet evaporation and the energy transport between the droplet, the free flow and the porous medium as well.

We presented an application of the new model to describe formation of multiple droplets at the interface. In a follow-on work, we will examine how formation, growth and detachment of a droplet might be influenced by the neighboring droplets, when multiple droplets emerge at the interface of a coupled free-flow–porous medium system. The current study is the first step of our project and gives a foundation for further developments.

Appendix A ANSYS Fluent model based on Volume of Fluid method

1.1 A.1 ANSYS Fluent model description

Multi-phase flows can be numerically calculated using either Euler-Lagrange or Euler-Euler approaches. The Volume of Fluid (VoF) method is a type of Euler-Euler approach and is a surface tracking technique. It is applicable if there are two or more immiscible fluids, when the interface position and behavior are important. The method uses only one set of momentum equations, which is shared by all fluids and the volume fraction of each fluid in each computational cell is tracked. The VoF method is a physically derived method, which directly solves the governing equations at each mesh element. The details of the VoF method implemented in ANSYS Fluent are described in the Theory Guide (Fluent 2009).

In the following, a brief summary of the main influencing model parts to describe droplet formation, growth and detachment in a channel flow is given.

1.1.1 A.1.1 Interface interpolation

ANSYS Fluent uses a geometric reconstruction of the interface between the phases. If a cell is near the interface between two phases, a geometric reconstruction scheme is applied that represents the interface using a piecewise-linear approach. As shown in Fig. 18, it is assumed that the interface between two fluids is linear within each cell and uses this linear shape for calculation of the advection of the fluids through the cell faces.

Fig. 18

Piecewise-linear approximation of the fluid interface in ANSYS Fluent in all cells in contact with the interface. Left: "real" interface shape. Right: linear-approximated interface shape (Fluent 2009)

Therefore, the position of the linear reconstruction is calculated based on the local volume fraction in each cell and its derivatives relative to the center of each partially-filled cell (interface cell). The algorithm then calculates the advective transport of the fluids through each cell face using the computed interface representation and the normal and tangential velocity distributions on that face. Afterwards, the volume fraction in each cell is recalculated using the balance of fluxes, which are calculated during the previous step.

Whenever a cell is completely filled with one fluid and no interface needs to be tracked there, the ANSYS’s standard interpolation schemes are obtained to calculate the fluxes across the cell faces.

1.1.2 A.1.2 Volume fraction equation

The tracking of interfaces between two phases is achieved by solving a continuity equation of the volume fraction of one of the phases. For the gas phase g, the equation is given exemplary by

$$\begin{aligned} \frac{1}{\rho _g} \left[ \frac{\partial }{\partial t} \left( \alpha _g \rho _g\right) + \nabla \cdot \left( \alpha _g \rho _g \mathbf{v}_g \right) \right] = \frac{1}{\rho _g} \left[ q_g + \sum _{p=1}^{n-1} \left( \dot{m}_{pg} - \dot{m}_{gp}\right) \right] \,, \end{aligned}$$

(A1)

with the gas density $\rho _g$, the face value of the volume fraction $\alpha _g$, the velocity vector $\mathbf{v}_g$, a source term $q_g$, which is usually zero, and the mass transfer from phase p to g, $\dot{m}_{pg}$, and vice versa, $\dot{m}_{pg}$, for a system with n phases.

The equation holds analogously to the gas phase for all occurring phases.

1.1.3 A.1.3 Discretization

The ANSYS Fluent time discretization is implemented with an explicit finite-difference interpolation scheme, which is applied to the volume fraction values that were computed at the previous time step

$$\begin{aligned} \frac{\alpha _g^{n+1}\rho _g^{n+1} - \alpha _q^n\rho _g^n}{\Delta t}V + \sum _f \left( \rho _g U_f^n \alpha _{g,f}^n\right) = \left[ q_g + \sum _{p=1}^{n-1} \left( \dot{m}_{pg} - \dot{m}_{gp}\right) \right] V \, . \end{aligned}$$

(A2)

Here, $n+1$ denotes the new time step and n is the index of the previous time step. The index f gives the cell face where the volume flux U normal to the face and the volume fraction $\alpha$ are evaluated. The volume of the cell is given by V.

In ANSYS Fluent the face value of the volume fraction $\alpha _{g,f}$ is computed from a first- or second order upwind scheme.

1.1.4 A.1.4 Fluid properties

The fluid properties appearing in the transport equations are determined by the occurring phases in the considered control volume. In a two-phase system, as it occurs in the liquid droplet in gas flow channel setup, we represent the phases by the subscripts g and l for gas and liquid, respectively. Here, we track the volume fraction of the gas phase. The density in each cell is given by the following formulation

$$\begin{aligned} \rho = \alpha _g \rho _g + \left( 1-\alpha _g\right) \rho _l\,. \end{aligned}$$

(A3)

All other fluid properties (e.g., viscosity $\mu$) are computed in this manner.

1.1.5 A.1.5 Momentum equation

An advantage of the VoF method is that a single momentum equation is solved throughout the whole domain. The resulting velocity field is shared among the phases. This is achieved with the dependency on the volume fractions of all phases and through the properties $\rho$ and $\mu$. The resulting momentum equation is

$$\begin{aligned} \frac{\partial}{\partial t} \left( \rho {\mathbf{v}}\right) + \nabla \cdot \left( \rho {\mathbf{v}} {\mathbf{v}}\right) = -\nabla p + \nabla \cdot \left[ \mu \left( \nabla {\mathbf{v}} + \nabla {\mathbf{v}}^T \right) \right] + \rho {\varvec{g}} + {\varvec{F}}_\text {ext} + {\varvec{F}}_{vol}\,. \end{aligned}$$

(A4)

Where p is the shared-field pressure, ${\varvec{g}}$ the gravity vector and ${\varvec{F}}_\text {ext}$ an external force vector. ${\varvec{F}}_{vol}$ is a volume force vector.

It is important to mention that a limitation of the shared-fields approximations occurs in case of large velocity differences locally between the phases. The accuracy of the velocities, which are computed near the interface, can be adversely affected.

1.1.6 A.1.6 Surface tension

As described in Sect. 2.1, the surface tension is a force, acting only at the surface and is required to maintain equilibrium in such instances. In ANSYS Fluent, the implemented surface tension model is the continuum surface force (CSF) model proposed by Brackbill et al (1992). The model adds the surface tension to the VoF calculation resulting in a source term in the momentum Eq. (A4). The CSF model includes a formulation where the surface curvature is computed from local gradients in the surface normal at the interface. Therefore, the surface normal ${n}$ is given by the gradient of the volume fraction of the phase i as:

$$\begin{aligned} {\varvec{n}} = \frac{\nabla \alpha _i}{\nabla \alpha _i}\,. \end{aligned}$$

(A5)

The curvature $\kappa$ is defined in terms of the divergence of the normal vector (Brackbill et al 1992)

$$\begin{aligned} \kappa = \nabla \cdot {\varvec{n}}\,. \end{aligned}$$

(A6)

The surface tension can be expressed in terms of the pressure jump across the surface. The force at the surface can be written as a volume force using the divergence theorem. This volume force is then used as the source term, which is added in the momentum equation (Eq. (A4)):

$$\begin{aligned} {\varvec{F}}_{vol} = \sum _{\text {pairs } ij, i<j} \sigma _{ij} \frac{\alpha _i \rho _i \kappa _j \nabla \alpha _j + \alpha _j\rho _j \kappa _i \nabla \alpha _i}{1/2 \left( \rho _i + \rho _j \right) }\,. \end{aligned}$$

This expression allows a smooth superposition of forces near cells, where more than two phases are present. Since, in this study, we only consider two fluid phases, the equation simplifies to:

$$\begin{aligned} {\varvec{F}}_{vol} = \sigma _{ij}\frac{\rho \kappa _i \nabla \alpha _i}{1/2 \left( \rho _i + \rho _j \right) }\,, \end{aligned}$$

(A7)

with the volume averaged density $\rho$ (Eq. (A3)). It can be seen that the surface tension source term is proportional to the average density in the cell.

1.1.7 A.1.7 Wall adhesion

To specify the equilibrium contact angle at the three-phase contact line (gas, liquid, solid), an option to specify a wall adhesion angle in conjunction with the surface tension model is available in ANSYS Fluent. The model is presented in detail in Brackbill et al (1992). The contact angle that the fluid is assumed to form with the wall is used to adjust the surface normal in cells near the wall, rather than imposing the boundary condition at the wall itself. This so-called dynamic boundary condition results in the adjustment of the curvature of the surface near the wall to fulfill the applied contact angle. (see Fig. 19).

The surface normal at the live cell next to the wall is given with the contact angle $\theta _\text {static}$ as:

$$\begin{aligned} {\varvec{n}} = {\varvec{n}}_l \cos {\theta _\text {static}} + {\varvec{t}}_l\sin {\theta _\text {static}}\,, \end{aligned}$$

(A8)

with the unit vectors normal and tangential to the wall ${n}_l$ and ${t}_l$, respectively. The combination of the contact angle with the normally calculated surface normal at the live cell next to the wall determine the local curvature of the surface, and this curvature is then used to adjust the body force term in the surface tension calculation.

When the droplet is in contact with the hydrophobic surface of the bottom channel wall, a constant contact angle is applied and a dynamic contact angle results from the position of the interface in the adjacent cells.

Fig. 19

Contact angle and resulting interface between the gas and liquid phases. In ANSYS Fluent a constant (static) contact angle $\theta _\text {static}$ is applied at the solid walls. The deformation of the droplet results from local stresses (Fluent 2009)

1.2 Limitations and assumptions of the VoF model

No contact angle hysteresis Surfaces are assumed to be perfectly smooth such that no contact angle hysteresis occurs. For the fluid-fluid-solid interface (three-phase contact line) an equilibrium contact angle is applied. However, the local forces on the droplet surface result in a drop deformation. Even though, the three-phase contact line is not pinned to the surface by contact angle hysteresis. The droplet is hold in its position by the filling throat such that the deformation and detachment of the droplet can be observed.

Fine Meshing The accuracy and numeric stability highly depends on the mesh used for the simulation. A fine mesh is required for adequate results of the multi-phase simulation.

Idealized geometry The considered geometry includes sharp edges. At these corners, the contact angle is not defined and the behavior of the fluid-fluid interface highly depends on the time step size and discretization. Therefore, no pressure inlet boundary condition can be applied for the liquid phase. However, in the considered setup with a mass flow inlet boundary condition, and having the main interest on the behavior of the droplet on the flat surface, this phenomenon is assumed to have a minor influence.

No phase change In the considered setup, no phase change is implemented. Effects of evaporation and condensation are assumed to be negligible.

Not conservative across the interface Using a VoF method, mass/volume loss occurs while tracking the whole mass on the system. This results from the interface interpolation and volume fraction equation (Eq. (A1)). The volume fraction equation is of hyperbolic type and the solution using an upwind scheme results in minor errors in the $\alpha$-field which causes an imbalance of mass/volume in the enclosed system. However, employing stabilization methods, this dissipation in $\alpha$ can be reduced to a minimum in a system with moderate flow velocities, as it is considered in this study. In addition, errors in the interface approximation can be reduced using a sufficiently fine mesh.