1 Introduction
Reservoir simuations are important tools for petroleum engineering, which have been applied to model underground flows and interations between reservoirs and wells to predict the well performance such as oil rates, water rates and bottom hole pressure. When a reservoir model is large enough, it may take severals days or even longer for the simulator to complete a single model. Therefore, the numerical methods and parallel techniques are essential to accelerate the simulations.
As a type of unconventional reservoirs, the reservoir with natural fractures and hydraulic fractures which are commonly seen in tight and shale oil and gas reservoir, the dualporosity/dualpermeability model and the multiple iteration contiua (MINC) model[19, 30] homogenize the fractures and use superpositioned cells to represent the fractures and the matrix. The multiple continuum approaches have been successfully employed in the unconventional reservoir problems. Wu et al.[28, 29] regarded the fractured vuggy rocks as a triple or multiplecontinuum medium with highly permeability and wellconnected fractures, lowpermeabilty rock matrix and varioussized vugs. Wu et al. also developed a hybrid multiplecontinuummedium modeling appoach to describe different types of fractures including hydrautic fractures, natural fracture network and micro fractures [27]. Jiang and Moinfar and their collaborators designed explicit fracture models coupled with MINC model to simulate the unconventional reservoirs [17, 10]. In reference [31], the authors presented a unified framework model to incorporate the known mechanisms and processes including gas adsorption and desorption, geomechanics effect, Klinkenberg or gasslippage effect and nonDarcy flow. They also used a hybrid fracturemodeling appoach to simulate the unconventional gas reservoirs. Jiang and Younis developed two alternative hybrid approaches to capture the effects of the multiscale fracture systems by combining the advantages of the multicontinuum and discretefracture/matrix representations [11].
The numerical methods for reservoir simulations have been explored by both mathematicians and engineers for decades. Killough et al. [12] studied local refinement techniques to improve the accuracy and reduce computational cost comparing to the global grid refinement. Those local refinement techniques are useful for the complex models such as insitu combustion. Dogru and his group [7] developed parallel simulators using structured and unstructured grids to handle faults, pinchouts, complex wells, polymer flooding in nonfractured and fractured reservoirs. Zhang et al. developed a generalpurpose parallel platform for largescale scientific applications which was designed using the adaptive element methods and the adaptive finite volumn methods [32, 33]. It has been applied to black oil simulation using discontinuous Galerkin methods[21]. Chen et al. studied finite element emtjods and finite different methods for black oil, compositional and thermal models [6]. They studied Newton methods, linear solvers and preconditioners as well. Chen and his collaborators also developed a parallel platform to support the newgeneration simulator development, such as black oil, compositional, thermal, polymer flooding models [16, 22, 23]. Wheeler studied discretization methods, linear solvers, preconditioner techniques and developed a parallel black oil simulator [26]. For the linear solver and preconditioners, many preconditioning methods have been proposed and applied to reservoir simulations such as constrained pressure residual (CPR) methods [14, 3], multistage methods [1], multiple level preconditioners [24], fast auxiliary space preconditioners (FASP) [9] and a family of parallel CPRlike methods[14].
In this paper, we present the mathematical model of the dualporosity reservoirs. After that, we give the numerical methods used in our simulation inlcuding the nonlinear solvers, the preconditioner and the linear solver and the parallel technique. Finally, we show the numerical results obtained from our simulator.
2 Mathematical Model
Darcy’s law describes the flow of a fluid through porous media, establishing the relationship between the volumtric flow rate and the pressure gradient:
(1) 
where is the permeability of a given reservoir, is the area in the flow direction, is the pressure difference, is the viscosity of the fluid and is the length of the reservoir. In three dimensional space, its differential form is
(2) 
Combining Darcy’s law and the mass conservation law for oil and water components, the twophase model is as follows:
(3) 
where and are the porosity and the permeabilities (in ,  and directions) respectively, is the potential, , , and () are the saturation, viscosity, relative permeability and production/injection rate, respectively, which satisfy the following constraints:
(4) 
where and are the pressure and the mass density of phase (). is the capillary pressure between oil phase and water phase, is the gravitational constant and is the reservoir depth.
For the dualporosity model, fluids can flow between matrix and fracture, as well as fracture and fracture, but there is no flow between matrices. There flow directions are described by Fig 1. The following equations are used to describe this model:
(5) 
where denotes the variables for the fracture and for the matrix. Note that here is a source term which represents the net addition of the fluid to the fracture from the matrix with the shape factor. In our simulator, the shape factor we use is based on the work of Warren and Root[25] and Gilman and Kazemi[8]. Note that denoted the value at or depending on which is the upsteam. The well rate at a perforation, , is modelled by the sinksource method,
(6) 
where is the well index at the perforation, is the bottomhole pressure of the well at the perforation, is the block pressure from oil phase at the perforation. For a producer or an injection well, its well rate is calculated as,
(7) 
The well index is calculated as follows:
(8) 
Here is a given value and is the facture between 0.0 and 1.0 specifying the fraction of a circle that the well models, is the wellbore radius, is a real number specifying the well skin factor, is the well effective radius calculated from
(9) 
with the geometric factor for the well element, the area perpendicular to the wellbore (e.g., for a vertical well) and as same above. The other option for the well effective radius is the following Peaceman’s form [18]
(10) 
where , and are the well effective radius in , and direction respectively, , and are the block sizes in , and direction respectively, , and are the permeabilities in , and direction respectively.
3 Numerical Methods and Parallel Computing
In this paper, the fully implicit method (FIM) is applie, and the oil phase pressure , water saturation and the well bottom hole pressure are unknowns. For the time differential and the space differential term, we use the backward Euler scheme and cellcentered finite difference method as discretization methods.
The system is highly nonlinear,
(11) 
where is a nonlinear map from to with . Here is the number of grid blocks and is the well size. In this nonlinear equation, the properties related to the saturation are strongly nonlinear while the properties related to the pressure are weakly nonlinear. In our simulation, we use Newton method (or inexact Newton method) to solve it.
3.1 Nonlinear Solver
The inexact Newton method can be regarded as an extension of the standard Newton method. Let be the Jacobian matrix of and is the correction between the last step approximate solution and the current step one. Its algorithm [5] is as follows:
(12) 
We can see that the algorithm is the same as Newton method except the choice of the tolerance . The standard Newton method usually applies a small constant such as , while the Algorithm 1 uses adaptive tolerance to avoid over solution. Common methods are defined as follows [5]:
(13) 
where is the residual of the th iteration.
3.2 Preconditioner and Linear Solver
For each Newton iteration, we need to solve a linear system which is very time comsuming. To improve the efficiency, we choose a proper preconditioner, CPRFPF method [14], for the model. Each grid block has four unknowns, pressures ( and ) and saturations ( and ) for matrix and fracture, and four equations. In our linear system, each unknown and each row are arranged cell by cell, in this case, the Jacobian matrix is a block matrix. It is also wellknown that block matrix has better convergence that pointwise matrix.
Before solving the linear system , a decoupling operation is applied, and an equivalent system is obtained:
(14) 
A proper decoupling method can improve linear solver dramatically. Many decoupling strategies have been proposed, such as QuasiIMPES[13] method and ABF[2] method. In this paper, a modified GaussJordan elimination method is applied.
3.3 Parallel Computing
The simulator is based on our parallel platform PRSI [16]
, which is developed using C language and MPI (Message Passing Interface). The platform has implemented many modules, such as grid generation, load balancing, well management, parallel input and output, distributed matrix and vector, linear solver and preconditioner, communication management and visualization. Based on the platform, physical modules, such as rock properties and rockfluid properties, are implemented. More details can be read from reference
[16].4 Numerical Results
Example 1
The grid dimension is with sizes 102.04 ft. in and directions and 100.0 ft. in direction from top to bottom. The depth of the top layer center is 2000.0 ft. The permeabilities for the matrix in , and directions are 100, 100, 10 mD respectively. The permeability for the fracture in , and directions are 395.85 mD. The reference porosity for the matrix is 0.1392. The reference porosity for the fracture is 0.039585. The rock compressibilities for the matrix and fracture are both 3e06 (1/psi). The reference pressure is 15.0 psi for both the matrix and the fracture. Component properties: densities of oil and water are 58.0 lbm/ft3 and 62.4 lbm/ft3 respectively. The reference pressure is 15 psi at which the oil formation volume factor is 1.036 RB/STB and the oil viscosity is 40.0 cp. The oil compressibility is const 1.313e5 l/psi. The initial conditions are as follows: initial pressure for the matrix is 2000 psi, initial pressure for the fracture is 1980 psi, and initial water saturations are 0.08 and 0.01 in matrix and fracture respectively. There are one injection well and two production wells. All of them are vertical. Injection well has maximum water injection rate 500.0 bbl/day, maximum bottom hole pressure 5.0e+4 psi, well index 200.0 with perforation at cell [5 1 1]. Both of the production well has maximum oil production rate 300.0 STB/day, minimum bottom hole pressure 15 psi with well radius 0.25ft. The perforation of Produer 1 is at cell [1 10 1] while the perforation of Produer 2 is at cell [10 10 1]. The simulation period is 800 days.
The results of oil production rate, bottomholepressure and water rate are shown in Fig. 2, 3, 4, from which we can see that the results from our simulator and from CMG IMEX match very well. This proves our methods and implementation are correct.
Example 2
The grid dimension is with mesh size 102 ft. in and directions and 100.0 ft. in direction from top to bottom. The depth of the top layer center is 2000.0 ft. The permeabilities for the matrix in , and directions are 100, 100, 100 mD respectively. The permeability for the fracture in , and directions are 395.85 mD. The reference porosity for the matrix is 0.1392. The reference porosity for the fracture is 0.039585. The rock compressibilities for the matrix and fracture are both 3e06 (1/psi). The reference pressure is 15.0 for both the matrix and the fracture. Component properties: densities of oil and water are 58.0 lbm/ft3 and 62.4 lbm/ft3 respectively. PVT: the reference pressure is 15 psi at which the oil formation volume factor is 1.036 RB/STB and the oil viscosity is 40.0 cp. The oil compressibility is const 1.313e5 l/psi. The initial conditions are as follows: initial pressure for the matrix is 800 psi, initial pressure for the fracture is 500 psi, and initial water saturations are 0.08 and 0.01 in matrix and fracture respectively. There are one injection well and two production wells. All of them are vertical. Injection well has maximum water injection rate 200.0 bbl/day, maximum bottom hole pressure 5.0e+4 psi, well index 200.0 with perforation at cell [5 5 1]. Both of the production well has maximum oil production rate 500.0 STB/day, minimum bottom hole pressure 15 psi with well radius 0.25ft. The perforation of Produer 1 is at cell [1 1 1] while the perforation of Produer 2 is at cell [10 10 1]. The simulation period is 1600 days.
The results of oil production rate, bottomholepressure and water rate are shown in Fig. 5, 6, 7. Again, these figures show that our simulator match CMG simulator.
Example 3
This example tests the scalability of our twophase dual porosity simulator by computing Example 1 with grid dimension .
As we introducerd before, we employ the GMRES linear solver and CPR preconditioner. Table 4 shows the numerical summaries, which show that our numerical methods are stable when increasing CPU cores. Figure 8 presents the scalability of our simulator, which demonstrates that the simulator and parallel implementation are scalable.
#MPIs  8  64 
#Time steps  68  68 
#Newton iterations  273  267 
#total linear iterations  641  623 
#total running time(s)  20009.71  2531.91 
5 Conclusion
This paper presents our work on development of twophase oilwater simulator for natural fractured reservoirs using dual porosity method. Effective numerical methods and parallel computing techniques are introduced. From the numerical experiments, we can see that our results match commercial simulator, CMG IMEX, and the simulatro has good parallel scalability and is capable of handling large scale reservoir simulation problems.
6 Acknowledgements
This work is partially supported by Department of Chemical Petroleum Engineering, University of Calgary, NSERC, AIEES, Foundation CMG, AITF iCore, IBM Thomas J. Watson Research Center, Frank and Sarah Meyer FCMG Collaboration Center, WestGrid (www.westgrid.ca), SciNet (www.scinetpc.ca) and Compute Canada Calcul Canada (www.computecanada.ca).
References
 [1] T. AlShaalan, H. Klie, A. Dogru, and M. Wheeler, Studies of robust two stage preconditioners for the solution of fully implicit multiphase flow problems, SPE Reservoir Simulation Symposium, 2009.

[2]
R. Bank, T. Chan, W. Coughran Jr., and R. Smith, The AlternateBlockFactorization procedure for systems of partial differential equations, BIT Numerical Mathematics, 29(4), 1989, 938954.
 [3] H. Cao, T. Schlumberger, A. Hamdi, J. Wallis, and H. Yardumian, Parallel scalable unstructured CPRtype linear solver for reservoir simulation, SPE Annual Technical Conference and Exhibition, 2005.
 [4] X. Cai and M. Sarkis, A restricted additive Schwarz preconditioner for general sparse linear systems, SIAM Journal on Scientific Computing, 21(2), 1999, 792797.
 [5] T. Chen, N. Gewecke, Z. Li, A. Rubiano, R. Shuttleworth B. Yong and X. Zhong, Fast Comutational Methods for Reservoir Flow Models, 2009.
 [6] Z. Chen, G. Huan and Y. Ma, Computational methods for multiphase flows in porous media, Vol. 2, Siam, 2006.
 [7] A. Dogru, L. Fung, U. Middya, T. AlShaalan, and J. Pita, A nextgeneration parallel reservoir simulation for the giant reservoirs, SPE/EAGE Reservoir Charaterization & Simulation Conference, 2009.
 [8] J.R. Gilman and H. Kazemi, Improvements in Simulation of Naturally Fractured Reservoirs, SPE10511PA, Vol. 23(04), 1983, 695707.
 [9] X. Hu, W. Liu, G. Qin, J. Xu, and Z. Zhang, Development of a fast auxiliary subspace preconditoner for numerical reservoir simulators, SPE Reservoir Charaterisation and Simulation Conference and Exhibition, 2011.
 [10] J. Jiang, Y. Shao, R. M. Younis, et al., Development of a multicontinuum multicomponent model for enhanced gas recovery and CO2 storage in fractured shale gas reservoirs, SPE Improved Oil Recovery Symposium, Society of Petroleum Engineers, 2014.
 [11] J. Jiang, R. Younis, et al., Hybrid coupled discretefracture/matrix and multicontimuum models for unconvenctionalreservoir simulation, SPE Journal 21(03), 2016 10091027.
 [12] J. Killough, D. Camilleri, B. Darlow, and J. Foster, Parallel reservoid simulation based on local grid refinement. SPE37978, SPE Reservoir Simulation Symposium., Dallas, 1997.
 [13] S. Lacroix, Y. Vassilevski, and M. Wheeler, Decoupling preconditioners in the implicit parallel accurate reservoir simulation (IPARS), Numerical linear algebra with applications, 8(8), 2001, 537549.
 [14] H. Liu, K. Wang, and Z. Chen, A family of constrained pressure residual preconditioners for parallel reservoir simulations, Numerical Linear Algebra with Applications, Vol. 23(1), 2016, 120146.
 [15] H. Liu, Dynamic load balancing on adaptive unstructured meshes, 10th IEEE International Conference on High Performace Computing and Communications, 2008.
 [16] H. Liu, K. Wang, Z. Chen, K. Jordan, J. Luo, and H. Deng, A parallel framework for reservoir simulators on distributedmemory supercomputers, SPE17645MS, SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition, 2022 October, Nusa Dua, Bali, Indonesia, 2015.
 [17] A. Moinfar, A. Varavei, K. Sepehrnoori, R. T. Johns, et al., Development of a coupled dual continuum and discrete fracture model for the simulation of unconventional reservoirs, SPE Reservoir Simulation Symposium, Society of Petroleum Engineers, 2013.
 [18] Peaceman, D. W., Interpretation of WellBlock Pressures in Numerical Reservoir Simulation with NonSquare Grid Blocks and Anisotropic Permeability, SPE Journal, 23(3), 1983, 531543.
 [19] K. Pruess, et al. A practical method for modeling fluid and heat flow in fractured porous media, Society of Petroleum Engineers Journal, Vol. 25(01), 1426 1985.
 [20] J. Wallis, R. Kendall, and T. Little, Constrained residual acceleration of conjugate residual methods, SPE Reservoir Simulation Symposium, 1985.
 [21] K. Wang, L. Zhang, Z. Chen, Development of discontinuous Galerkin methods and a parallel simulator for reservoir simulation, SPE176168MS, SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition, 2022 October, Nusa Dua, Bali, Indonesia, 2015.
 [22] K. Wang, H. Liu, and Z. Chen, A scalable parallel black oil simulator on distributed memory parallel computers, Journal of Computational Physics, Vol. 301, 1934.
 [23] K. Wang, H. Liu, J. Luo, and Z. Chen, Parallel simulation of fullfield polymer flooding, The 2nd IEEE International Conference on High Performance and Smart Computing, 2016.
 [24] B. Wang, S. Wu, Q. Li, H. Li, C. Zhang, and J. Xu, A multilevel preconditioner and its shared memory implementation for new generation reservoir simulation, SPE172988MS, SPE Large Scale Computing and Big Data Challenges in Reservoir Simulation Conference and Exibition, 1517 September, Istanbul, Turkey, 2014.
 [25] J.E. Warren and P.J. Root, The Behavior of Naturally Fractured Reservoirs, SPE426PA, Vol. 3(02), 1963. 245255.
 [26] M. Wheeler, Advanced techniques and algorithms for reservoir simulation, II: The multiblock approch in the integrated parallel accurate reservoir simulation (IPARS), the IMA Volumes in Mathematics and its Applications, Springer New York, 919, 2002.
 [27] Y.S Wu, N. Li, C. Wang, Q. Ran, J. Li, J. Yuan, et al., A multiplecontinuum model for simulation of gas production from shale gas reservoirs, SPE Reservoir Characterization and Simulation Conference and Exhibition, Society of Petroleum Engineers, 2013.
 [28] Y.S. Wu, G. Qin, R.E.Ewing, Y. Efendiev, Z. Kang, Y. Ren, et al., A multiplecontinuum approach for modeling multiphase flow in naturally fractured vuggy petroleum reservoirs: International Oil & Gas Conference and Exhibition in China, Society of Petroleum Engineers, 2006.
 [29] Y.S. Wu, Y. Di, Z. Kang, P. Fakcharoenphol, A multiplecontinuum model for simulating singlephase and multiphase flow in naturally fractured vuggy reservoirs, Journal of Petroleum Science and Engineering, 78(1), 2011, 1322.
 [30] Y.S. Wu, K. Pruess, et al., A multipleporosity method for simulation of naturally fractured petroleum reservoir, SPE Reservoir Engineering 3(01), 1998, 327336.
 [31] Y.S. Wu, J. Li, D. Ding, C. Wang, Y. Di, et al., A generalized framework model for the simulation of gas production in unconventional gas reservoirs, SPE Journal 19(05), 2014, 845857.
 [32] L. Zhang, A parallel algorithm for adaptive local refinement of tetrahedral meshes using bisection, Numer. Math.: Theory, Methods and applications, Vol. 2, 6589 2009.
 [33] L. Zhang, T. Cui and H. Hui, A set of symmetric quadrature rules on triangles and tetrahedra, J. Comput. Math. 2009, 27(1), 8996.
 [34] CMG Ltd, IMEX User Guid, Version 2014, 2014.
Comments
There are no comments yet.