A Unified Algorithm for the Young–Laplace Method Applied to Porous Media

Abstract

Young–Laplace equation-based algorithms for simulating capillary-driven flow have been used in porous media research for more than two decades. However, the lack of a uniform mathematical description hinders a wider application of these algorithms, as well as impede their comparison. After conducting a detailed review of the most important publications in the area, we propose a unified algorithm. This resulting framework is capable of handling four distinct physical situations: drainage and imbibition with either compressible or incompressible displacement fluid. Additionally, there is no restriction regarding the geometry or the initial fluid distribution used. The proposed algorithm can simulate variable or mixed wettability, even for imbibition, which has not already been described in the literature. Finally, at the end of the manuscript, we provide an efficient open-source C +  + code for the proposed unified algorithm and also many examples of its use.

Appendices

Appendix 1. Mathematical Morphology Applied to YLM

Figure 9 illustrates how the region occupied by the invading fluid is found using Mathematical Morphology (MM) operations. In the YLM, it is considered that the image is exclusively formed by the union between the solid and pore ($P$) related pixels and that the MM operations will not expand the image beyond its original dimensions.

Fig. 9

Examples of the fundamental MM operations. a The porous region P is shown in beige. b and c The beige region represents the set of points resulting from the operations of erosion and dilation, respectively. d The final result

In MM, the opening operation ($\circ$) is defined using the combination of two fundamental operations: erosion and dilation. In the erosion operation, denoted by $\ominus$, a structuring element with a disk shape $D$ is applied over the region $P$ (Fig. 9a), resulting in region $E=P\ominus D$, as shown in Fig. 9b. Region $E$ comprises all the points in $P$, where the structuring element can be centered without overlapping the solid region. Notice that although no segment can overlap over the solid region, it is allowed for part of the structuring element to be outside of the boundary image.

Figure 9c exemplifies the dilation operation of $E$ by a structuring element $D$ (denoted by $E\oplus D$), resulting in region $O$ with all the points intersected by the disk $D$ as its center moves around the $E$ region. It is important to notice that the result includes only points originally belonging to $P$; otherwise, the final image would be expanded.

Lastly, Fig. 9c is also the region $O$ obtained directly as a result of the opening operation of $P$ by the structuring element $D$ ($O=P\circ D$), which comprises all the points of $D$, as this element moves around $P$ without extending over the solid region. As the case shown in Fig. 9b, the structuring element is allowed to lay over regions outside the image. Moreover, as in Fig. 9c, the resulting region $O$ contains only points belonging initially to $P$. Figure 9d shows the final results after relabeling the porous region.

As one can realize from Fig. 9, the erosion of $P$ by $D$ followed by the dilation by $D$ is identical to the opening of $P$ by $D$, i.e.,

$$O=P\;\circ\; D=E\oplus D=(P\ominus D)\oplus D$$

(4)

Actually, it is possible to prove that it represents a general result, which is of great practical importance because it allows performing the opening operation using the combination of dilation and erosion operations. Furthermore, dilation can be obtained by applying erosion to the negative of an image. The negative is retrieved by swapping the image foreground and background, which in our case are porous and solid regions, respectively. Denoting the negative by an overbar, one can show that

$$E\oplus D=\overline{\left(\overline E\ominus D\right)}$$

(5)

The right-hand side operation is shown in Fig. 10a–d, which gives the same result already displayed in Fig. 9, which depicts the operation of the left-hand side. Using Eq. (5), the opening operation in Eq. (4) is written solely in terms of erosions

$$O=\overline{\left(\overline E\ominus D\right)}=\overline{\left[\overline{\left(P\ominus D\right)}\ominus D\right]}$$

(6)

Fig. 10

The opening operation of P is performed by two erosion operations. This is possible because the negative of the dilation is equal to the erosion of the negative. a Results from the negative of the E region in Fig. 9b. b and c The beige region represents the set of points resulting from the operations of erosion and the negative of the N region, respectively. d The final result

As previously mentioned, to obtain an efficient algorithm, the erosion operations must be performed with the aid of the Euclidian Distance Transform (EDT). The result of the EDT is a map with the distance of each pore pixel to the closest solid pixel, and it can be computed using a very efficient and parallel algorithm proposed by Saito and Toriwaki [37].

If the structuring element is a circle of radius $R$, then the erosion operation is given by the set of pixels whose distance to the solid is larger than $R$ or, equivalently, ${{\text{EDT}}}_{ij}^{2}>{R}^{2}$ where ${{\text{EDT}}}_{ij}^{2}$ indicates the squared distance of the pixel located at the image position $ij$ to the closest solid pixel [39]. Using the squared distance is advantageous since it allows for the computation to be handled solely by integers. Figure 11 depicts how the ${\text{EDT}}$ can be used to compute the opening operation of the porous region by a circle. The procedure starts from the ${\text{EDT}}$ result, shown in Fig. 11a. As explained, the values are the squared distance to the background (solid). By taking the negative of this result, i.e., setting it as background (Fig. 11b), the ${\text{EDT}}$ is once again applied from which a new erosion is performed. The negative result of this last operation (Fig. 11c) gives the opening of the porous region by the circle (Fig. 11d).

Fig. 11

Opening operation of the porous region by a circle of 7 pixels in diameter (i.e., ${R}^{2}={3.5}^{2}=12.25$) using the EDT (see Fig. 2). a EDT values used to erode the porous region. b EDT values used to erode the negative image of a. c The negative image of b is equivalent to the opening operation of the porous region directly. d Final image after opening operation of the porous region by a circle of 7 pixels in diameter

Figure 12 illustrates the bottom-up drainage through a narrow channel. The non-wetting fluid source is represented by the dark blue square (Fig. 12a). By the opening operation, and consequently, by the Young–Laplace equation, the upper region could be occupied by the non-wetting fluid; however, the circle cannot cross the constriction for this radius (Fig. 12b). From the morphological standpoint, this question is addressed by discarding the upper region since it is not connected to the source.

Fig. 12

a A bottom-up drainage through a narrow channel with the wetting fluid in red, and the non-wetting fluid source represented by the dark blue square. b The opening operation identifies two regions (beige color) that could be occupied by c the invading non-wetting fluid (light blue), but the upper region is not connected to the fluid source (dark blue pixel), and hence cannot be filled in the current pressure step

The next step is the identification of the connection between the regions can be performed by using the algorithm of He et al. [40]. This technique assigns a label to each unconnected group of pixels, allowing us to ignore the pixels with labels different from that of the sources (Fig. 12c). This methodology already considers the possibility of multiple sources.

By using the MM notation discussed above, the YLM can be expressed as [4]:

$$s\left(r\right)=\frac{\sum {F}_{{{\text{I}}}_{{\text{v}}}}(P\;\circ \;D)}{\sum P}$$

(7)

where $s(r)$ is the saturation of the invading fluid for a given circle radius ($r$, and $D\equiv D(r)$), and the summation symbol represents the operation of counting the pixels of each region. As for the operator ${F}_{{I}_{{\text{v}}}}$, it denotes the process of filtering the pixels connected with the invading fluid source.

As described in Sect. 2, imbibition can be obtained in a similar way to the drainage processes. But actually, some implementation details distinguish these two processes. These differences were noticed by the earlier authors of the method [12,13,14], and they were discussed in Sect. 3.

Appendix 2. Main Improvements for Drainage Simulation with the YLM

Since the majority of the applications of the YLM are limited to drainage simulations, this process received greater attention from researchers over the years compared to imbibition. Consequently, important improvements were proposed in the literature for the drainage process, and here, we describe two significant changes to the original procedure summarized by Eq. (7).

Hilpert and Miller [4] identified an error in Eq. (7) which led to an overestimation of the saturation value in some situations. As displayed in Fig. 13a, for a drainage process of the wetting fluid (red) from the bottom to the top in a selected geometry with two large pore regions separated by a short channel (like an ink bottle geometry), the opening operation results in the two regions being connected (Fig. 13c). Nevertheless, during a real invasion, the non-wetting fluid (blue) would not reach the upper section for the value of pressure related to the circle radius indicated in Fig. 13b, since it cannot cross the channel. This would only be possible if the entire circle could cross the channel.

Fig. 13

a A drainage process of the wetting fluid (red) from the bottom to the top in a selected geometry with two large pore regions separated by a short channel. Common to both cases, the erosion operations are performed in b. c, e, and g The original YLM drainage algorithm as proposed by Hazlett [12] and also by Magnani [13] and Magnani et al. [14]. d, f, and h The correction of Hilpert and Miller [4]. The non-wetting invading fluid (blue) source is represented by a dark blue square

To avoid this error, Hilpert and Miller [4] proposed a simple correction (Fig. 13d) that consists of identifying the regions connected to the fluid source between the erosion (Fig. 13b) and the dilation (Fig. 13f) processes instead of after the dilation (Fig. 13e). Figure 13g, h shows the difference between the original YLM drainage and Hilpert’s correction YLM drainage, respectively. Considering this correction, Eq. (7) is replaced by

$$s\left(r\right)=\frac{\sum \left\{\left[{F}_{{{\text{I}}}_{{\text{v}}}}(P\ominus D)\right]\oplus D\right\}}{\sum P}$$

(8)

The drainage process described by Eq. (8) is limited to the simulation of the invasion of a fully non-wetting fluid, or equivalently, the value for the contact angle in Fig. 1a is always zero. Schulz et al. [5] presented a small modification to this algorithm, allowing different values to be set. The concept used by the authors is depicted in Fig. 14. It shows that if a circle with a different radius is used (${R}^{\prime}$ for erosion and $R$ for dilation); then, a final interface with a radius of curvature ${R}^{\prime}$ and contact angle $\theta$, given by

$${R}^{\prime}=R\;{\text{cos}}\theta$$

(9)

is achieved. Including this improvement, Eq. (8) is modified to Eq. (2) in which the structuring element ${D}^{\prime}$ is associated with a circle with a radius ${R}^{\prime}$. Even though this generalization enables more realistic applications of the YLM, a fact that it is not yet widely known by the research community, as recent works still point out the inability to change the contact angle [10, 41].

Fig. 14

The sequence of MM operations proposed by Schulz et al. [5] to enable drainage simulations with different contact angles. The disks $D^{\prime}$ and $D$ are associated with ${R}^{\prime}$ for erosion and $R$ for dilation, respectively

One main advantage of this approach is concerning to its implementation easiness. Since just the radius of the erosion needs to be modified, it is easily incorporated into an existing YLM code. Moreover, as the erosion is computed pixel by pixel using the ${\text{EDT}}$ result, a case of mixed wettability could be readily performed by assigning an individual contact angle to each pixel [5].

On the other hand, since the porous media surface is, in general, curved, the contact angle will depend on this curvature. Also, if the larger circle used in the dilation extends beyond the wall to a near pore, this region can be incorrectly assigned as occupied by the invading fluid. Although, for it to happen, this region must also have a connection with the fluid source. This situation is shown in Fig. 14, where the circle with a radius ${R}^{\prime}$ extends to a bottom channel.

As mentioned by the authors, in petrophysics applications, these two drawbacks do not present a serious problem if $\theta \lesssim 40^\circ$ as the rock samples generally have low porosity and a large number of small pores.

As shown in Sect. 4, an appropriate membrane makes Schulz’s proposed modification suitable even for imbibition.