Subscriber access provided by TUFTS UNIV
Article
The investigation of dynamic contact angle using a direct numerical simulation method Guangpu Zhu, Jun Yao, Lei Zhang, Hai Sun, and Aifen Li Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.6b02543 • Publication Date (Web): 18 Oct 2016 Downloaded from http://pubs.acs.org on October 23, 2016
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
Langmuir is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Langmuir
The investigation of dynamic contact angle using a direct numerical simulation method Guangpu Zhu, Jun Yao*, Lei Zhang*, Hai Sun, Aifen Li School of Petroleum Engineering, China University of Petroleum, Qingdao 266580, China
Abstract A large amount of residual oil, which exists as isolated oil slugs, remains trapped in reservoirs after water flooding. Numerous numerical studies are performed to investigate the fundamental flow mechanism of oil slugs to improve flooding efficiency. Dynamic contact angle models are usually introduced to simulate an accurate contact angle and meniscus displacement of oil slug under a high capillary number. Nevertheless, in the oil slug flow simulation process, it is unnecessary to introduce the dynamic contact angle model due to a negligible change of meniscus displacement after using the dynamic contact angle model when the capillary number is small. Therefore, a critical capillary number should be introduced to judge whether the dynamic contact model should be incorporated into simulations. In this study, a direct numerical simulation method is employed to simulate the oil slug flow in a capillary tube at the pore scale. The position of the interface between water and oil slug is determined by using the phase field method. The capacity and accuracy of the model is validated by a classical benchmark: a dynamic capillary filling process. Then different dynamic contact angle models and the factors which affect the dynamic contact angle are analyzed. The meniscus displacements of oil slug with dynamic contact angle and static contact angle are respectively obtained during simulations and the relative error between them are calculated automatically. Relative error limit has been defined to 5%, beyond which the dynamic contact angle model needs to be incorporated into the simulation to approach the realistic displacement. Thus, the desired critical capillary number can be determined. A three dimensional universal chart of critical capillary number, which functions as static contact angle and viscosity ratio, is given to provide a guideline for oil slug simulation. Also, a fitting formula is presented for ease of use. Keywords: oil slug, dynamic contact angle, direct numerical simulation, phase field method, critical capillary number.
1 Introduction Slug flow in a capillary tube is one of the simplest and most important two-phase flow patterns that attract considerable research interest due to numerous practical applications in the petroleum industry.1 A large amount of residual oil remains trapped in reservoirs after water flooding,2 and the residual oil is thought to exist as isolated oil slugs held in place by the combination of capillary forces and the resistance of pore walls of the porous reservoirs. To improve the efficiency of water flooding, it’s crucial to comprehend the fundamental mechanism of mobilizing oil slugs in pore structure.3, 4 Numerous experimental studies of oil slug in a capillary tube under high capillary numbers, that are several orders higher than the range of capillary numbers encountered in flow through porous media in oil reservoirs, have been reported in literatures and great achievements have been made in this field.3 Furthermore, experiments were carried out on the mobilization and very slow flow of oil slug in a capillary tube, as well as the relationship between pressure drop and flow rate was investigated. However, the experimental
ACS Paragon Plus Environment
Langmuir
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
findings are limited by the available experimental conditions and facilities.5 With this consideration, the numerical approach is used as a powerful tool to investigate fundamental mechanisms of oil slug at the pore scale and demonstrates significant advantages over experimental studies. Several numerical studies based on computational fluid dynamics have been conducted to study the mobilization of oil slug and show good agreement with the experiment. Nevertheless, the experimentally observed contact angle generally differs from its static contact angle (SCA). Also, the meniscus displacement of oil slug obtained in the numerical simulation (using a static contact angle) differs from that of actual experiment.6-11 It is necessary to introduce the dynamic contact angle model to explain this phenomenon. The traditionally used static contact angle is a unique property of the material system, and the relationship between static contact angle and surface tension of solid/vapor and solid/liquid interface can be depicted by Young’s equation for an ideal solid surface which is flat, rigid, perfectly smooth, and chemically homogeneous.9, 12, 13 The contact angle at a moving contact line observed experimentally is called the dynamic contact angle,14-16 which is not only a material property and usually differs from its equilibrium value or static contact angle. The main parameter used to quantify the dynamic wetting is the relative velocity at which the liquid moves along the wall, called moving contact line velocity. Two main theories have been proposed based on numerous studies on dynamic contact angle. One of these theories, known as the hydrodynamic theory, emphasizes energy dissipation due to viscous flow within the wedge of liquid near the moving contact line.17-19 Three relevant length scales are introduced in this theory, namely microscopic scale, mesoscopic scale and macroscopic scale. Short-range intermolecular forces govern the microscopic angle thereby retain its static value. In a contrasting viewpoint provided by the molecular-kinetic theory,7, 20-22 the motion of the contact line is determined by the statistical dynamics of the molecules within a three-phase zone where the solid, liquid and gas phases meet. Therefore, the microscopic contact angle is velocity-dependent and identical with the experimentally observed angle. Only two length scales are contained in this theory, including the molecular scale where the energy dissipation occurs and the macroscopic scale where its effects are seen. Numerous dynamic contact angle models are proposed based on the two theories and experiments, such as the Cox model,23 the Blake model21 and the Jiang model.21, 24, 25 These models are predominately used for capillary phenomena, such as oil slug flow in a capillary tube, meniscus rise in a vertical or horizontal tube and droplet spreading, as well as coating.24, 25 Most of these models need adjustable parameters to match the experimental dynamic contact angle and some of them cannot be applied for the entire range of contact angle and other parameters. Therefore, to approach the realistic value of contact angle and obtain exact meniscus displacement, it’s necessary to incorporate the dynamic contact angle model, which considering all the effects that influence the capillary phenomenon, into the computational fluid dynamic (CFD) simulation of oil slug in a capillary tube. The introduction of dynamic contact angle model makes it possible to obtain accurate contact angle and meniscus displacement. For some capillary geometry25 (e.g. oil slug displacement in a capillary tube and meniscus rising in a horizontal tube), the meniscus displacement is a primary concern and the accuracy of contact angle can be overlooked. In some cases, (e.g. large contact angle and large capillary number), the difference of meniscus displacement is negligible between dynamic contact angle model and static contact angle. Hence, it is unnecessary to incorporate the dynamic contact angle into computational fluid dynamics simulation. Also, introduction of
ACS Paragon Plus Environment
Page 2 of 18
Page 3 of 18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Langmuir
dynamic contact angle increases the nonlinearity and computation cost of numerical model, which may lead to non-convergence. However, with the increasing in capillary number, the difference between simulated meniscus displacement by using dynamic contact angle and that by using static contact angle cannot be neglected, and the meniscus simulated by using a static contact angle may be incorrect. Therefore, a critical value or a limit, as a function of capillary number or some other factors, should be defined for the pore scale simulation of oil slug in a capillary tube. The critical value should be able to provide a reference about whether the dynamic contact model should be incorporated in the computational fluid dynamics simulation. To the authors’ knowledge, no research has been reported to determine the critical capillary number. The critical value is supposed to be defined on the bases of large amount of pore scale simulation26, 27 of oil slug in a capillary tube. In this study, direct simulation method is used in the pore scale simulation. Pore scale flow simulation is conducted to explain the phenomena at the macro scale or identify macroscopic properties of porous media.28-31 Several approaches have been developed for pore scale simulation of multiphase flows in the porous media. In this study a direct numerical simulation method (DNS) are used to explore the critical value. The main advantage of the DNS method is the natural and consistent way that the flow models can be extended to model processes other than just the flow.32, 33 In conjunction with high-performance computing, well established numerical methods used in computational fluid dynamics, such as finite volume and finite element methods, have also become practical for direct numerical simulation of flow in the complex geometry of heterogeneous pore space.30 However, the physics of the interface movement are not included in Navier-Stokes (N-S) equations. Consequently, for multiphase flow, N-S equations should be coupled with an interface tracking method. The Level set method, Volume of fluid method (VOF) and phase field method are the three commonly used methods for interface tracking. The advantage of the Level set model is that one can perform numerical computations involving curves and surfaces on a fixed Cartesian grid without having to parameterize these objects.30, 34 However, Level set method often significantly violates mass conservation. The VOF method is a Eulerian method that uses a volume fraction indicator to determine the locations of the interfaces of different phases in all cells of a computational domain.30, 35 VOF method makes it very easy to follow shapes that change topology, but it is less elegant than the Level set method in terms of computational accuracy. Also, the phase field method treats the interface as a physically diffuse thin layer. The interface differs from other fixed-grid methods, where the interface is shaped conceptually but regularized numerically by spreading the interface force over a volume. Phase-field parameters are introduced to distinguish the phases and interface, and are governed by the Cahn-Hilliard nonlinear diffusion equation.36-38 Moreover, the phase field method can guarantee mass conservation in the simulation. Considering numerical instability arising at the interface region when the interfacial tension becomes a dominant factor if using the VOF method or Level set method for pore scale simulation, phase field method is applied to track the two phase interface during our simulation. In this study, a direct numerical simulation method is employed to simulate the oil slug flow in a capillary tube at the pore scale. The position of the interface between water and oil slug is determined by using the phase field method. The capacity and accuracy of the model is validated by a classical benchmark: a dynamic capillary filling process. Different dynamic contact angle models and the factors which affect the dynamic contact angle are analyzed. The meniscus displacements of oil slug with dynamic contact angle and static contact angle are respectively
ACS Paragon Plus Environment
Langmuir
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 18
obtained during simulations and the relative error between them are calculated automatically. Relatively error limit has been defined to 5%, beyond which the dynamic contact angle needs to be incorporated into the simulation to approach the experimentally observed displacement. Thus, the desired critical capillary number can be determined. A 3D universal chart of critical capillary number, which functions as static contact angle and viscosity ratio is given to provide some guidance for oil slug simulation. Also, a fitting formula is presented for ease of use.
2 Numerical simulation 2.1 Governing equations The continuity equation satisfies the conservation of mass for fluid flow. For an incompressible and Newtonian fluid, the continuity equation is given by27 ∇ ⋅u = 0 (1) The N-S equation represents the conservation of momentum for the fluid flow. For an incompressible fluid, it can be written as34 T ∂u ρ + ( u ⋅∇ ) u = ∇ ⋅ − pI + µ ∇u + ( ∇u ) + Fst (2) ∂t Here u is the flow velocity of the fluid, m/s; ρ denotes the density of the fluid, kg/m3; I is the unit vector; t is the time, s; µ is the viscosity of fluid, mPa·s. Fst is the term representing interfacial tension farce, N/m. The gravity is neglected since the entire fluid volume is assumed to be at the same gravitational potential.
(
)
2.2 Phase field method During the flow process of oil slug, the interface evolves according to the flow. To track the interface of two immiscible fluids accurately, the phase field method is introduced. The interface differs from the Level set method where the interface is shaped conceptually but treated as a physically diffuse thin layer.37, 38 The dynamics of the thin layer are affected by convection and diffusion. In some high resolution problems, the layer must be adequately thin to track the real physical interface, which can be achieved via local grid refinement and adaptive meshing technique. The properties in each phase are constant but vary continuously across the mixing layer. Compared to other interface tracking method, the phase field method is more suitable for numerical simulation. The phase field variable is introduced to distinguish the phases and interface that obey a Cahn–Hilliard type of convection-diffusion equation,39 ∂φ γλ ∂t + u ⋅∇φ = ∇ ⋅ ε 2 ∇ψ ψ = −∇ ⋅ ε 2 ∇φ + (φ 2 − 1) φ
(3)
Here, φ denotes a nondimensional phase-field variable, that is equal to 1 in one phase and -1 in the other and varies continuously in the two phase interface region; ε is the interface thickness, m; γ is the mobility, which denotes the mobility velocity of the interface under a unit driven force, m3·s/kg; λ is the magnitude of the mixing energy, N. The surface tension acting on the fluid-fluid interface can be written as40 (4) Fst = G∇φ where, G is the chemical potential defined as
ACS Paragon Plus Environment
Page 5 of 18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Langmuir
φ (φ 2 − 1) λ 2 = ψ G = λ −∇ φ + ε2 ε2
(5)
Note that the fluid density and viscosity in the above equation are defined as functions of the phase field variable, as follows,37, 40 ρ=
1+ φ 1−φ ρn + ρw 2 2
(6)
µ=
1+ φ 1− φ µn + µw 2 2
(7)
where the subscripts (w) and (n) denote the wetting phase and non-wetting phase, respectively. 2.3 Dynamic contact angle model The theoretical and experimental studies on dynamic contact wetting have produced a number of dynamic contact angle models, including hydrodynamic models,41 molecular kinetic models,21 experimental models,24 as well as hybrid models.9 All of these models are expressed as function of contact line speed.24 Various dynamic contact models are presented in this section. Based on theoretical analysis, Cox18, 23 gives an essential relationship between the contact line velocity and the dynamic contact angle, −1 1 Q1 Q2 −1 Ca = ln ( ε ) − + g (θ D ) − g (θ S ) + O ln (ω −1 ) f (θ D ) f (θ S )
3
(8)
where Ca=µv/σ is the capillary number, v is the velocity of contact line, m/s, σ is the surface tension, N/m; θD is the dynamic contact angle, θS is the static contact angle; Q1, Q2 are constants associated with the outer flow and with the wall slip; ω is a small dimensionless parameter associated with the micro-region of the contact line, f (·) and g (·) are functions obtained in the reference.42 Most of the existing empirical dynamic contact angle models are simplification of the equation (8). The most widespread dynamic contact angle is given by Hoffman-Voinov-Tanner law for small capillary number and it is called the Cox model in this study,43 θ D3 − θ S3 = 144Ca
(9)
Jiang et al18, 25 give an empirical correlation of dynamic contact angle, which depends on the contact line velocity, and this model fits the data equally well, cos θ D − cos θ S = − tanh ( 4.96Ca 0.702 ) (10) 1 + cos θ S Kalliadasis24 proposed a dynamic contact angle model as tan θ D = 7.48Ca1/3 − 3.28α 0.04 Ca 0.293
(11)
where α=10-8. Bracke et al22 distil an operational equation based on the dynamic contact angle measurement when a continuous solid stripe is drawn into a large pool of liquid. The model shows its universality on the wetting process of a vertical solid plate and the spreading of a liquid drop on a solid. The model is shown as follows, cos θ D = cos θ S − 2 (1 + cos θ S ) Ca 0.5
Kistler18 proposed a dynamic contact model as,
ACS Paragon Plus Environment
(12)
Langmuir
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
0.706 Ca + fHI cos θ = 1 − 2 tanh 5.16 D 1 + 1.31( Ca + fHI )0.99 0.706 θS 1 fHI0.706 = (1 + 1.31fHI0.99 ) tanh −1 1 − cos 5.16 2
Page 6 of 18
(13)
where fHI is the inverse function of the Hoffman’s function fHoff (·) defined as 0.706 x f Hoff ( x ) = arccos 1 − 2 tanh 5.16 0.99 1 + 1.31x
(14)
The models mentioned above belong to typical hydrodynamic and empirical models, which attribute the changes in the experimentally observed dynamic contact angle to the viscous bending of the two-phase interface within a mesoscopic region. Whereas, according to the contrasting viewpoint provided by the molecular kinetic theory, the motion of the contact line is determined by the molecules statistical dynamics within a three-phase zone. The model is based on the idea that the driving force for the contact line to move in a given direction is equal to the out-of-balance surface tension force when equilibrium is disturbed.9 The resulting equation for the contact line velocity is,9, 20 σ ( cos θ S − cos θ D ) d 2 v = 2κ 0 λ sinh 2kBT
(15)
If the argument of the sinh function is small, equation (16) can be reduced to its linear form cos θ D = cos θ S −
k BTv σκ 0 d 3
(16)
where κ0 is the characteristic frequency, kB is the Boltzmann constant, T is the absolute temperature, d is the distance between adsorption sites on the solid surface.
Despite fundamentally different physics between hydrodynamic and molecular kinetic models, both models have been shown to be reasonably effective in describing the experimentally observed behavior of the dynamic contact angle. Recently, a point of view has been increasingly focused that both wetting-line friction and viscous dissipation play a part in determining the dynamic contact angle. Hence, a hybrid model is proposed by Gennes et al,44 as is shown in equation (17), cos θ D = cos θ S −
v kBT 6u L + ln σ κ 0 d 3 θ D Lm
(17)
where L and Lm are appropriately chosen macroscopic and microscopic length scales, respectively. In this study, only hydrodynamic models and empirical models are applied in our simulations due to the difficulty in determining the parameters of the molecular kinetic models. 2.4 Numerical solution The above equations are solved using the CFD module of COMSOL Multiphysics (commercial software, which can directly implement partial differential equations and initial-boundary conditions). The mobility tuning parameter χ and interface thickness are set to 1 and half the maximum element size, respectively. An initial time step of 1×10-7 s is selected and the time step is allowed to automatically adjust based on the backward differentiation formula (BDF). A convergence criterion of 0.001 is specified to control the iterative solution process.
ACS Paragon Plus Environment
Page 7 of 18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Langmuir
PARDISO is used as a system solver in this work.
3 Results and Discussion 3.1 Validation of numerical model This section uses a simple benchmark to validate our numerical model. The capillary filling is used to assess whether the proposed multiphase model can simulate moving contact line problems, which originated from the pioneering works of Washburn and Lucas.39, 45 The wetting fluid filling into a horizontal capillary tube is controlled by capillary force and viscous force. As already remarked in many works, if the gravity and inertial effects are neglected, the relation for the predicting the position of the meniscus can be expressed as, 6 r
σ cos(θ ) = µw x + µn ( L − x )
dx dt
(18)
Where σ is surface tension, N/m; θ is the contact angle; r is the diameter of capillary tube, m; µ w and µn are the dynamic viscosities of the wetting and non-wetting phase, respectively, Pa·s; L is the length of capillary tube, m. The equation (18) can be transformed into the integration form,
µ w − µn 2
(x
2
− x02 ) + µn L ( x − x0 ) =
rσ cos(θ ) t 6
(19)
A domain, as shown in Figure 1, with the dimension of 700 µm×10 µm, can be divide into two parts. In the middle portion, the boundaries of the capillary tube are non-slip and wetting, with a given contact angle θ=π/3. The length of the capillary tube taken as L=300 µm. The left and right portions of the domain have top and bottom periodic boundary conditions so as to support a perfect flat two-phase interface, mimicking a “infinite reservoir”. The length of left and right parts of the domain are Lleft=Lright=200 µm. The initial length of wetting phase in the right part of the domain Lw=130 µm. The periodic boundary conditions in x direction are also imposed at the both ends of the domain to ensure total conservation of mass inside the system.
Figure 1 Schematic of capillary filling
The densities of fluids used in the simulations are ρw=ρn=1000 kg/m3, the surface tension σ is 0.04 N/m. The viscosity ratio M is 0.1 (µn=1 mPa·s, µw=10 mPa·s). In Figure 2, we show the behavior of meniscus position x(t), as a function of time for a given contact angle, a given surface tension and a given viscosity ratio. The numerical solutions agree quite well with the analytical solutions, which indicates the accuracy and capacity of our model in dealing with the contact line problem.
ACS Paragon Plus Environment
Langmuir
60
The position of meniscus (µm)
50
analytical solution numerical solution
40
30
20
10
0 0.0
-4
2.0x10
4.0x10
-4
6.0x10
-4
8.0x10
-4
1.0x10
-3
t (s)
Figure 2 The numerical and analytical distance of capillary filling
3.2 The analysis of dynamic contact angle models Models based on the molecular kinetic theory, as outlined above, manifest the feasibility in describing the experimentally observed behavior of dynamic contact angle in few cases that are quite typical. However, parameters used in these models are not fixed and are impacted by experiment conditions. Hence, it’s difficult to determine these parameters without experimentally observed data. In this study, only hydrodynamic theory based models and empirical models are taken into consideration, namely the Cox model, the Jiang model, the Kalliadasis model and the Bracke model, as well as the Kistler model. It can be observed that the five dynamic contact models are predominately affected by the capillary number and static contact angle. Before applying them to the oil slug flow simulation in a capillary tube, we will explore the potential effects of capillary number and static contact angle towards the dynamic contact angle for better understanding. 3.2.1 Static contact angle To investigate the effect of static contact angle, the dynamic contact angles of five models are plotted as a function of static contact angle at a constant capillary number, as shown in Figure 3. 90 80
Jiang Cox Bracke Kalliadasis Kistler Static
70
dynamic contact angle （ °）
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 8 of 18
60 50 40 30 20 10 0 10
20
30
40
50
60
70
static contact angle (°)
(a)
ACS Paragon Plus Environment
80
Page 9 of 18
100
dynamic contact angle (°)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Langmuir
80
60
Jiang Cox Bracke Kalliadasis kistler Static
40
20
10
20
30
40
50
60
70
80
static contact angle (°)
(b) Figure 3 Dynamic contact angles with various static contact angles: (a) Ca=5×10-4; (b) Ca=1.25×10-2
Figure 3 demonstrates that the dynamic contact angle calculated by the Kalliadasis model is a constant at a given capillary number, which indicates that this model is independent of the specified static contact angle. The discrepancies between the Kalliadasis model and other models are remarkable, except in some special situations (e.g. low static angle and low capillary number). The Kalliadasis model has nothing to do with the property of material surface, which may be unrealistic in most real situations. Thus, it will not be applied in our later simulations. It’s observed that the difference between dynamic contact angle and static contact angle increases as capillary number increases, and decreases as the static contact angle decreases. This phenomenon is consistent with the experimentally observed results,7, 46 which provide guidance for determining the critical capillary number in the later section. 3.2.2 The effect of dynamic contact angle Based on the above analysis, the difference between the dynamic contact angle and the static contact angle is obvious at a high capillary number. Also, the static contact angle, which reflects the properties of material surface, contributes to the difference. Therefore, to investigate the effect of the dynamic contact angle on the mobilization of oil slug in a capillary tube at the pore scale, the four dynamic contact angle models are employed in the simulations in a high capillary number with various static contact angles.
Figure 4 Schematic of oil slug in a capillary tube
As shown in Figure 4, the length of the capillary tube L=400 µm, the radius of the tube is r=20 µm, an oil slug is located at the position of 50 µm from the inlet, the initial length of the oil slug Lo=50 µm. The pressure boundary is set at the inlet. The pressure of inlet is 6000 Pa, and a pressure of 0 Pa is set at the outlet. Both the viscosities of continuous phase µw and discontinuous phase µo are 20 mPa·s. The densities of continuous phase ρw and discontinuous phase ρo are 1000 kg/m3 and 870 kg/m3. The interface tension between oil and water is 0.04 N/m. Therefore, the capillary number (Ca=µwUin/σ, Uin refers to the average inlet velocity) of the whole flow system is
ACS Paragon Plus Environment
Langmuir
about 1.25×10-2, which is a considerable high capillary number encountered in flow through porous media in oil reservoirs to ensure the big difference between dynamic contact angle and static contact angle. Two different initial static contact angles of 20° and 60° are set at the wall, respectively. The meniscus displacement is recorded during the simulation and the results are shown in Figure 5.
meniscus displacement (µm)
250
200
SCA Cox Bracke Jiang Kistler
150
100
50
0 0.0000
0.0002
0.0004
0.0006
0.0008
0.0010
0.006
0.008
0.010
t (s)
(a) 250
200
meniscus displacement (µm)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 18
150
SCA Cox Bracke Jiang Kistler
100
50
0 0.000
0.002
0.004
t (s)
(b)
Figure 5 The meniscus displacement of oil slug: (a) Ca=1.25×10-2, θs=20°; (b) Ca=1.25×10-2, θs=60°; the results of Kistler and Bracke model are behind other curves
Figure 5 demonstrates the meniscus displacement of oil slug as a function of time when different static contact angles are applied. It can be observed that the meniscus displacement without considering dynamic contact angle is obviously larger than that with considering dynamic angle when the static contact angle is 20°. The oil slug moves approximately 252.5 µm in 0.01 s when static contact angle is set at the wall. The simulated meniscus displacement almost has no difference with considering various dynamic contact models, and the oil slug moves approximately 217 µm in 0.01 s. The relative error between meniscus displacements with and without considering dynamic contact angle model is 17.5%. When the static contact angle is 60°, the oil slug moves approximately 248.8 µm and 245 µm with and without considering dynamic
ACS Paragon Plus Environment
Page 11 of 18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Langmuir
contact models, respectively. The relative error between them is 1.94%, which can be overlooked. We can conclude that the effect of dynamic contact angle towards meniscus displacement is related to the original static contact angle, and the differences between various dynamic contact models are negligible in terms of meniscus displacement. Therefore, it is necessary to incorporate the dynamic contact angle model into the simulation when original static contact angle is 20°, which makes the simulation result approach the realistic meniscus displacement. Whereas, the result without considering the dynamic contact angle model may already be close to the actual value when the static contact angle is 60°. So now comes the question, when does the dynamic contact models need to be incorporated into the simulation of oil slug mobilization. 3.3 The exploration of critical capillary number Based on the above discussion, we know that the dynamic contact angle and meniscus displacement are predominately controlled by the capillary number and the original static contact angle. The difference between meniscus displacements with and without considering dynamic contact angle model widens with the increasing of the capillary number at a constant static contact angle. The dynamic contact angle can be neglected in the situation when the difference is small and the meniscus is the predominate concern. Nevertheless, once the difference reaches to a level, the dynamic contact model must be incorporated into the simulation of oil slug flow, and the corresponding capillary number is the desired critical capillary number. It is necessary to explore the critical capillary number for a constant static angle, which provides extremely important guidance for the simulation of oil slug mobilization. In this section, the capillary number for various static contact angles at a constant viscosity ratio will be explored. The capillary tube at the pore scale used in the simulations is shown in Figure 4, and an oil slug locates at the same position. An initial static contact angle of 40° is chosen as an example to show how to determine a critical capillary number. The viscosities of non-wetting phase (injecting fluid) and wetting phase (oil slug) are 105.64 mPa·s and 20 mPa·s, respectively. Thus, the viscosity ratio M is 5.282. A wide range of pressures are set at the inlets to acquire different capillary numbers. Under a constant pressure difference, the simulated meniscus displacements dDCA and dSCA can be obtained with and without considering dynamic contact angle, respectively. Then the relative error ε can be expressed as, ε=
d SCA − d DCA × 100% d SCA
(20)
Nine groups of simulations are conducted at various pressure differences, and corresponding relative errors are calculated by equation (20), as shown in Figure 6. In this study, as per customary practice, we defined the relative error of 5% as the limit, beyond which the dynamic contact angle needs to be considered. Thus, the desired critical capillary number can be determined.
ACS Paragon Plus Environment
Langmuir
5.4 5.3
relative error(%)
5.2 5.1 5.0
(-1.49, 5%)
4.9 4.8 4.7 4.6 4.5 -1.60
-1.55
-1.50
-1.45
-1.40
logCa
Figure 6 Relative errors with various capillary numbers: the main function of blue arrow is show that the relative error of 5% (y-coordinate) corresponds to the critical capillary number -1.49 (x-coordinate)
Figure 6 shows the relative errors versus capillary number. The critical capillary number in logarithmic form, which corresponds to the relative error of 5%, can be obtained by using the interpolation method and the interpolation result is -1.49. For a simulation, if the capillary number is larger than -1.49, then a dynamic contact angle model needs to be set at the wall to get an exact meniscus displacement. Otherwise, the result obtained by using a static contact angle is enough due to the relative error is small. To obtain the critical capillary numbers at various initial static contact angles, a large number of simulations are conducted and the corresponding relative error is calculated. The results are shown in Figure 7. -1.0
-1.5
M=5.282
-2.0
logCa
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 12 of 18
-2.5
-3.0
-3.5 20
25
30
35
40
45
50
static contact angle(°)
Figure 7 The critical capillary numbers with different static contact angles Figure 7 gives the critical capillary numbers for different static contact angles. It can be observed that the critical capillary number increases with the increasing in static contact angle. Figure 7 provides a tool to judge whether a dynamic contact angle model needs to be incorporated into the simulation according to the capillary number and the static contact angle. For a particular simulation, a point can be located in Figure 7 with a capillary number and a static contact angle. If the point located in the shadow area, dynamic contact angle is not necessary to set at the wall of a capillary tube, and the result based on the initial static contact angle is accurate enough. Otherwise, it is essential to take the
ACS Paragon Plus Environment
Page 13 of 18
dynamic contact angle into the simulation to approach the realistic result. However, Figure 7 is incomplete as a judgment criterion based on the capillary number and the static contact angle for the following reason. A capillary number only reflects the moving contact line velocity, the viscosity of continuous phase and surface tension, and an initial static contact angle reflects the property of material surface, which are not enough to contain all the necessary information of a two phase flow system in a capillary tube. Viscosity ratio is a significant factor for a two phase flow system, which should be taken into account. To investigate the potential effect of viscosity ratio towards the critical capillary number, a series of simulations are conducted with a wide range of viscosity ratio and a chosen static contact angle of 35°. The results are shown in Figure 8. -1.0 -1.2 -1.4 -1.6
logCa
-1.8 -2.0 -2.2 -2.4 -2.6 -2.8 -3.0 -1.0
-0.5
0.0
0.5
1.0
logM
Figure 8 The critical number with different viscosity ratio (35°) Figure 8 demonstrates that the viscosity ratio can strongly affect the critical capillary number, and the critical capillary number increases with the increasing in the viscosity ratio for a constant static contact angle. Figure 9 shows the critical capillary numbers as a function of viscosity ratio with various static contact angles. It can be observed that the critical capillary number increases with the increasing of static contact angle, this phenomenon can be explained by Figure 3, in which the difference between dynamic contact angle and static contact angle decreases with the increasing of static contact angle at a constant capillary number. -1.0 -1.5
-2.0
-2.5
logCa
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Langmuir
-3.0 SCA=30° SCA=32.5° SCA=35°
-3.5
-4.0
-4.5 -1.0
-0.5
0.0
0.5
logM
ACS Paragon Plus Environment
1.0
Langmuir
Figure 9 The capillary numbers as a function of viscosity with various static contact angle Figure 8 and Figure 9 demonstrate that the critical capillary number is determined by the initial static contact angle of material and the viscosity ratio of immiscible two phases. Therefore, the capillary number and the static contact angle is not enough to judge whether a dynamic contact angle model should be taken into account. A three dimensional universal chart of critical capillary number, which functions as static contact angle and viscosity ratio should be given based on about 500 group simulations. The chart is shown in Figure 10.
logCa
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 14 of 18
Figure 10 The critical capillary number versus viscosity ratio and static contact angle For ease of use, a fitting formula based on the data in Figure 10 is provided in this study with R2=0.9734.
( log Ca )critical
= −2.872 ×10−2 θ S lg M − 2.207θ S 2 + 1.784 lg M + 0.2559θ S − 8.627
(21)
For a particular flow condition of oil slug in a capillary tube, the static contact angle θs, the viscosity ratio M (including the viscosity of each phase µw and µo), surface tension σ and characteristic flow velocity Uin are known. Then the critical capillary number can be calculated by using the fitting formula. If the calculated capillary number of flow system smaller than the critical capillary number, then it is unnecessary to incorporate the dynamic contact angle model into the simulation. Otherwise, we should employ an appropriate dynamic contact model to approach the realistic results.
However, this research is incomplete for a few reasons. Firstly, the upper limit of logarithm of capillary number in this study is -1, which means the 3D chart presented here is suitable for oil slug flow in reservoirs, but cannot extend to the other industrial application directly. Secondly, the range of static contact angle provided is limited and it could be extend based on more simulations. Thirdly, the dynamic contact models used in the simulation are hydrodynamic and empirical models, the molecular kinetic theory based model should be incorporate into our study despite the difficulty in determining parameters.
Conclusions In this work, a direct numerical simulation method together with phase field method are employed to simulate the oil slug flow in a capillary tube at the pore scale. The critical capillary numbers as functions of static contact angle and viscosity ratio of two phase are provided to judge whether a dynamic contact angle model needs to be incorporated into the oil slug simulation. Also,
ACS Paragon Plus Environment
Page 15 of 18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Langmuir
a fitting formula of the logarithm of critical capillary number is presented for ease of use.
Author information Corresponding authors *E-mail:
[email protected] *E-mail:
[email protected] Author contribution Jun Yao and Lei Zhang contributed equally to this work. Notes The authors declare no competing financial interest.
Acknowledgements This work was financially supported by the National Natural Science Foundation of China (51490654, 51504276), the China Postdoctoral Science Foundation (2015M580621) and National Science and Technology Major Project (2016ZX05011-001). References 1.
Dai, L.; Zhang, Y. Experimental study of oil–water two-phase flow in a capillary model. J. Pet.
Sci. Eng. 2013, 108, 96-106. 2.
Hou, J. Network modeling of residual oil displacement after polymer flooding. J. Pet. Sci. Eng.
2007, 59, 321-332. 3.
Dong, M.; Fan, Q.; Dai, L. An experimental study of mobilization and creeping flow of oil slugs
in a water-filled capillary. Transp. Porous Media 2009, 80, 455-467. 4.
Ahmed, W. H. Experimental investigation of air–oil slug flow using capacitance probes, hot-film
anemometer, and image processing. Int. J. Multiphase Flow 2011, 37 (8), 876-887. 5.
Zhu, G.-p.; Yao, J.; Sun, H.; Zhang, M.; Xie, M.-j.; Sun, Z.-x.; Lu, T. The numerical simulation of
thermal recovery based on hydraulic fracture heating technology in shale gas reservoir. J. Nat. Gas Sci. Eng. 2016, 28, 305-316. 6.
Heshmati, M.; Piri, M. Experimental investigation of dynamic contact angle and capillary rise in
tubes with circular and noncircular cross sections. Langmuir 2014, 30, 14151-14162. 7.
Blake, T.; Shikhmurzaev, Y. Dynamic wetting by liquids of different viscosity. J. Colloid Interface
Sci. 2002, 253, 196-202. 8.
Kaveh, N. S.; Rudolph, E.; Van Hemert, P.; Rossen, W.; Wolf, K.-H. Wettability evaluation of a
CO2/water/bentheimer sandstone system: Contact angle, dissolution, and bubble size. Energy Fuels 2014, 28, 4002-4020. 9.
Blake, T. D. The physics of moving wetting lines. J. Colloid Interface Sci. 2006, 299, 1-13.
10. Jadhunandan, P.; Morrow, N. R. Effect of wettability on waterflood recovery for crude-oil/brine/rock systems. SPE Reservoir Eng. 1995, 10, 40-46. 11. Ramé, E.; Garoff, S. On identifying the appropriate boundary conditions at a moving contact line: an experimental investigation. J. Fluid Mech. 1991, 230, 97-116. 12. Quéré, D. Wetting and roughness. Annu. Rev. Mater. Res. 2008, 38, 71-99. 13. Chow, T. Wetting of rough surfaces. J. Phys.: Condens. Matter 1998, 10, L445. 14. Wang, X.; Lee, D.; Peng, X.; Lai, J. Spreading dynamics and dynamic contact angle of non-Newtonian fluids. Langmuir 2007, 23, 8042-8047.
ACS Paragon Plus Environment
Langmuir
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
15. Popescu, M. N.; Ralston, J.; Sedev, R. Capillary rise with velocity-dependent dynamic contact angle. Langmuir 2008, 24, 12710-12716. 16. Slovin, M. R.; Shirts, M. R. Identifying Differences and Similarities in Static and Dynamic Contact Angles between Nanoscale and Microscale Textured Surfaces Using Molecular Dynamics Simulations. Langmuir 2015, 31 (29), 7980-7990. 17. Hansen, R. J.; Toong, T. Dynamic contact angle and its relationship to forces of hydrodynamic origin. J. Colloid Interface Sci. 1971, 37, 196-207. 18. Šikalo, Š.; Wilhelm, H.-D.; Roisman, I.; Jakirlić, S.; Tropea, C. Dynamic contact angle of spreading droplets: experiments and simulations. Phys. Fluids 2005, 17, 062103. 19. Huh, C.; Mason, S. The steady movement of a liquid meniscus in a capillary tube. J. Fluid Mech. 1977, 81, 401-419. 20. De Ruijter, M. J.; Blake, T.; De Coninck, J. Dynamic wetting studied by molecular modeling simulations of droplet spreading. Langmuir 1999, 15, 7836-7847. 21. Blake, T.; Haynes, J. Kinetics of liquidliquid displacement. J. Colloid Interface Sci. 1969, 30, 421-423. 22. Bracke, M.; De Voeght, F.; Joos, P. The kinetics of wetting: the dynamic contact angle. Prog. Colloid Polym. Sci. 1989, 79, 142. 23. Cox, R. The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow. J. Fluid Mech. 1986, 168, 169-194. 24. Saha, A. A.; Mitra, S. K. Effect of dynamic contact angle in a volume of fluid (VOF) model for a microfluidic capillary flow. J. Colloid Interface Sci. 2009, 339, 461-480. 25. Fries, N.; Dreyer, M. The transition from inertial to viscous flow in capillary rise. J. Colloid Interface Sci. 2008, 327, 125-128. 26. Zhao, X.; Blunt, M. J.; Yao, J. Pore-scale modeling: Effects of wettability on waterflood oil recovery. J. Pet. Sci. Eng. 2010, 71, 169-178. 27. Auset, M.; Keller, A. A. Pore‐scale processes that control dispersion of colloids in saturated porous media. Water Resour. Res. 2004, 40. 28. Yao, J.; Sun, H.; Fan, D.-y.; Wang, C.-c.; Sun, Z.-x. Numerical simulation of gas transport mechanisms in tight shale gas reservoirs. Pet. Sci. 2013, 10, 528-537. 29. Zhang, L.; Kang, Q.; Yao, J.; Gao, Y.; Sun, Z.; Liu, H.; Valocchi, A. J. Pore scale simulation of liquid and gas two-phase flow based on digital core technology. Sci. China: Technol. Sci. 2015, 58, 1375-1384. 30. Carlo, M. Review of Direct Methods to Simulate Multiphase Flow. 2010. 31. Chen, L.; Luan, H.-B.; He, Y.-L.; Tao, W.-Q. Pore-scale flow and mass transport in gas diffusion layer of proton exchange membrane fuel cell with interdigitated flow fields. Int. J. Therm. Sci. 2012, 51, 132-144. 32. Tryggvason, G.; Bunner, B.; Esmaeeli, A.; Juric, D.; Al-Rawahi, N.; Tauber, W.; Han, J.; Nas, S.; Jan, Y.-J. A front-tracking method for the computations of multiphase flow. J. Comput. Phys. 2001, 169 (2), 708-759. 33. Peszynska, M.; Trykozko, A. Pore-to-core simulations of flow with large velocities using continuum models and imaging data. Comput. Geosci. 2013, 17, 623-645. 34. Gunde, A. C.; Bera, B.; Mitra, S. K. Investigation of water and CO 2 (carbon dioxide) flooding using micro-CT (micro-computed tomography) images of Berea sandstone core using finite element simulations. Energy 2010, 35, 5209-5216.
ACS Paragon Plus Environment
Page 16 of 18
Page 17 of 18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Langmuir
35. Liovic, P.; Liow, J.-L.; Rudman, M. A volume of fluid (VOF) method for the simulation of metallurgical flows. ISIJ Int. 2001, 41, 225-233. 36. Boyer, F.; Lapuerta, C.; Minjeaud, S.; Piar, B.; Quintard, M. Cahn–Hilliard/Navier–Stokes model for the simulation of three-phase flows. Transp. Porous Media 2010, 82, 463-483. 37. Yue, P.; Feng, J. J.; Liu, C.; Shen, J. A diffuse-interface method for simulating two-phase flows of complex fluids. J. Fluid Mech. 2004, 515, 293-317. 38. Yue, P.; Zhou, C.; Feng, J. J.; Ollivier-Gooch, C. F.; Hu, H. H. Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing. J. Comput. Phys. 2006, 219, 47-67. 39. Diotallevi, F.; Biferale, L.; Chibbaro, S.; Lamura, A.; Pontrelli, G.; Sbragaglia, M.; Succi, S.; Toschi, F. Capillary filling using lattice Boltzmann equations: the case of multi-phase flows. Eur. Phys. J. Special Top. 2009, 166, 111-116. 40. Villanueva, W.; Amberg, G. Some generic capillary-driven flows. Int. J. Multiphase Flow 2006, 32, 1072-1086. 41. Hocking, L.; Rivers, A. The spreading of a drop by capillary action. J. Fluid Mech. 1982, 121, 425-442. 42. Moffatt, H. Viscous and resistive eddies near a sharp corner. J. Fluid Mech. 1964, 18, 1-18. 43. Hoffman, R. L. A study of the advancing interface. I. Interface shape in liquid—gas systems. J. Colloid Interface Sci. 1975, 50, 228-241. 44. Brochard-Wyart, F.; De Gennes, P. Dynamics of partial wetting. Adv. Colloid Interface Sci. 1992, 39, 1-11. 45. Liu, H.; Valocchi, A. J.; Werth, C.; Kang, Q.; Oostrom, M. Pore-scale simulation of liquid CO 2 displacement of water using a two-phase lattice Boltzmann model. Adv. Water Res. 2014, 73, 144-158. 46. Xiao, Y.; Yang, F.; Pitchumani, R. A generalized analysis of capillary flows in channels. J. Colloid Interface Sci. 2006, 298, 880-888.
ACS Paragon Plus Environment
Langmuir
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
TOC graphic
ACS Paragon Plus Environment
Page 18 of 18
