Mathematical Problems in Engineering
Volume 2017 (2017), Article ID 6373957, 12 pages
https://doi.org/10.1155/2017/6373957

Research Article

Deformation Behavior between Hydraulic and Natural Fractures Using Fully Coupled Hydromechanical Model with XFEM

Fei Liu,¹ Zhifeng Luo,² Yu Sang,¹ Liqiang Zhao,² and Changlin Zhou¹

¹Engineering Technology Research Institute of Southwest Oil & Gas Field Company, PetroChina, Chengdu 610017, China
²State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu 610500, China

Correspondence should be addressed to Fei Liu; liufei10628@aliyun.com

Received 12 March 2017; Revised 20 April 2017; Accepted 7 May 2017; Published 7 June 2017

Academic Editor: Eric Feulvarch

Copyright © 2017 Fei Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

There has been a growing consensus that preexisting natural fractures play an important role during stimulation. A novel fully coupled hydromechanical model using extended finite element method is proposed. This directly coupled scheme avoids the cumbersome process during calculating the fluid pressure in complicated fracture networks and translating into an equivalent nodal force. Numerical examples are presented to simulate the hydraulic fracture propagation paths for simultaneous multifracture treatments with properly using the stress shadow effects for horizontal wells and to reveal the deformation response and interaction mechanism between hydraulic induced fracture and nonintersected natural fractures at orthotropic and nonorthotropic angles. With the stress shadow effects, the induced hydraulic flexural fracture deflecting to wellbore rather than transverse fracture would be formed during the progress of simultaneous fracturing for a horizontal well. The coupled hydromechanical simulation reveals that the adjacent section to the intersection is opened and the others are closed for orthogonal natural fracture, while the nonorthogonal natural fracture is activated near the intersection firstly and along the whole section with increasing perturbed stresses. The results imply that the induced hydraulic fracture tends to cross orthotropic natural fracture, while it is prior to being arrested by the nonorthotropic natural fracture.

1. Introduction

Unconventional reservoirs, such as shale gas, have an extremely tight matrix, but well-developed natural fractures [1]. Volume stimulation drives these natural fractures to expand continuously and introduce shear slip in brittle rocks, thus forming complicated networks of intersecting natural fractures and artificial fractures that further increase the stimulated volume [2]. Therefore, the volume stimulation method has become one of the key technologies in unconventional oil and gas development [3]. However, simulation hydraulic fracturing in the presence of natural fracture network is a challenging task, owing to the complex interactions between fracturing fluid, rock matrix, and rock interfaces, as well as the mechanical behaviors before, during, and after the intersections between hydraulic fractures and preexisting natural fractures [4–6].

The morphology of complicated fracture networks created by volume fracturing can be estimated by determining the mechanical mechanism and propagation behavior of hydraulically induced fractures before, during, and after intersecting with natural fractures. In laboratory experiments, a true triaxial fracturing simulation system can be used in combination with high-energy CT scanning, acoustic emission, and microscopic X-ray chromatography to characterize the intersection behavior between hydraulically induced fractures and natural fractures during the fracturing of artificial rock or reservoir outcrop rock samples. These methods also reveal the influences of in situ stress state, horizontal principle stress difference, approach angle, natural fracture interface properties, injection fluid amount, and fracturing fluid viscosity [7–11]. Scholars have developed a series of criteria to describe the intersection behavior between hydraulic fractures and natural fractures. For example, Renshaw and Pollard proposed a criterion for the intersection between hydraulic fractures and orthogonal fractures (R-P criterion) [12]. Gu and Weng expanded the R-P criterion to encompass nonorthogonal conditions (G-W criterion) by considering that the approach angles of hydraulic fractures and natural fractures are between 0° and 90° [13]. Meanwhile, Chen et al. considered the dip angle of natural fractures and expanded the R-P criterion to three-dimensional conditions [14]. Dahi-Taleghani and Olson investigated the tensile property and shearing behavior of natural fractures while hydraulic fractures are approaching but before they come into contact [15]. All of these criteria only consider the mechanical behavior caused by the interaction between hydraulic fractures and natural fractures, without accounting for the influence of the hydraulic effect induced by the injection of fracturing fluid.

Volume fracturing presents a complicated hydromechanical coupling problem, with two special conditions. Firstly, there would be fluid flow within the fracture network and matrix rock and fluid exchange or fluid leak between the fracture and the surrounding porous media, because of the pressure difference between the fracture network and the matrix rock during fracturing fluid pumping. Secondly, although the net pressure drives the initiation and propagation of fractures, the increased pore fluid pressure in the matrix surrounding the fractures would impact fracture propagation. There were considerable differences between piston-like displacement, single phase and multiphase leakage behaviors for nonintersected fractures, and matrix rock using extended finite element method (XFEM) [16]. Based on the vital important role of leakage, this paper promoted the multiphase fluid flow model to fully coupled hydromechanical model for fractured porous media to provide a numerical tool for simulating hydraulic fracturing of fractured porous media. This decision is enlightened on the strong discontinuity of displacement at either side of fractures and the weak discontinuity of pore fluid pressure, as XFEM represents the most commonly used method in simulating discontinuous media. This paper also discusses the deformation features and mechanical mechanisms before the intersection between hydraulically induced fractures and orthogonal and nonorthogonal natural fractures.

2. Fully Coupled Hydromechanical Model

The fractured porous media were composed of matrix rock and macro fractures as a whole. Although there were plenty of micro fractures, which show performances similar to those of macro fractures, developed in matrix rock, the fractured porous media were considered to be constructed by an intact matrix and macro fracture networks for simplicity, as displayed in Figure 1.

Figure 1: Schematic diagram of hydromechanical coupling for fractured porous media.

The strong form, weak form, and the associated discretization form of governing equations for hydromechanical coupling inside both matrix domain and fracture domain are demonstrated in this section.

2.1. Strong Form

Figure 2 shows a bounded control body () with an outer boundary of .

Figure 2: Control body and boundary conditions of fractured porous media.

The effective stress equilibrium equation for porous media is written as

According to the law of mass conservation, the continuity equations for the wetting and nonwetting phase fluid flow in the matrix domain are expressed as

Similarly, the continuity equations for the wetting and nonwetting phase fluid flow inside the fracture domain are given by

The displacement and pressure of the wetting and nonwetting phase initially are provided according to the initial conditions

The displacement and pressure of the wetting and nonwetting phase at the boundary are provided by the forced boundary conditions

The stress and mass flow rate of the wetting and nonwetting phase at the boundary are provided by the natural boundary conditions

The fluid pressure is applied on the fracture surface

The pore fluid pressure at the fracture boundary is continuous, with fluid exchange occurring between fractures and the surrounding matrix:

2.2. Weak Form

The weak form of the governing equations can be obtained by combining the weighted residual method with the associated boundary conditions. With regard to fluid flow in fractures, since the fracture width is much smaller than its length, the fluid pressure variation along its width can be neglected (see Figure 3). The areal integral inside the fracture domain can be transformed into the linear integral along the fracture surface . The fluid exchange items or leak between the fractures and the surrounding matrix rock can then be vanished by adding the equivalent integrals of the matrix rock and fractures, which have the same absolute value but opposite signs. Thus, it becomes unnecessary to explicitly express the fluid exchange items or leak [16–18]:

Figure 3: Schematic diagram of fluid flows in a fracture.

2.3. Discretization Form

XFEM introduces Heaviside step functions and branch functions to enrich the shape function, which make it easy to deal with typical strong discontinuous problems like fractures without grid matching and remeshing. The XFEM displacement field for multiple fractures is [19]

In fractured porous media, the natural fractures are the main channels for seepage. The pore fluid pressure field is continuous at the fracture surface, while the derivative (i.e., flow rate) is discontinuous. This presents a typical weak discontinuous problem. Similar to the interpolation of displacement field, the pressure field can be interpolated as [16, 18]

The discretization form of the fully coupled hydromechanical model for fractured porous media can be obtained by regarding the trial functions and in (9)–(11) as the variation of displacement and pressure interpolation functions, as shown in

The stiffness matrix , coupling matrix , flow matrix , compressibility matrix , saturation matrix , equivalent nodal force vector , and equivalent node flow rate vector are assembled matrixes for each element matrix labeled by the degrees of freedom, where each element matrix is expressed as

Equation (14) is implicitly solved in the time domain using the generalized trapezoidal method. The implicit pressure explicit saturation method (IMPES) is used to treat the multiphase fluid flow. This strong nonlinear multiphase fluid hydromechanical coupling problem can be solved by incorporating the Newmark scheme [20].

3. Numerical Simulation

3.1. Simultaneous Fracturing for Horizontal Wells

The cluster perforation and simultaneous fracturing are key technologies for volume stimulation of horizontal wells. The propped hydraulic fracture induces additional stresses in the surrounding rock. This stress interference between fractures can complicate the propagation of hydraulic fractures to form complex fracture networks. Crosby et al. [21] and Rahman et al. [22] investigated the propagation behavior of two isolated hydraulic fractures induced by the simultaneous fracturing of a horizontal well.

This study simulated Crosby’s DC-8 specimen, as Figure 4 shows. Modeling was performed using material properties and injection parameters similar to those employed during the laboratory test DC-8. The specimen was meshed into uniform 80 × 80 structured bilinear quadrilateral elements, with the propagation step set at a constant value of 0.01 m. The coupling process was neglected during this simulation because the synthetic fracture test material was cast by silica flour and acrylic emulsion to provide the specimen with a negligible leak coefficient, 3.8 × 10⁻⁸ , and an ultralow permeability, 8 × 10⁻⁵ mD. Therefore, only the first item in (14) was used. Equation (9) states that the fluid pressure inside the fracture, , should be calculated, which presents a very complicated problem, especially for fractures that are not straight. Moreover, the fluid pressure inside the fracture must be transferred into an equivalent nodal force using the linear integral form, which requires integral points on the fracture surface [23]. We set the net pressure inside the fractures at 3 MPa.

Figure 4: Dimensions of the artificial rock specimen.

The maximum circumferential stress criterion is selected to determine the propagation condition and crack direction. The stress intensity factors of hydraulic fracture were calculated with the interaction integral

The first item of the right side can be calculated easily. The second item of the right side was derived explicitly:

Figure 5 shows comparisons between the XFEM simulation results, experimental results, and HYFRANC finite element simulation results.

Figure 5: Comparison between the experimental results and XFEM numerical simulation results of the staged fracturing of a horizontal well.

The results of XFEM simulation are more consistent with the experimental result in the upper fracture than in the lower fracture. This can be largely attributed to the fact that the lower fracture was initiated preferentially during the experiment. Therefore, the propagation path was relatively straight. However, the two fractures were considered to initiate at the same time in the simulation. Then, the lower fracture was forced to deflect towards the wellbore for stress shadow effect. In comparison, the HYFRANC simulation results [21] and the XFEM simulation results show high consistency with the experimental results. Moreover, the simulation efficiency improved significantly.

The stress distribution at the sixth propagation step (Figure 6) shows that the -axis stress, , between the two fractures increased from 5 MPa to about 7.62 MPa (compressive), the -axis stress, , increased from 7 MPa to about 8.06 MPa (compressive), and the shear stress, , was lower than 0.0002 MPa. In other words, the stress state between the two fractures changed significantly. The maximum horizontal nominal stress direction changed from the -axis to the -axis. Perforation fracturing between two fractures tends to cause the hydraulic fracture to propagate along the -axis, intersect with other fractures, and form a fracture network.

(a)-axis stress (MPa)

(b)-axis stress (MPa)

(c) planar shearing stress (0.1 kPa)

Figure 6: The displacement distribution of the sixth propagation step.

3.2. Coupled Hydromechanical Behavior for Hydraulic Fracturing

Figure 7 shows a 2D fractured porous media control domain. The physical model has a length and width of 150 m, which contains three fractures. Table 1 lists detailed parameters for the model. The hydraulic fracture, which is parallel to the global -axis, has a length of 50 m and an initial width of 1 mm. The orthogonal natural fracture, which is parallel to the global -axis, has a length of 40 m and an initial width of 0.1 mm. The nonorthogonal natural fracture, which lies at an included angle of 45° to the global -axis, has a length of 30 m and an initial width of 0.1 mm. The initial stress is zero, while the pore fluid pressure inside both the fractures and the matrix rock is 50 MPa. The displacements of the upper, lower, right, and left sides are all restricted (), and they are therefore impermeable boundaries. In reality, the length, width, height, and net pressure would be increased with increasing injection time, and the stress intensity factors and induced stresses would be strengthened to force the deformation of natural fracture, subsequently. For simplicity, the fractures were presumed as nonpropagating cracks to investigate the deformation behavior between hydraulic and natural fractures before intersection. However, the net pressure and width of fractures were solved with the fully coupled hydromechanical model. Although this assumption was not coinciding with the reality, the numerical simulation can represent the deformation behavior between hydraulic and natural fractures to some extent.

Table 1: Model parameters for the seepage-stress coupling analysis.

Item
Value

Young’s modulus, (GPa)
30

Poisson’s ratio,
0.25

Biot coefficient,
0.85

Rock permeability, (mD)
0.1

Rock porosity,
0.08

Bulk modulus of rock, (GPa)
134.6

Coefficient between real and ideal fractures,
1.3

Bulk modulus of wetting phase, (GPa)
2.2

Viscosity of wetting phase, (mPa⋅s)
1

Liquid injection rate, (m³/min)
10

Reservoir thickness, (m)
30

Density of wetting phase, (kg/m³)
1000

Volume factor of wetting phase,
1.05

Figure 7: Control domain of fractured porous media for the hydromechanical coupling analysis.

A uniform 80 × 79 bilinear quadrilateral mesh was used for calculation, with the intention of locating the injection well point at the element boundary, while the fracture location was independent of the generated mesh.

The -axis displacement, , distribution is shown in Figure 8. The -axis displacement, , distribution is shown in Figure 9. The displacement fields demonstrate that the injected fluid initially flows along the first fracture, forcing the fluid pressure inside the fracture to increase. Meanwhile, the fluid filtrated into the surrounding matrix rock, leading to an increase of pore fluid pressure in the matrix rock. The increase in fluid pressure drives the deformation of fractured porous media, generating displacement along the -axis and the -axis. Since the contact condition is not considered here, parts of the other fracture surfaces may come into contact and intersect with each other, while other parts become staggered.

(a) = 60 s

(b) = 600 s

Figure 8: -axis displacement at different moments (mm).

(a) = 60 s

(b) = 600 s

Figure 9: -axis displacement at different moments (mm).

Fracture width represents the key parameter in hydraulic fractures. Equation (12) allows the displacement of the two surfaces of a fracture at any position to be obtained. The fracture width can be calculated by subtracting the positive-side displacement from the negative-side displacement. The width decreases with increasing distance from the injecting well point, and it increases with increasing injection duration, as Figure 10 shows. The fracture width is distributed symmetrically. The orthogonal natural fracture remains closed at first under the effect of stress induced by the left crack tip of a hydraulic fracture. As the fluid pressure inside the hydraulic fracture increases, the crack-tip-induced stress increases and forces the center section to open, while the faraway section remains closed, as shown in Figure 11. This is one of the reasons why a hydraulic fracture intersects a natural fracture at a large intersection angle. In contrast, the stress induced by a hydraulic fracture drives the nonorthogonal natural fracture to open. As Figure 12 shows, only a small part of the initial fracture remains closed, while the whole fracture opens finally.

Figure 10: Width section of the hydraulic fracture.

Figure 11: Width section of the orthogonal natural fracture.

Figure 12: Width section of the nonorthogonal natural fracture.

4. Conclusion

Hydraulic fracturing is an interactive process incorporating both mechanical and hydraulic effects. The hydromechanical coupling model and associated solution scheme proposed under the framework of XFEM can effectively avoid tedious calculations of both the fluid pressure distribution in fractures and the transformation of equivalent nodal forces. This model also considers the influence of pore pressure variation in the rock matrix on fracture propagation.

The mechanical module can simulate the fracture propagation pathway of the simultaneous fracturing for a horizontal well. Under the influence of stress perturbation between hydraulic fractures, not the straight transverse fractures but the flexural fractures deflecting toward the wellbore were created for the staged fracturing of a horizontal well. The stress shadow effect introduced by the propped hydraulic fracture reduces the horizontal stress differential between fractures or even leads to stress inversion. Therefore, it is beneficial to enhance the complexity of fracture networks.

The induced stresses added to the crack tip of the hydraulic fracture keep most of the orthogonally intersected natural fractures closed, with only a small part opening due to the induced stress and fracturing fluid leakage. The stress concentration at the crack tip of the hydraulic fracture fully opens most of the nonorthogonal intersected natural fractures as time passes. The hydraulic fracture tends to traverse the orthogonal natural fracture and continue to propagate, while it is more easily arrested by the nonorthogonal natural fracture.

Nomenclature

: | Overall stress
: | Effective stress
: | Fluid density
: | Fluid density of wetting phase
: | Fluid density of nonwetting phase
: | Gravity acceleration vector
: | Unit vector
: | Biot coefficient
ϕ: | Porosity
: | Bulk modulus of rock particle
: | Bulk modulus of wetting phase
: | Bulk modulus of nonwetting phase
: | Saturation of wetting phase in matrix rock
: | Saturation of nonwetting phase in matrix rock
: | Saturation of wetting phase in fracture
: | Saturation of nonwetting phase in fracture
: | Pore fluid pressure
: | Capillary pressure of matrix rock
: | Capillary pressure of fracture
: | Wetting phase fluid pressure in matrix rock
: | Nonwetting phase fluid pressure in matrix rock
: | Initial wetting phase fluid pressure in matrix rock
: | Initial nonwetting phase fluid pressure in matrix rock
: | Wetting phase fluid pressure at the boundary
: | Nonwetting phase fluid pressure at the boundary
: | Wetting phase fluid pressure in fracture
: | Nonwetting phase fluid pressure in fracture
: | Permeability tensor of matrix rock
: | Fracture permeability
: | Relative permeability of wetting phase in matrix rock
: | Relative permeability of nonwetting phase in matrix rock
: | Relative permeability of wetting phase in fracture
: | Relative permeability of nonwetting phase in fracture
: | Viscosity of wetting phase
: | Viscosity of nonwetting phase
: | Ground source sink term of wetting phase
: | Bulk coefficient of wetting phase
: | Displacement vector
: | Boundary displacement vector
: | Initial displacement
: | Strain tensor
: | Boundary force tensor
: | Flow rate of wetting phase at the boundary
: | Flow rate of nonwetting phase at the boundary
: | Unit normal vector at the boundary
: | Value difference between two sides of fracture
: | Transformation matrix between local coordinate and global coordinate
: | Unit normal vector in local coordinate system of fracture surface
: | Unit tangential vector in local coordinate system of fracture surface
: | Interpolating shape function
: | Heaviside function
: | Ramp function
: | Branch function for displacement
: | Junction function
: | Signed distance function
: | Branch function of pore pressure
: | Interpolating function matrix of displacement
: | Interpolating function matrix of pore pressure
K: | Stiffness matrix
L: | Coupling matrix
C: | Compressibility matrix
S: | Saturation matrix
B: | Strain matrix
: | Equivalent nodal force vector
: | Elasticity matrix of rock
: | Equivalent nodal flow rate vector
: | Fracture width.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors are grateful for the research support of the National Natural Science Foundation of China (NSFC, no. 51404207) and the Postdoctoral Program of Postdoctoral Work Centre, Southwest Oil & Gas Field Company (no. 20150304-08).

References

J. F. W. Gale, R. M. Reed, and J. Holder, “Natural fractures in the Barnett Shale and their importance for hydraulic fracture treatments,” AAPG Bulletin, vol. 91, no. 4, pp. 603–622, 2007.

Z. Chen, C. Xue, T. Jiang et al., “Proposals for the application of fracturing by stimulated reservoir volume (SRV) in shale gas wells in China,” Natural Gas Industry, vol. 30, no. 10, pp. 30–32, 2010.

Q. Wu, Y. Xu, X. Wang et al., “Volume fracturing technology of unconventional reservoirs: Connotation, optimization design and implementation,” Petroleum Exploration and Development, vol. 39, no. 3, pp. 352–358, 2012.

P. Fu, S. M. Johnson, and C. R. Carrigan, “An explicitly coupled hydro-geomechanical model for simulating hydraulic fracturing in arbitrary discrete fracture networks,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 37, no. 14, pp. 2278–2300, 2013.

L. Zhao, F. Liu, P. Wang et al., “A review of creation and propagation of complex hydraulic fracture network,” Oil and Gas Geology, vol. 35, no. 4, pp. 562–569, 2014.

M. Sheng and G. Li, “Extended finite element modeling of hydraulic fracture propagation,” Engineering Mechanics, vol. 31, no. 10, pp. 123–128, 2014.

T. Guo, S. Zhang, Z. Qu, T. Zhou, Y. Xiao, and J. Gao, “Experimental study of hydraulic fracturing for shale by stimulated reservoir volume,” Fuel, vol. 128, pp. 373–380, 2014.

J. Zhou, M. Chen, Y. Jin, and G.-Q. Zhang, “Analysis of fracture propagation behavior and fracture geometry using a tri-axial fracturing system in naturally fractured reservoirs,” International Journal of Rock Mechanics and Mining Sciences, vol. 45, no. 7, pp. 1143–1152, 2008.

M. Chen, J. Zhou, Y. Jin et al., “Experimental study on fracturing features in naturally fractured reservoir,” Acta Petrolei Sinica, vol. 29, no. 3, pp. 431–434, 2008.

J. Zhou, M. Chen, Y. Jin et al., “Experiment of propagation mechanism of hydraulic fracture in multi-fracture reservoir,” Journal of China University of Petroleum, vol. 32, no. 4, pp. 51–59, 2008.

S. Mohammad, Experimental and numerical study of interaction of a pre-existing natural interface and an induced hydraulic fracture [Ph.D. thesis], Ph.D. Curtin University, 2012.

C. E. Renshaw and D. D. Pollard, “An experimentally verified criterion for propagation across unbounded frictional interfaces in brittle, linear elastic materials,” International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, vol. 32, no. 3, pp. 237–249, 1995.

H. Gu, X. Weng, and J. Lund, “Hydraulic fracture crossing natural fracture at non-orthogonal angles: a criterion and its validation and application,” SPE Production & Operations, vol. 27, no. 1, pp. 20–26, 2012.

W. Chen, Y. Jin, M. Chen et al., “A criterion for identifying hydraulic fractures crossing natural fractures in 3D space,” Petroleum Exploration and Development, vol. 41, no. 3, pp. 336–340, 2014.

A. Dahi-Taleghani and J. E. Olson, “Numerical modeling of multistranded-hydraulic-fracture propagation: accounting for the interaction between induced and natural fractures,” SPE Journal, vol. 16, no. 3, pp. 575–581, 2011.

P. Liu, F. Liu, C. She et al., “Multi-phase fracturing fluid leakoff model for fractured reservoir using extended finite element method,” Journal of Natural Gas Science and Engineering, vol. 28, pp. 548–557, 2016.

T. Mohammadnejad and A. R. Khoei, “An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model,” Finite Elements in Analysis and Design, vol. 73, pp. 77–95, 2013.

F. Liu, L.-q. Zhao, P.-l. Liu, Z.-f. Luo, N.-y. Li, and P.-s. Wang, “An extended finite element model for fluid flow in fractured porous media,” Mathematical Problems in Engineering, Article ID 604212, Art. ID 604212, 10 pages, 2015.

T. T. Yu, Extended Finite Element Method: Theory, Application and Program, Science Press, Beijing, China, 2014.

R. W. Lewis and B. A. Schrefler, The Finite Element Mothod in the Static And Dynamic Deformation and Consolidation of Porous Media, John Wiley & Sons, New York, NY, USA, 1998.

D. G. Crosby, M. M. Rahman, M. K. Rahman, and S. S. Rahman, “Single and multiple transverse fracture initiation from horizontal wells,” Journal of Petroleum Science and Engineering, vol. 35, no. 3-4, pp. 191–204, 2002.

M. M. Rahman, M. M. Hossain, D. G. Crosby, M. K. Rahman, and S. S. Rahman, “Analytical, numerical and experimental investigations of transverse fracture propagation from horizontal wells,” Journal of Petroleum Science and Engineering, vol. 35, no. 3-4, pp. 127–150, 2002.

C. Annavarapu, R. R. Settgast, S. M. Johnson, P. Fu, and E. B. Herbold, “A weighted Nitsche stabilized method for small-sliding contact on frictional surfaces,” Computer Methods in Applied Mechanics and Engineering, vol. 283, pp. 763–781, 2015.