Understanding the capillary behavior using the extended reduced similar geometry method

Author's Accepted Manuscript

Understanding the capillary behavior using the extended reduced similar geometry method Bei Wei, Jian Hou, Kang Zhou, Bo Yu

www.elsevier.com/locate/ces

PII: DOI: Reference:

S0009-2509(14)00706-4 http://dx.doi.org/10.1016/j.ces.2014.11.052 CES12016

To appear in:

Chemical Engineering Science

Received date: 16 May 2014 Revised date: 5 October 2014 Accepted date: 21 November 2014 Cite this article as: Bei Wei, Jian Hou, Kang Zhou, Bo Yu, Understanding the capillary behavior using the extended reduced similar geometry method, Chemical Engineering Science, http://dx.doi.org/10.1016/j.ces.2014.11.052 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Understanding the capillary behavior using the extended reduced similar geometry method Bei Wei, Jian Hou, Kang Zhou, Bo Yu Author affiliations: Jian Hou: 1. School of Petroleum Engineering, China University of Petroleum, Qingdao, Shandong, 266580, China Bei Wei and Kang Zhou: School of Petroleum Engineering, China University of Petroleum, Qingdao, Shandong, 266580, China Bo Yu: School of Science, China University of Petroleum, Qingdao, Shandong, 266580, China

Correspondence information: (The two corresponding authors contributed the ideas from different aspects) Corresponding author1: Jian Hou Affiliation: School of Petroleum Engineering, China University of Petroleum Address: School of Petroleum Engineering, China University of Petroleum, No. 66 Changjiang West Road, Economic& Technological Development Zone, Qingdao, Shandong, 266580, China E-mail: [email protected] Tel: 86-546-8395660 Fax: 86-546-8395660 Corresponding author2: Kang Zhou Affiliation: School of Petroleum Engineering, China University of Petroleum Address: School of Petroleum Engineering, China University of Petroleum, No. 66 Changjiang West Road, Economic& Technological Development Zone, Qingdao, Shandong, 266580, China E-mail:[email protected] Tel: 86-18766422967 Abstract Pore scale modeling is used to understand the transport phenomena in porous media. Capillary behavior is an important property in pore scale research. We extend the reduced similar geometry (RSG) method to tangential polygons and simplify the results to concise equations. Based on the RSG method results, the effects of the shape similarity, the shape factor and the contact angle on the capillary behavior are investigated. A new and more accurate parameter is proposed to predict the threshold radius during primary drainage. This parameter is demonstrated to predict the threshold radius accurately for arbitrary polygon shapes and real rock shapes. It also performs better than the traditional predictor when predicting the capillary behavior for concave shapes. Keywords: pore scale modeling, capillary behavior, RSG method, threshold radius 1. Introduction The phenomena of flow through porous media are ubiquitous in natural and artificial materials. Some examples include groundwater flow, underground oil flow, the microcirculation of blood in animals, the translocation of water in plants, and gas migration in a packed bed. Usually, the flow 1

involves multiple fluids. Therefore, pore scale modeling of multi-phase flow is important in understanding the flow mechanism at the microscale. The shape of the pore space in porous media is irregular and the connectivity is complex. On the one hand, it is possible to directly simulate flow based on the pore-space image using numerical Stokes solvers [Al-Omari and Masad, 2004; Zaretskiy et al., 2010] or the lattice Boltzmann (LB) method [Iliev et al., 2013]. However, these approaches can be computationally expensive. As a ‘top-down’ approach, the Stokes solvers discretize nonlinear partial differential equations and perform an extensive mesh generation. The method must also handle boundary conditions appropriately in the complex porous media [Wolf-Gladrow, 2000; Kutay et al., 2006]. The LB method is a ‘bottom-up’ approach and can handle the complex geometry more efficiently using a bounceback boundary. However, the LB model requires a “stream and collision” step at every node and can also be computationally expensive [Pan et al., 2004].On the other hand, the pore space can be described as a network of pores connected by narrow regions called throats using an idealized geometry. Next, the Hagen-Poiseuille equation is used to model the flow in the pore network [Piri and Blunt, 2005; Nejad Ebrahimi et al., 2013]. This simplification improves the calculation speed and effectively reduces calculation time. Usually, only a simple system of linear algebraic equations in the pore network model must be solved. Every pore or throat has its own shape in the pore network and a few shapes have been used to characterize the pore cross-sections. Circles, arbitrary triangles and squares (C-T-S) are the general abstraction of pore cross-sections [Oren et al., 1998; Patzek, 2001; Blunt, 2001]. Higher-order polygons are also popular. Helland et al. [2008] proposed a regular n-cornered star shape based on the C-T-S characterization. Some shapes with curved sides have also been used. MAN et al. [2000] studied the electrical resistivity and capillary pressure characteristics using a grain boundary pore (GBP) shape. In this study, we focus on polygonal shapes. Once the pore shape is decided, the capillary behavior in the tubes can be determined. Particularly in the case of multi-phase flow, the capillary behavior is the basis of the capillary pressure curve and relative permeability curve [Joekar-Niasar et al., 2013; Das et al., 2006]. An important property of the capillary behavior is the value of the threshold pressure (or capillary entry pressure), i.e., the minimum pressure difference required for the non-wetting fluid to penetrate a capillary filled with a wetting fluid. Although a numerical algorithm can be used to study the threshold pressure (Nassehi et al., 2011), the simulation may be complicated and the results cannot usually be reduced to a concise formula. The MS-P method, which was originally proposed by Mayer and Stowe [1965] and further developed by Princen [1969a, 1969b, 1970], has been widely used to analytically estimate the threshold pressure. The method is based on an energy balance and takes the principle of free energy minimization into account. Moreover, it is valid for different wetting conditions [Ma et al., 1996; Blunt, 1997; Jia, 2007]. Based on the MS-P method, Oren et al. [1998] derived the threshold pressures of triangular tubes in different flow conditions and Lago et al. [2001] derived the threshold pressure in a capillary with a polygonal cross-section. The group led by Helland has worked on the threshold pressure from 2D rock images [Frette et al., 2010; Zhou et al., 2014]. However, both of these equations are complicated. Meanwhile, starting with the MS-P method, Mason and Morrow [1991] proposed the reduced similar geometry (RSG) method. The RSG method was used to study the capillary behavior of a perfectly wetting liquid in irregular triangular tubes. Jia et al. [2007] extended the RSG method to irregular triangular capillaries with arbitrary contact angles. The RSG method provides a concise 2

analytical expression of threshold pressure; therefore, it is simple to analyze the physical meaning. Additionally, the RSG results are the same as the results by Oren and Lago. We focus on a quasi-static primary drainage process controlled by capillary pressure in the following discussion, i.e., the flow rate is very low and the capillary number is small. Therefore, phase saturation changes slowly with time and is regarded as an equilibrium condition. The discussion of the quasi-static drainage provides a theoretical basis for the dynamic process. With respect to the dynamic process, usually only the viscous force is coupled in a dynamic pore network. While Das et al. and his co-authors have done some systemic experimental work and numerical analyses to study the dynamic effects on capillary pressure, and a dynamic coefficient is introduced to describe the dynamic effects [Mirzaei et al., 2007; Das et al., 2012; Hanspal et al., 2013]. Because the analytical expression for the threshold pressure in the previous works is complicated, it is difficult to find a good predictor for the threshold pressure. Lindquist [2006] studied three parameters using computational geometry theory and found that the hydraulic radius is an excellent predictor of the threshold radius (and hence the threshold pressure) during primary drainage, while the inscribed radius and the area-equivalent radius over-predict the threshold radius. In this study, we extend the RSG method to tangential polygons. Next, the threshold pressure equation derived from the RSG method is reformed. A concise, more accurate parameter with some physical meaning is developed to predict the threshold pressure. Arbitrary polygons and real rock shapes are used to verify the new predictor. 2. MS-P Method and the RSG Method 2.1 MS-P method When the non-wetting fluid invades the wetting fluid in a non-circular capillary, the wetting phase will be displaced from the central region, leaving some residual in the corners. The invading meniscus is referred to as the main terminal meniscus (MTM) while the menisci formed in the cross-sectional plane are referred to as arc menisci (AMs) [Mason and Morrow, 1991]. As the drainage continues, the curvature of the AMs increases gradually until it reaches a threshold drainage curvature, at which the total fraction of the wetting phase is lowered to a specific level. The average number of the mobile menisci will also be conserved during the drainage, such that the birth of new menisci is compensated by the death of other menisci (Panfilov et al. 2005). In the absence of gravity, the MS-P method assumes the curvature of the AMs is the same as that of the MTM. Fig. 1 shows drainage in a straight triangular tube.

Fig. 1 Drainage in a triangular capillary: (a) lateral view; the curvature of the MTM is Cd and (b) sectional view; the curvature of the AMs is equal to that of the MTM (Mason and Morrow, 1990).

The capillary pressure Pc is calculated as follows:

Pc = Pn − Pw , 3

(1)

where the subscripts w and n refer to the wetting and non-wetting phases, respectively. If the AMs are displaced a small distance dx, the energy balance equation is written as follows: Pc Aeff dx = ( Lnw γ nw + Lnsγ ns − Lnsγ ws )dx,

(2)

where the subscript s indicates the solid, Aeff is the effective area occupied by the non-wetting phase, Lns is the length of the solid wall in contact with the non-wetting phase, Lnw is the perimeter of the AMs, and γ is the interfacial tension. When the contact angle is θ, Young’s equation yields the following: γns − γws = γnw cos θ ,

(3)

and the energy balance equation can be simplified as follows: Pc L + Lns cos θr Leff = nw = . γnw Aeff Aeff

(4)

If the radius of the threshold drainage curvature is rd (threshold radius), the Young-Laplace equation yields the following: Pc =

γnw . rd

(5)

The capillary pressure and the threshold drainage radius are then obtained by solving Eq. (4) and Eq. (5); the wetting phase saturation is calculated based on the geometrical characteristics of the cross-section. 2.2 RSG method Mason and Morrow [1991] applied the RSG method first in their analysis of perfectly wetting triangular tubes. This method greatly simplifies the analysis using a geometric construction by reducing the triangle sides and removing the central portion occupied by the non-wetting phase. As shown in Fig. 2(a), the center of the triangle is occupied by the non-wetting phase. The non-wetting phase is then broken into pieces, as shown in Fig 2(b), and the AMs are circle arcs with the same radius. Finally, a new triangle, which is similar to the original, is constructed, after the central portion is removed, as shown in Fig. 2(c). The AMs form a complete circle in the reduced triangle in the perfect wetting situation.

Fig. 2 Sketch of the RSG method (Mason and Morrow,1990).

Jia et al. [2007] extended the RSG method to irregular triangular capillaries with arbitrary contact angles. Considering the contact angle, the condition for the existence of the wetting phase in a corner is as follows: ε π (6) θ+ < . 2 2 4

where θ is the contact angle and ε is the angle of the corner, as shown in Fig. 3.

Fig. 3 Transversal section of a corner filled with the wetting phase (Jia et al.,2007).

3. Extended RSG Method The main advantage of the RSG method is that it constructs a reduced shape similar to the original. To ensure that the reduced shape is similar to the original, we consider the tangential polygons in which there must be an inscribed circle. Based on the models above, we derive the threshold pressure equation as concisely as possible. Fig. 4 shows an n-sided polygon and the reduced geometry. The reduced shape is similar to the original one and the similarity ratio is k. Denote the area of the original shape as A, the perimeter as P, the area of the reduced shape as A’, and the perimeter as P’. They are then related as follows:  P ' = kP, A ' = k 2 A .   A ' = 0.5P ' rd cos θ

(7)

For each corner, we define a parameter δi as follows:

0, εi > π − 2θ δi =  , i = 1, 2...n − 1, n. 1, εi ≤ π − 2θ

(8)

This is a flag that indicates whether the wetting phase occupies the corner or not.

Fig. 4 Transversal section of an n-sided polygon and the reduced geometry; the blue sections indicate the wetting phase. 5

The MS-P method describes the case as follows:

 Aeff = A − A '+ mrd2 .   Leff = ( P − P ') cos θ + 2mrd n

(9)

n

2π − ∑ (1 − δi )(π − εi ) − ∑ 2δi θ

n n ε i =1 i =1 m = ∑ δi cos θ sin θ + π + ∑ (1 − δi ) cos θ 2 cot i 2π 2 i =1 i =1 
 

  term 2

term1

(10)

term3

where m is a parameter decided by the shape and the contact angle that indicates the completeness of the circle in the reduced shape. The parameter m is equal to π in the perfect wetting condition, indicating that there is a complete circle in the reduced shape. Eq. (9) is a key step in our model because it links Aeff and Leff through m; term1, term2 and term3 in Eq. 10 are represented by the green sections, the red sections and the yellow sections, respectively, in Fig. 4 to calculate the Aeff or Leff .

Using Eq. (4) to Eq. (10), the threshold radius is solved as follows:

rd =

P cos θ ± 4 Am P2 2( cos 2 θ − m) 4A

(11)

If the threshold radius is given by a valid root (the radius is smaller than that of the inscribed circle), then it can be reformed as follows: P cos θ m = + rd 2G G

(12)

where G is the shape factor G = A / P2 Eq. (12) not only considers the contact angles but also has a concise form similar to that found by Mason and Morrow. It can be further simplified as follows: 1 cos θ m / π cos θ m/π = + = + rd 2rh rarea rin rarea

(13)

where rh is the hydraulic radius rh = A / P , rarea is the area-equivalent radius rarea = A / π , rin is the radius of the inscribed circle and rin = 2rh for the tangential polygons. Eq. (12) can also be written as follows:

1 cos θ 4mG cos θ + 4mG = + = rd rin rin 2 rin where

(14)

4mG ranges from 0 to 1 and its value reaches 1 when the shape is a circle and the contact 6

angle is 0. Eq. (13) and Eq. (14) are important to understand the capillary behavior further. The threshold wetting phase saturation can also be expressed in terms of m as follows:

Sw = [ A '− mrd2 ] / A = (

2G cos 2 θ m − ) )2 ( 4G 2 G cos θ + 4Gm

(15)

4. Discussion The shapes with same inscribed radius and shape factor tend to have similar capillary behavior according to Eq. (14) and Eq. (15). Especially in the perfect wetting case, whether the shapes are triangles or higher-order polygons, the capillary behavior is the same. However, even though the shape factors are the same when the contact angle is not 0, the capillary behavior may vary to some extent. This is because shapes with the same shape factor are not unique; therefore, the value of m may be different. The effect of m is essentially the effect of the similarity of the shapes when the contact angle and the shape factor are fixed. In this section, the influence of the shape similarity, the shape factor and the contact angle on the capillary entry pressure and the threshold saturation is discussed. 4.1 The construction of quadrangles with the same inscribed circle Quadrangles with the same inscribed radius are constructed to study the effect of m. The quadrangles are constructed using the following method. A unit circle is set in the polar coordinate system with its center at (0, 0). A random point, outside the circle, in the first quadrant x1, with coordinates (r1, θ1), is selected. Similarly, a random point in the third quadrant x3 is selected with coordinates (r3, θ3). Subsequently, the tangent line of the circle through the selected point is drawn. There are four tangent lines. The tangent equations are obtained geometrically. The four tangent lines intersect to form a quadrangle, as shown in Fig. 5. The vertexes of the quadrangle are calculated using coordinate operations. The other geometric parameters of the polygon, such as the perimeter, the area and the angles are consequently calculated. With the quadrangles thus constructed, we can calculate the capillary behavior for different wetting conditions.

Fig. 5 The construction of quadrangles with the same inscribed circle.

4.2 The effect of m on the capillary behavior Define the normalized capillary entry pressure Pcnorm as follows: 7

Pcnorm =

rin r Pc = in = cos θ + 4mG γnw rd

(16)

Fig. 6 The effect of m on the capillary entry pressure.

The normalized capillary entry pressures for the quadrangles for different wetting conditions are shown in Fig. 6. The capillary entry pressure decreases as the contact angle increases for shapes with same shape factor and increases as the shape factor increases when the contact angle is fixed. Usually, shapes with the same shape factor have slightly different normalized capillary entry pressure profiles for a specific wetting condition. The difference is caused by m and it will increase with the contact angle and decrease with the shape factor. The maximum difference is at approximately 0.1 when the contact angle is 60°. Therefore, the difference can be ignored (0.1/ Pcmax =5%).

Fig. 7 (a) The effect of m on the threshold wetting saturation and (b) the value of m for different wetting conditions corresponding to (a); different colors are different values of m.

The threshold wetting phase saturations for the quadrangles in different wetting condition are shown in Fig. 7; the values of m are also indicated. Both the wetting phase saturation and m decrease as the contact angle increases for shapes with same shape factor and they decrease as the shape factor increases when the contact angle is fixed. However, the wetting phase saturations of the shapes with the same shape factor may vary in certain wetting conditions. This difference is also caused by m and it increases as the contact angle increases and decreases as the shape factor increases. However, the difference of the saturation is more pronounced than in the case of the capillary entry pressure when the contact angle is large. The results in Fig. 7(a) show that the difference is as high as 0.05 (0.05/ 8

S wmax =10%) when the contact angle is 60°. In summary, the similarity of the shapes influences the capillary behavior primarily for shapes with the same inscribed radius and the same shape factor in the same wetting conditions. The difference is controlled by the value of m. Saturation is more sensitive to the value of m than the capillary entry pressure is. Additionally, the effect of m increases as the contact angle increases. 4.3 The further comprehension of m The value of m for different wetting conditions is shown in Fig. 8. First, m decreases as both the shape factor and the contact angle increase. The contact angle influences m much more than the shape factor does. Second, the data are more scattered when the contact angle is large, especially for shapes with small shape factors. However, the curve converges better when the shape factor increases sufficiently. Finally, the data in Fig. 6 and Fig. 7 show that as m decreases, the capillary entry pressure and the wetting phase saturation also decrease, i.e., there is less capillarity in the corners.

Fig. 8 The value of m for different wetting conditions.

The threshold radius configurations in a triangular pore space are shown in Fig. 9 to demonstrate the physical meaning of m. The value of m is 3.1416, 2.8698, 2.2854 and 1.3261. First, the threshold radius increases as the contact angle increases, i.e., the space occupied by the wetting phase in the corners decreases gradually, finally disappearing. More details of the process are shown in Fig. 10. The inscribed circle center and the vertexes of the polygon are used to construct three segments, i.e., the angular bisectors (the blue lines in Fig. 10). The centers of the sectors (red lines in Fig. 10) in the corners move along the segments when the contact angle increases. Second, the sector in the corner intersects the extended lines of the polygon side when there is no wetting phase in the corner (see Fig. 9 (d)). Finally, the sectors in the corners splice into a complete circle, while the parts that lie outside the polygon (the green sections in Fig. 9) do not appear in the reduced shape. Instead of π, m is introduced to calculate the area of the remaining parts of the circle in the reduced shape, as follows: Area _ outside + mrd 2 = πrd 2

(17)

where Area _ outside is the area of the parts that lie outside of the polygon. According to Eq. 8, the corner with the small angle holds the wetting phase easily, which is 9

confirmed in Fig. 9. The interior angles of the triangle are 72°, 70° and 38°, and the wetting phase fills only the corner with the smallest angle when the contact angle is 60°. Another key point is that the corner with smaller angle holds more wetting phase in any wetting condition. In a polygon, the average of the interior angles is as follows: Ang = ( n − 2)180 / n

(18)

such that a higher-order polygon tends to have larger interior angles and thus less capillarity in the corners. That is the reason that there is more capillarity in the corners of triangles. However, there will be more capillarity in the corners of the convex ones when the polygon is concave. This is because the concave angle (>180°) in the polygon will reduce the other angles sharply, e.g., the 4-sided concave quadrangle has more capillarity in the corners than the convex quadrangle or the triangle.

Fig. 9 Threshold radius configurations in a triangular pore space; the blue sections indicate the wetting phase.

Fig. 10 Changes of the threshold radius in a triangular pore space. 10

5. Application The RSG method is a special case of the MS-P method. The conclusions and equations from our RSG model apply only to tangential polygons. In this section, we describe the application of the RSG method to a practical problem and several general cases. 5.1 The new predictor of the threshold radius In the perfect wetting case Eq. (13) is written as follows:

，

1 1 1 1 1 = + = + rd 2rh rarea rin rarea

(19)

The threshold radius is half of the harmonic average of inscribed radius (or two times the hydraulic radius) and area-equivalent radius. This equation is similar to an equation of electric resistance in a parallel circuit. Eq. (19) can explain the conclusion by Lindguist [2006] easily. No matter which one (the inscribed radius or the area-equivalent radius) is selected as the predictor of the threshold radius, Eq. (19) would lose one term and hence be smaller, causing the threshold radius to inevitably be over-predicted by the selected predictor.

Fig. 11 Surface plots of m, the shape factor and the contact angle.

With contact angles considered, Eq. (19) is modified with weighting factors to become Eq. (13). However, the calculation of m in Eq. (13) needs information about the shape geometry, which would be time-consuming. Therefore, we calculate m for many quadrangles in different wetting conditions and fit the relation between m, the shape factor and the contact angle as follows: m = p00 + p10G + p01θ + p20G 2 + p11Gθ + p02θ 2

(20)

where G is the shape factor, θ is the contact angle, and p00=3.161, p10=3.034, p01 = 0.0008066, p20 = -38.9, p11=-0.1084 and p02 =-0.0004336. If m<0, m is set equal to 0. The surface plot of Eq. 20 is shown in Fig. 11. The data show that the value of m decreases as both the shape factor and the contact angle increase. According to Eq. 13 and Eq. 20, we can predict the threshold radius easily. The new parameter thus found is as follows:

11

 cos θ m/π + )  rnew = 1/ ( 2rh rarea  (21)   r ' = 1/ ( cos θ + m / π )  new rin rarea  Because there is not always an inscribed circle tangent to all sides in arbitrary polygons, the inscribed radius can be extended to the radius of the maximum circle inside the polygon. Therefore, we still need to look for the maximum circle when we use rnew ' to predict the threshold radius. For convenience, we usually select rnew to express the new predictor because it only requires the area, the perimeter and the contact angle as the input parameters. 5.2 From the specific to the general 5.2.1 Verification for arbitrary n-sided polygons The new predictor is obtained for tangential polygons; whether it is accurate enough for arbitrary shapes is unknown. Lago et al. [2001] derived the equation of threshold pressure in capillaries with polygonal cross sections. We choose arbitrary n-sided polygons to test the accuracy of rnew and the true threshold radius is calculated referring to Lago’s method. A method to generate arbitrary n-sided polygons is introduced first. The basic idea is to select n points in a unit circle as the vertices of arbitrary polygons. The specific steps are the following. 1) Divide the interval [0,2π] into n equal parts of a unit circle. Select a radian value αi in every interval as follows:

αi =

2π 2π (i − 1) + rand , i=1, 2, 3…n, where rand ∈ (0,1) and is random. n n

2). Let ri = rand , i=1,2,3…n, where rand ∈ (0,1) and is a random number. 3). Let xi = ri cos αi , yi = ri sin αi , i=1, 2, 3…n. Choose the coordinates (xi, yi) as the vertex coordinates of the polygons; then the n-sided polygons are uniquely determined. Arbitrary polygons generated according to this algorithm are in a unit circle. The shapes can be either concave or convex. The shapes are irregular and are not necessarily tangential polygons. After the polygons are constructed, the coordinates of every vertex are known. Then, the geometric parameters of the arbitrary polygons are calculated by coordinate operations. First, the perimeter is calculated according to the distance formula between two points. Second, the n-sided (n > 3) polygons are divided into n-2 triangles and the area is calculated according to Heron's formula. Finally, whether the interior angle is concave is judged by the relation between point and line, and the angle is obtained by calculating the included angle of the vectors. For example, in Fig. 12, the area calculation of quadrangle ABCD is found by calculating the areas of the triangles ABC and ACD. Inserting point C into the straight line equation of BD, we obtain a positive value, indicating that the interior angle BCD is concave.

12

Fig. 12 Schematic showing the calculation of the geometric parameters in a quadrangle.

Let n=4,5,6 and generate 5000 polygons randomly. Then calculate the true threshold radius and threshold pressure for each shape. Fig. 13 shows the distribution of the radii when the contact angles are 0° and 45°. The hydraulic radius tends to underestimate the threshold radius and the area-equivalent area radius overestimates it. The new predictor more closely predicts the threshold radius distribution when compared with the hydraulic radius distribution. Both the hydraulic radius and the new predictor predict the threshold radius more accurately when the contact angle is 0°. However, the hydraulic radius becomes gradually less accurate as the contact angle increases. The relative error when the contact angle is 45° is shown in Fig. 14. The relative error of the new parameter is less than 10% compared with the threshold radius while the error in the hydraulic radius is 30%.

Fig. 13 Different radius distributions when the contact angle is 0° and 45°.

The new predictor is proven to be more accurate from another aspect also. Fig. 15 shows the relation between the threshold radius predictors and the threshold pressure. In theory, there is a one-to-one correspondence between the threshold radius and the threshold pressure. The curve relating the new predictor to the threshold pressure is less scattered than the same curve for the hydraulic radius and threshold radius, indicating that the new predictor is more accurate than the hydraulic radius.

13

Fig. 14 The relative error of the threshold radius predictor when the contact angle is 45°.

Fig. 15 The relation between the threshold radius predictors and threshold pressure when the contact angle is 45°.

In summary, the new predictor works well for arbitrary shapes, regardless if the shape is convex or concave, regular or irregular. Therefore, our method can predict the capillary entry pressure accurately using only the area, the perimeter and the contact angle as input parameters. 5.2.2 Verification using a real rock sample Many methods have been used to calculate the fluid configuration of the cross-section during primary drainage based on the pore image. A semi-analytical model for the computation of the capillary entry pressures and fluid configurations in pore spaces from 2D rock images was proposed by Frette et al. [2010], and the configurations were calculated for Bentheimer sandstone pores. Lindquist [2006] used computational geometry theory to calculate the configurations for a Fontainebleau core sample image. The capillary entry pressure and the threshold radius of the real shape can be calculated by the MS-P method according to the image configuration. Configurations of real shapes from different papers are selected to verify the proposed threshold radius predictor [Frette et al., 2010; Helland et al., 2008; Lindquist, 2006]. The contact angle is 0 and 15 shapes are chosen. Aeff , Leff , rin and other geometric parameters are measured by image processing techniques; some examples are shown in Fig. 16. The comparisons of the threshold radii calculated from the configurations with different threshold radius predictors are shown in Fig 17. The results show that both the proposed predictors are more accurate and the hydraulic radius method is less accurate. 14

Helland et al. [2008] used the dimensionless hydraulic radius to judge the concave cross-sections. The dimensionless hydraulic radius is defined as H = rh / rin . If H<0.5, then the shape is non-convex; otherwise, the shape is convex. The four data points inside the ellipse in the Fig. 17 correspond to the four shapes in Fig. 16. These four shapes are all concave, which may be the reason why the hydraulic radius fails to predict the threshold radius accurately.

Fig. 16 Fluid configurations of real shapes during primary drainage; the blue sections indicate the wetting phase (Frette et al., 2010; Helland et al., 2008).

Fig. 17 Different threshold radius predictors in real rock shapes.

6. Conclusions We extend the RSG method to tangential polygons. By introducing an intermediate variable m, the threshold radius equation not only considers the contact angles and the shape information but also has a concise form. The parameter m indicates the completeness of the circle in the reduced shape. It is determined by the shape geometry and the contact angle. The effect of m indicates the effect of the similarity of the shapes when the contact angle and the shape factor are fixed. This effect increases with the increase of the contact angle and the decrease of the shape factor. The value of m influences the phase saturation more than it influences the capillary entry pressure. The maximum relative 15

difference of the threshold wetting phase saturation caused by m is approximately 0.05, while that of the normalized threshold capillary entry pressure is approximately 0.1. A more accurate parameter considering the contact angle and shape information is proposed as a threshold radius predictor. Even though it is derived from convex polygons with inscribed circles, it is proven to apply to irregular polygonal shapes and real rock shapes also. The new predictor is better than the hydraulic radius in predicting the capillary behavior in concave shapes. Moreover, the parameter has some physical meaning. It is half of the harmonic average of the inscribed radius (or two times the hydraulic radius) and the area-equivalent radius with the contact angle as the weighting factor. This study contributes to the understanding of capillary behavior and is an addition to the theory of capillarity. Further work on the imbibition process, the dynamic capillary pressure and pore network modeling can be performed based on these results. Acknowledgments The authors greatly appreciate the financial support of the National Natural Science Foundation of China (Grant no. 10972237), the Important National Science and Technology Specific Projects of China (Grant no. 2011ZX05011), the Natural Science Foundation for Distinguished Young Scholars of Shandong Province, China (Grant no. JQ201115), the Program for New Century Excellent Talents in University (Grant no. NCET-11-0734), the Fundamental Research Funds for the Central Universities (Grant no. 13CX05007A), and the Program for Changjiang Scholars and Innovative Research Team in University (Grant no. IRT1294).

References

‐

‐

1. Al Omari, A., and Masad, E.: Three dimensional simulation of fluid flow in X ray CT images of porous media. Int. J. Numer. Anal. Met. 28(13), 1327-1360 (2004). 2. Blunt, M.J.: Pore level modeling of the effects of wettability. SPE J. 2(4), 494-510 (1997). 3. Blunt, M.J.: Flow in porous media—pore-network models and multiphase flow. Curr. Opinion Colloid Interface Sci. 6(3), 197-207 (2001). 4. Das, D. B., Mirzaei, M., and Widdows, N.: Non-uniqueness in capillary pressure–saturation–relative permeability relationships for two-phase flow in porous media: interplay between intensity and distribution of random micro-heterogeneities. Chem. Eng. Sci. 61(20), 6786-6803 (2006). 5. Das, D. B., and Mirzaei, M.: Dynamic effects in capillary pressure relationships for two-phase flow in porous media: experiments and numerical analyses. AIChE J., 58(12), 3891–3903 (2012). 6. Frette, O. I., and Helland, J. O.: A semi-analytical model for computation of capillary entry pressures and fluid configurations in uniformly-wet pore spaces from 2D rock images. Adv. Water Resour. 33(8), 846-866 (2010). 7. Helland, J.O., Ryazanov, A.V., and van Dijke, M.I.J.: Characterization of Pore Shapes for Pore Network Models. In : Proceedings of the 11th European Conference on the Mathematics of Oil Recovery (ECMOR XI), Bergen, Norway, September 2008. 8. Hanspal, N. S., Allison, B. A., Deka, L., and Das, D. B.: Artificial neural network (ANN) modeling of dynamic effects on two-phase flow in homogenous porous media. J. Hydroinform., 15(2), 540-554 (2013). 16

9. Iliev, O., Efendiev, Y., Li, J., Calo, V., and Brown, D.: Multiscale lattice Boltzmann method for flow simulations in highly heterogenous porous media. In SPE Reservoir Characterization and Simulation Conference and Exhibition. Society of Petroleum Engineers, September, 2013. 10. Jia, P., Dong, M., and Dai, L.: Threshold pressure in arbitrary triangular tubes using RSG concept for all wetting conditions. Colloids Surf. A: Physicochem. Eng. Aspects 302(1), 88-95 (2007). 11. Joekar-Niasar, V., Doster, F., Armstrong, R. T., Wildenschild, D., and Celia, M. A.: Trapping and hysteresis in two-phase flow in porous media: A pore-network study. Water Resour. Res. 49(7), 4244-4256 (2013). 12. Kutay, M.E., Aydilek, A.H. and Masad, E., 2006. Laboratory validation of lattice Boltzmann method for modeling pore-scale flow in granular materials. Comput. Geotech. 33(8), 381-395. 13. Lago, M., and Araujo, M.: Threshold pressure in capillaries with polygonal cross-section. J. Colloid Interface Sci. 243(1), 219-226 (2001). 14. Lindquist, W.B.: The geometry of primary drainage. J. Colloid Interface Sci. 296(2), 655-668 (2006). 15. Mayer, R.P., and Stowe, R.A.: Mercury porosimetry—breakthrough pressure for penetration between packed spheres. J. Colloid Interface Sci. 20(8), 893-911 (1965). 16. Mason, G., and Morrow, N.R.: Capillary behavior of a perfectly wetting liquid in irregular triangular tubes. J. Colloid Interface Sci. 141(1), 262-274 (1991). 17. Ma, S., Mason, G., and Morrow, N.R.: Effect of contact angle on drainage and imbibition in regular polygonal tubes. Colloids Surf. A: Physicochem. Eng. Aspects 117(3), 273-291 (1996). 18. Man, H.N., and Jing, X.D.: Pore network modelling of electrical resistivity and capillary pressure characteristics. Transp. Porous Media 41(3), 263-285 (2000). 19. Mirzaei, M., and Das, D. B.: Dynamic effects in capillary pressure–saturations relationships for two-phase flow in 3D porous media: Implications of micro-heterogeneities. Chem. Eng. Sci., 62(7), 1927-1947 (2007). 20. Nassehi, V., B. Das, D., M. T. A. Shigidi, I., and J. Wakeman, R.: Numerical analyses of bubble point tests used for membrane characterisation: model development and experimental validation. Asia-Pac. J. Chem. Eng., 6(6), 850–862 (2011). 21. Nejad Ebrahimi, A., Jamshidi, S., Iglauer, S., and Boozarjomehry, R. B.: Genetic algorithm-based pore network extraction from micro-computed tomography images. Chem. Eng. Sci. 92, 157-166 (2013). 22. Oren, P.E., Bakke, S., and Arntzen, O.J.: Extending predictive capabilities to network models. SPE J. 3(4), 324-336 (1998). 23. Princen, H.M.: Capillary phenomena in assemblies of parallel cylinders: I. Capillary rise between two cylinders. J. Colloid Interface Sci. 30(1), 69-75 (1969a). 24. Princen, H.M.: Capillary phenomena in assemblies of parallel cylinders: II. Capillary rise in systems with more than two cylinders. J. Colloid Interface Sci. 30(3), 359-371 (1969b). 25. Princen, H.M: Capillary phenomena in assemblies of parallel cylinders: III. Liquid columns between horizontal parallel cylinders. J. Colloid Interface Sci. 34(2), 171-184 (1970). 26. Patzek, T.W.: Verification of a complete pore network simulator of drainage and imbibition. SPE J. 6(2), 144-156 (2001). 27. Patzek, T.W., Silin, D.B.: Shape factor and hydraulic conductance in noncircular capillaries: I. one-phase creeping flow. J. Colloid Interface Sci. 236(2), 295-304 (2001). 17

‐

28. Pan, C., Hilpert, M., and Miller, C.T.: Lattice Boltzmann simulation of two-phase flow in porous media. Water Resour. Res. 40 (1), W01501 (2004). 29. Panfilov, M., and Panfilova, I.: Phenomenological meniscus model for two-phase flows in porous media. Transp. Porous Media, 58(1-2), 87-119 (2005). 30. Piri, M., and Blunt, M.J.: Three-dimensional mixed-wet random pore-scale network modeling of two-and three-phase flow in porous media. I. Model description. Phys. Rev. E 71(2), 026301 (2005). 31. Wolf-Gladrow, D.A., 2000. Lattice-gas cellular automata and lattice Boltzmann models: an introduction (No. 1725). Springer, 2000. 32. Zaretskiy Y., Geiger, S., Sorbie, K. and Förster, M.: Efficient flow and transport simulations in reconstructed 3D pore geometries. Adv. Water Resour. 33(12), 1508-1516 (2010). 33. Zhou, Y., Helland, J. O., and Hatzignatiou, D. G.: Computation of three-phase capillary entry pressures and arc menisci configurations in pore geometries from 2D rock images: A combinatorial approach. Adv. Water Resour. 69, 49-64 (2014).

1. We extend the reduced similar geometry (RSG) method to tangential polygons 2. The effect of the shape similarity on the capillary behavior is investigated. 3. A new parameter is proposed to predict the capillary entry pressure.

18

Graphical Abstract (for review)