Anisotropy of hydraulic properties of the subsurface geologic formations is an essential feature that has been established as a consequence of the different geologic processes that undergo during the longer geologic time scale. With respect to subsurface reservoirs, in many cases, anisotropy plays significant role in dictating the direction of flow that becomes no longer dependent only on driving forces like the pressure gradient and gravity but also on the principal directions of anisotropy. Therefore, there has been a great deal of motivation to consider anisotropy into the subsurface flow and transport models. In this dissertation, we present subsurface flow modeling in single and dual continuum anisotropic porous media, which include the single-phase groundwater flow coupled with the solute transport in anisotropic porous media, the two-phase flow with gravity effect in anisotropic porous media, and the natural gas flow in anisotropic shale reservoirs. We have employed the multipoint flux approximation (MPFA) method to handle anisotropy in the flow model. The MPFA method is designed to provide correct discretization of the flow equations for general orientation of the principal directions of the permeability tensor. The implementation of MPFA method is combined with the experimenting pressure field approach, a newly developed technique that enables the solution of the global problem breaks down into the solution of multitude of local problems. The numerical results of the study demonstrate the significant effects of anisotropy of the subsurface formations. For the single-phase groundwater flow coupled with the solute transport modeling in anisotropic porous media, the results shows the strong impact of anisotropy on the pressure field and the migration of the solute concentration. For the two-phase flow modeling with gravity effect in anisotropic porous media, it is observed that the buoyancy-driven flow, which emerges due to the density differences between the phases, migrates upwards and the anisotropy aligns the flow directions closer to the principal direction of anisotropy. Lastly, for the gas flow modeling in anisotropic shale reservoirs, we observe that anisotropy affects the pressure fields and the velocity fields of the matrix and fracture systems as well as the production rate and cumulative production. It is observed from the results that all of the anisotropic cases produce higher amount of gas compared to isotropic case during the same production time. Furthermore, we have also examined the performance of MPFA with respect to mixed finite element (MFE) method over the lowest-order Raviart-Thomas (RT0) space and the first-order Brezzi-Douglas-Marini (BDM1) space. From the comparison of the numerical results we observe that MPFA method show very good agreement with the BDM1 than RT0. In terms of numerical implementation, however, MPFA method is easier than BDM1 and it also offers explicit discrete fluxes that are advantageous. Combining MPFA with the experimenting pressure field approach will certainly adds another advantage of implementing MPFA method as compared with RT0 and BDM1. Moreover, the computational cost (CPU cost) of the three different methods are also discussed.

In this research we consider developing a reservoir simulator capable of simulating complex coupled poromechanical processes on massively parallel computers. A variety of problems arising from petroleum and environmental engineering inherently necessitate the understanding of interactions between fluid flow and solid mechanics. Examples in petroleum engineering include reservoir compaction, wellbore collapse, sand production, and hydraulic fracturing. In environmental engineering, surface subsidence, carbon sequestration, and waste disposal are also coupled poromechanical processes. These economically and environmentally important problems motivate the active pursuit of robust, efficient, and accurate simulation tools for coupled poromechanical problems. Three coupling approaches are currently employed in the reservoir simulation community to solve the poromechanics system, namely, the fully implicit coupling (FIM), the explicit coupling, and the iterative coupling. The choice of the coupling scheme significantly affects the efficiency of the simulator and the accuracy of the solution. We adopt the fixed-stress iterative coupling scheme to solve the coupled system due to its advantages over the other two. Unlike the explicit coupling, the fixed-stress split has been theoretically proven to converge to the FIM for linear poroelasticity model. In addition, it is more efficient and easier to implement than the FIM. Our computational results indicate that this approach is also valid for multiphase flow. We discretize the quasi-static linear elasticity model for geomechanics in space using the continuous Galerkin (CG) finite element method (FEM) on general hexahedral grids. Fluid flow models are discretized by locally mass conservative schemes, specifically, the mixed finite element method (MFE) for the equation of state compositional flow on Cartesian grids and the multipoint flux mixed finite element method (MFMFE) for the single phase and two-phase flows on general hexahedral grids. While both the MFE and the MFMFE generate cell-centered stencils for pressure, the MFMFE has advantages in handling full tensor permeabilities and general geometry and boundary conditions. The MFMFE also obtains accurate fluxes at cell interfaces. These characteristics enable the simulation of more practical problems. For many reservoir simulation applications, for instance, the carbon sequestration simulation, we need to account for thermal effects on the compositional flow phase behavior and the solid structure stress evolution. We explicitly couple the poromechanics equations to a simplified energy conservation equation. A time-split scheme is used to solve heat convection and conduction successively. For the convection equation, a higher order Godunov method is employed to capture the sharp temperature front; for the conduction equation, the MFE is utilized. Simulations of coupled poromechanical or thermoporomechanical processes in field scales with high resolution usually require parallel computing capabilities. The flow models, the geomechanics model, and the thermodynamics model are modularized in the Integrated Parallel Accurate Reservoir Simulator (IPARS) which has been developed at the Center for Subsurface Modeling at the University of Texas at Austin. The IPARS framework handles structured (logically rectangular) grids and was originally designed for element-based data communication, such as the pressure data in the flow models. To parallelize the node-based geomechanics model, we enhance the capabilities of the IPARS framework for node-based data communication. Because the geomechanics linear system is more costly to solve than those of flow and thermodynamics models, the performance of linear solvers for the geomechanics model largely dictates the speed and scalability of the coupled simulator. We use the generalized minimal residual (GMRES) solver with the BoomerAMG preconditioner from the hypre library and the geometric multigrid (GMG) solver from the UG4 software toolbox to solve the geomechanics linear system. Additionally, the multilevel k-way mesh partitioning algorithm from METIS is used to generate high quality mesh partitioning to improve solver performance. Numerical examples of coupled poromechanics and thermoporomechanics simulations are presented to show the capabilities of the coupled simulator in solving practical problems accurately and efficiently. These examples include a real carbon sequestration field case with stress-dependent permeability, a synthetic thermoporoelastic reservoir simulation, poroelasticity simulations on highly distorted hexahedral grids, and parallel scalability tests on a massively parallel computer. / text

Submitted by Fabio Sobreira Campos da Costa (fabio.sobreira@ufpe.br) on 2016-07-01T15:05:14Z No. of bitstreams: 2 license_rdf: 1232 bytes, checksum: 66e71c371cc565284e70f40736c94386 (MD5) Souza_Tese_2015_09_22.pdf: 8187999 bytes, checksum: 664629aed28d692dce410fefbfe793dc (MD5) / Made available in DSpace on 2016-07-01T15:05:14Z (GMT). No. of bitstreams: 2 license_rdf: 1232 bytes, checksum: 66e71c371cc565284e70f40736c94386 (MD5) Souza_Tese_2015_09_22.pdf: 8187999 bytes, checksum: 664629aed28d692dce410fefbfe793dc (MD5) Previous issue date: 2015-09-22 / Anp / Sob certas hipóteses simplificadoras, o modelo matemático que descreve o escoamento de água e óleo em reservatórios de petróleo pode ser representado por um sistema não linear de Equações Diferenciais Parciais composto por uma equação elíptica de pressão (fluxo) e uma equação hiperbólica de saturação (transporte). Devido a complexidades na modelagem de ambientes deposicionais, nos quais são incluídos camadas inclinadas, canais, falhas e poços inclinados, há uma dificuldade de se construir um modelo que represente adequadamente certas características dos reservatórios, especialmente quando malhas estruturadas são usadas (cartesianas ou corner point). Além disso, a modelagem do escoamento multifásico nessas estruturas geológicas incluem descontinuidades na variável e instabilidades no escoamento, associadas à elevadas razões de mobilidade e efeitos de orientação de malha. Isso representa um grande desafio do ponto de vista numérico. No presente trabalho, uma formulação fundamentada no Método de Volumes Finitos é estudada e proposta para discretizar as equações elíptica de pressão e hiperbólica de saturação. Para resolver a equação de pressão três formulações robustas, com aproximação dos fluxos por múltiplos pontos são estudadas. Essas formulações são abeis para lidar com tensores de permeabilidade completos e malhas poligonais arbitrárias, sendo portanto uma generalização de métodos mais tradicionais com aproximação do fluxo por apenas dois pontos. A discretização da equação de saturação é feita com duas abordagens com característica multidimensional. Em uma abordagem mais convencional, os fluxos numéricos são extrapolados diretamente nas superfícies de controle por uma aproximação de alta resolução no espaço (2ª a 4ª ordem) usando uma estratégia do tipo MUSCL. Uma estratégia baseada na Técnica de Mínimos Quadrados é usada para a reconstrução polinomial. Em uma segunda abordagem, uma variação de uma esquema numérico Verdadeiramente Multidimensional é proposto. Esse esquema diminui o efeito de orientação de malha, especialmente para malhas ortogonais, mesmo embora alguma falta de robustez possa ser observada pra malhas excessivamente distorcidas. Nesse tipo de formulação, os fluxos numéricos são calculados de uma forma multidimensional. Consiste em uma combinação convexa de valores de saturação ou fluxo fracionário, seguindo a orientação do escoamento através do domínio computacional. No entanto, a maioria dos esquemas numéricos achados na literatura tem aproximação apenas de primeira ordem no espaço e requer uma solução implícita de sistemas algébricos locais. Adicionalmente, no presente texto, uma forma modificada desses esquemas “Verdadeiramente” Multidimensionais é proposta em um contexto centrado na célula. Nesse caso, os fluxos numéricos multidimensionais são calculados explicitamente usando aproximações de alta ordem no espaço. Para o esquema proposto, a robustez e o caráter multidimensional também leva em conta a distorção da malha por meio de uma ponderação adaptativa. Essa ponderação regula a característica multidimensional da formulação de acordo com a distorção da malha. Claramente, os efeitos de orientação de malha são reduzidos. A supressão de oscilações espúrias, típicas de aproximações de alta ordem, são obtidas usando, pela primeira vez no contexto de simulação de reservatórios, uma estratégia de limitação multidimensional ou Multidimensional Limiting Process (MLP). Essa estratégia garante soluções monótonas e podem ser usadas em qualquer malha poligonal, sendo naturalmente aplicada em aproximações de ordem arbitrária. Por fim, de modo a garantir soluções convergentes, mesmo para problemas tipicamente não convexos, associados ao modelo de Buckley-Leverett, uma estratégia robusta de correção de entropia é empregada. O desempenho dessas formulações é verificado com a solução de problemas relevantes achados na literatura. / Under certain simplifying assumptions, the problem that describes the fluid flow of oil and water in heterogeneous and anisotropic petroleum reservoir can be described by a system of non-linear partial differential equations that comprises an elliptic pressure equation (flow) and a hyperbolic saturation equation (transport). Due to the modeling of complex depositional environments, including inclined laminated layers, channels, fractures, faults and the geometrical modeling of deviated wells, it is difficult to properly build and handle the Reservoir Characterization Process (RCM), particularly by using structured meshes (cartesian or corner point), which is the current standard in petroleum reservoir simulators. Besides, the multiphase flow in such geological structures includes the proper modeling of water saturation shocks and flow instabilities associated to high mobility ratios and Grid Orientation Effects (GOE), posing a great challenge from a numerical point of view. In this work, a Full Finite Volume Formulation is studied and proposed to discretize both, the elliptic pressure and the hyperbolic saturation equations. To solve the pressure equation, we study and use three robust Multipoint Flux Approximation Methods (MPFA) that are able to deal with full permeability tensors and arbitrary polygonal meshes, making it relatively easy to handle complex geological structures, inclined wells and mesh adaptivity in a natural way. To discretize the saturation equation, two different multidimensional approaches are employed. In a more conventional approach, the numerical fluxes are extrapolated directly on the control surfaces for a higher resolution approximation in space (2nd to 4th order) by a MUSCL (Monotone Upstream Centered Scheme for Conservation Laws) procedure. A least squares based strategy is employed for the polynomial reconstruction. In a second approach, a variation of a “Truly” Multidimensional Finite Volume method is proposed. This scheme diminishes GOE, especially for orthogonal grids, even though some lack of robustness can be observed for extremely distorted meshes. In this type of scheme, the numerical flux is computed in each control surface in a multidimensional way, by a convex combination of the saturation or the fractional flow values, following the approximate wave orientation throughout the computational domain. However, the majority of the schemes found in literature is only first order accurate in space and demand the implicit solution of local conservation problems. In the present text, a Modified Truly Multidimensional Finite Volume Method (MTM-FVM) is proposed in a cell centered context. The truly multidimensional numerical fluxes are explicitly computed using higher order accuracy in space. For the proposed scheme, the robustness and the multidimensional character of the aforementioned MTM-FVM explicitly takes into account the angular distortion of the computational mesh by means of an adaptive weight, that tunes the multidimensional character of the formulation according to the grid distortion, clearly diminishing GOE. The suppression of the spurious oscillations, typical from higher order schemes, is achieved by using for the first time in the context of reservoir simulation a Multidimensional Limiting Process (MLP). The MLP strategy formally guarantees monotone solutions and can be used with any polygonal mesh and arbitrary orders of approximation. Finally, in order to guarantee physically meaningful solutions, a robust “entropy fix” strategy is employed. This produces convergent solutions even for the typical non-convex flux functions that are associated to the Buckley-Leverett problem. The performance of the proposed full finite volume formulation is verified by solving some relevant benchmark problems.