Resolving Wave Propagation in Anisotropic Poroelastic Media Using Graphical Processing Units (GPUs)

Biot's equations describe the physics of hydromechanically coupled systems establishing the widely recognized theory of poroelasticity. This theory has a broad range of applications in Earth and biological sciences as well as in engineering. The numerical solution of Biot's equations is challenging because wave propagation and fluid pressure diffusion processes occur simultaneously but feature very different characteristic time scales. Analogous to geophysical data acquisition, high resolution and three dimensional numerical experiments lately redefined state of the art. Tackling high spatial and temporal resolution requires a high-performance computing approach. We developed a multi- graphical processing units (GPU) numerical application to resolve the anisotropic elastodynamic Biot's equations that relies on a conservative numerical scheme to simulate, in a few seconds, wave fields for spatial domains involving more than 1.5 billion grid cells. We present a comprehensive dimensional analysis reducing the number of material parameters needed for the numerical experiments from ten to four. Furthermore, the dimensional analysis emphasizes the key material parameters governing the physics of wave propagation in poroelastic media. We perform a dispersion analysis as function of dimensionless parameters leading to simple and transparent dispersion relations. We then benchmark our numerical solution against an analytical plane wave solution. Finally, we present several numerical modeling experiments, including a three-dimensional simulation of fluid injection into a poroelastic medium. We provide the Matlab, symbolic Maple, and GPU CUDA C routines to reproduce the main presented results. The high efficiency of our numerical implementation makes it readily usable to investigate three-dimensional and high-resolution scenarios of practical applications.ISSN:2169-9313ISSN:0148-0227ISSN:2169-935

Numerical modelling of thermo-hydromechanical (THM) in deforming porous media for subsurface systems

The study of multiphase flow and heat flow in partially saturated porous media is important in environmental geomechanics engineering because of its relevance to consolidation of porous media in unsaturated zone. A numerical model which describes the thermo-hydro-mechanical (THM) coupled problems in deformable porous material with two-phase flow has been developed. The relationships between capillary pressure, saturation of water and relative permeabilities of water and gas, proposed by Brooks and Corey was used. An extended study of the numerical model, based on the COMES-GEO code was conducted recently to solve unsaturated problems in local condition of Kg. Puteh wellfield, Kota Bharu. This site is a potential shallow aquifer which contribute to the largest groundwater supply in Kota Bharu, Kelantan. Some numerical investigation on the proposed formulation is discussed with illustrative example problems to demonstrate solution procedures and validating of the model

Development of robust and efficient solution strategies for coupled problems

Det er mange modeller i moderne vitenskap hvor sammenkoblingen mellom forskjellige fysiske prosesser er svært viktig. Disse finner man for eksempel i forbindelse med lagring av karbondioksid i undervannsreservoarer, flyt i kroppsvev, kreftsvulstvekst og geotermisk energiutvinning. Denne avhandlingen har to fokusområder som er knyttet til sammenkoblede modeller. Det første er å utvikle pålitelige og effektive tilnærmingsmetoder, og det andre er utviklingen av en ny modell som tar for seg flyt i et porøst medium som består av to forskjellige materialer. For tilnærmingsmetodene har det vært et spesielt fokus på splittemetoder. Dette er metoder hvor hver av de sammenkoblede modellene håndteres separat, og så itererer man mellom dem. Dette gjøres i hovedsak fordi man kan utnytte tilgjengelig teori og programvare for å løse hver undermodell svært effektivt. Ulempen er at man kan ende opp med løsningsalgoritmer for den sammenkoblede modellen som er trege, eller ikke kommer frem til noen løsning i det hele tatt. I denne avhandlingen har tre forskjellige metoder for å forbedre splittemetoder blitt utviklet for tre forskjellige sammenkoblede modeller. Den første modellen beskriver flyt gjennom deformerbart porøst medium og er kjent som Biot ligningene. For å anvende en splittemetode på denne modellen har et stabiliseringsledd blitt tilført. Dette sikrer at metoden konvergerer (kommer frem til en løsning), men dersom man ikke skalerer stabiliseringsleddet riktig kan det ta veldig lang tid. Derfor har et intervall hvor den optimale skaleringen av stabiliseringsleddet befinner seg blitt identifisert, og utfra dette presenteres det en måte å praktisk velge den riktige skaleringen på. Den andre modellen er en fasefeltmodell for sprekkpropagering. Denne modellen løses vanligvis med en splittemetode som er veldig treg, men konvergent. For å forbedre dette har en ny akselerasjonsmetode har blitt utviklet. Denne anvendes som et postprosesseringssteg til den klassiske splittemetoden, og utnytter både overrelaksering og Anderson akselerasjon. Disse to forskjellige akselerasjonsmetodene har kompatible styrker i at overrelaksering akselererer når man er langt fra løsningen (som er tilfellet når sprekken propagerer), og Anderson akselerasjon fungerer bra når man er nærme løsningen. For å veksle mellom de to metodene har et kriterium basert på residualfeilen blitt brukt. Resultatet er en pålitelig akselerasjonsmetode som alltid akselererer og ofte er svært effektiv. Det siste modellen kalles Cahn-Larché ligningene og er også en fasefeltmodell, men denne beskriver elastisitet i et medium bestående av to elastiske materialer som kan bevege seg basert på overflatespenningen mellom dem. Dette problemet er spesielt utfordrende å løse da det verken er lineært eller konvekst. For å håndtere dette har en ny måte å behandle tidsavhengigheten til det underliggende koblede problemet på blitt utviklet. Dette leder til et diskret system som er ekvivalent med et konvekst minimeringsproblem, som derfor er velegnet til å løses med de fleste numeriske optimeringsmetoder, også splittemetoder. Den nye modellen som har blitt utviklet er en utvidelse av Cahn-Larché ligningene og har fått navnet Cahn-Hilliard-Biot. Dette er fordi ligningene utgjør en fasefelt modell som beskriver flyt i et deformerbart porøst medium med to poroelastiske materialer. Disse kan forflytte seg basert på overflatespenning, elastisk spenning, og poretrykk, og det er tenkt at modellen kan anvendes i forbindelse med kreftsvulstmodellering.There are many applications where the study of coupled physical processes is of great importance. These range from the life sciences with flow in deformable human tissue to structural engineering with fracture propagation in elastic solids. In this doctoral dissertation, there is a twofold focus on coupled problems. Firstly, robust and efficient solution strategies, with a focus on iterative decoupling methods, have been applied to several coupled systems of equations. Secondly, a new thermodynamically consistent coupled system of equations is proposed. Solution strategies are developed for three different coupled problems; the quasi-static linearized Biot equations that couples flow through porous materials and elastic deformation of the solid medium, variational phase-field models for brittle fracture that couple a phase-field equation for fracture evolution with linearized elasticity, and the Cahn-Larché equations that model elastic effects in a two-phase elastic material and couples an extended Cahn-Hilliard phase-field equation and linearized elasticity. Finally, the new system of equations that is proposed models flow through a two-phase deformable porous material where the solid phase evolution is governed by interfacial forces as well as effects from both the fluid and elastic properties of the material. In the work that concerns the quasi-static linearized Biot equations, the focus is on the fixed-stress splitting scheme, which is a popular method for sequentially solving the flow and elasticity subsystems of the full model. Using such a method is beneficial as it allows for the use of readily available solvers for the subproblems; however, a stabilizing term is required for the scheme to converge. It is well known that the convergence properties of the method strongly depend on how this term is chosen, and here, the optimal choice of it is addressed both theoretically and practically. An interval where the optimal stabilization parameter lies is provided, depending on the material parameters. In addition, two different ways of optimizing the parameter are proposed. The first is a brute-force method that relies on the mesh independence of the scheme's optimal stabilization parameter, and the second is valid for low-permeable media and utilizes an equivalence between the fixed-stress splitting scheme and the modified Richardson iteration. Regarding the variational phase-field model for brittle fracture propagation, the focus is on improving the convergence properties of the most commonly used solution strategy with an acceleration method. This solution strategy relies on a staggered scheme that alternates between solving the elasticity and phase-field subproblems in an iterative way. This is known to be a robust method compared to the monolithic Newton method. However, the staggered scheme often requires many iterations to converge to satisfactory precision. The contribution of this work is to accelerate the solver through a new acceleration method that combines Anderson acceleration and over-relaxation, dynamically switching back and forth between them depending on a criterion that takes the residual evolution into account. The acceleration scheme takes advantage of the strengths of both Anderson acceleration and over-relaxation, and the fact that they are complementary when applied to this problem, resulting in a significant speed-up of the convergence. Moreover, the method is applied as a post-processing technique to the increments of the solver, and can thus be implemented with minor modifications to readily available software. The final contribution toward solution strategies for coupled problems focuses on the Cahn-Larché equations. This is a model for linearized elasticity in a medium with two elastic phases that evolve with respect to interfacial forces and elastic effects. The system couples linearized elasticity and an extended Cahn-Hilliard phase-field equation. There are several challenging features with regards to solution strategies for this system including nonlinear coupling terms, and the fourth-order term that comes from the Cahn-Hilliard subsystem. Moreover, the system is nonlinear and non-convex with respect to both the phase-field and the displacement. In this work, a new semi-implicit time discretization that extends the standard convex-concave splitting method applied to the double-well potential from the Cahn-Hilliard subsystem is proposed. The extension includes special treatment for the elastic energy, and it is shown that the resulting discrete system is equivalent to a convex minimization problem. Furthermore, an alternating minimization solver is proposed for the fully discrete system, together with a convergence proof that includes convergence rates. Through numerical experiments, it becomes evident that the newly proposed discretization method leads to a system that is far better conditioned for linearization methods than standard time discretizations. Finally, a new model for flow through a two-phase deformable porous material is proposed. The two poroelastic phases have distinct material properties, and their interface evolves according to a generalized Ginzburg–Landau energy functional. As a result, a model that extends the Cahn-Larché equations to poroelasticity is proposed, and essential coupling terms for several applications are highlighted. These include solid tumor growth, biogrout, and wood growth. Moreover, the coupled set of equations is shown to be a generalized gradient flow. This implies that the system is thermodynamically consistent and makes a toolbox of analysis and solvers available for further study of the model.Doktorgradsavhandlin

Finite volume discretization for poroelastic media with fractures modeled by contact mechanics

A fractured poroelastic body is considered where the opening of the fractures is governed by a nonpenetration law, whereas slip is described by a Coulomb‐type friction law. This physical model results in a nonlinear variational inequality problem. The variational inequality is rewritten as a complementary function, and a semismooth Newton method is used to solve the system of equations. For the discretization, we use a hybrid scheme where the displacements are given in terms of degrees of freedom per element, and an additional Lagrange multiplier representing the traction is added on the fracture faces. The novelty of our method comes from combining the Lagrange multiplier from the hybrid scheme with a finite volume discretization of the poroelastic Biot equation, which allows us to directly impose the inequality constraints on each subface. The convergence of the method is studied for several challenging geometries in 2D and 3D, showing that the convergence rates of the finite volume scheme do not deteriorate when it is coupled to the Lagrange multipliers. Our method is especially attractive for the poroelastic problem because it allows for a straightforward coupling between the matrix deformation, contact conditions, and fluid pressure.publishedVersio

A Study of Interwell Interference and Well Performance in Unconventional Reservoirs Based on Coupled Flow and Geomechanics Modeling with Improved Computational Efficiency

Completion quality of tightly spaced horizontal wells in unconventional reservoirs is important for hydrocarbon recovery efficiency. Parent well production usually leads to heterogeneous stress evolution around parent wells and at infill well locations, which affects hydraulic fracture growth along infill wells. Recent field observations indicate that infill well completions lead to frac hits and production interference between parent and infill wells. Therefore, it is important to characterize the heterogeneous interwell stress/pressure evolutions and hydraulic fracture networks. This work presents a reservoir-geomechanics-fracturing modeling workflow and its implementation in unconventional reservoirs for the characterization of interwell stress and pressure evolutions and for the modeling of interwell hydraulic fracture geometry. An in-house finite element model coupling fluid flow and geomechanics is first introduced and used to characterize production-induced stress and pressure changes in the reservoir. Then, an in-house complex fracture propagation model coupling fracture mechanics and wellbore/fracture fluid flow is used for the simulation of hydraulic fractures along infill wells. A parallel solver is also implemented in a reservoir geomechanics simulator in a separate study to investigate the potential of improving computational efficiency. Results show that differential stress (DS), parent well fracture geometry, legacy production time, bottomhole pressure (BHP) for legacy production, and perforation cluster location are key parameters affecting interwell fracture geometry and the occurrence of frac hits. In general, transverse infill well fractures are obtained in scenarios with large DS and small legacy producing time/BHP. Non-uniform parent well fracture geometry leads to frac hits in certain cases, while the assumption of uniform parent well fracture half-lengths in the numerical model could not capture the phenomenon of frac hits. Perforation cluster locations along infill wells do not play an important role in determining whether an infill well hydraulic fracture is transverse, while they are important for the occurrence of frac hits. In addition, the implementation of a parallel solver, PETSc, in a fortran-based simulator indicates that an overall speedup of 14 can be achieved for simulations with one million grid blocks. This result provides a reference for improving computational efficiency for geomechanical simulation involving large matrices using finite element methods (FEM)

A-posteriori error estimation in axisymmetric geotechnical analyses

In this paper, an a-posteriori error estimator suitable for use in axisymmetric geotechnical analyses has been developed. The consolidation superconvergent patch recovery with equilibrium and boundaries (CSPREB) method, developed for plane-strain coupledconsolidation problems, is extended to axisymmetric analyses. The use of pore pressures in the error estimator was found to improve results when predicting consolidation. Collapse loads under undrained soil conditions are known to be over-predicted due to “locking”, especially in axial symmetry where there are further displacement constraints. The proposed solution technique reduced locking slightly, but could not eliminate it, as it is inherent in the displacement formulation for lower order elements

Geomechanics-Reservoir Modeling by Displacement Discontinuity-Finite Element Method

There are two big challenges which restrict the extensive application of fully coupled geomechanics-reservoir modeling. The first challenge is computational effort. Consider a 3-D simulation combining pressure and heat diffusion, elastoplastic mechanical response, and saturation changes; each node has at least 5 degrees of freedom, each leading to a separate equation. Furthermore, regions of large p, T and σ′ gradients require small-scale discretization for accurate solutions, greatly increasing the number of equations. When the rock mass surrounding the reservoir region is included, it is represented by many elements or nodes. These factors mean that accurate analysis of realistic 3-D problems is challenging, and will so remain as we seek to solve larger and larger coupled problems involving nonlinear responses. To overcome the first challenge, the displacement discontinuity method is introduced wherein a large-scale 3-D case is divided into a reservoir region where Δp, ΔT and non-linear effects are critical and analyzed using FEM, and an outside region in which the reservoir is encased where Δp and ΔT effects are inconsequential and the rock may be treated as elastic, analyzed with a 3D displacement discontinuity formulation. This scheme leads to a tremendous reduction in the degrees of freedom, yet allows for reasonably rigorous incorporation of the reactions of the surrounding rock. The second challenge arises from some forms of numerical instability. There are actually two types of sharp gradients implied in the transient advection-diffusion problem: one is caused by the high Peclet numbers, the other by the sharp gradient which appears during the small time steps due to the transient solution. The way to eliminate the spurious oscillations is different when the sharp gradients are induced by the transient evolution than when they are produced by the advective terms, and existing literature focuses mainly on eliminating the spurious spatial temperature oscillations caused by advection-dominated flow. To overcome the second challenge, numerical instability sources are addressed by introducing a new stabilized finite element method, the subgrid scale/gradient subgrid scale (SGS/GSGS) method

A Study of Interwell Interference and Well Performance in Unconventional Reservoirs Based on Coupled Flow and Geomechanics Modeling with Improved Computational Efficiency

Completion quality of tightly spaced horizontal wells in unconventional reservoirs is important for hydrocarbon recovery efficiency. Parent well production usually leads to heterogeneous stress evolution around parent wells and at infill well locations, which affects hydraulic fracture growth along infill wells. Recent field observations indicate that infill well completions lead to frac hits and production interference between parent and infill wells. Therefore, it is important to characterize the heterogeneous interwell stress/pressure evolutions and hydraulic fracture networks. This work presents a reservoir-geomechanics-fracturing modeling workflow and its implementation in unconventional reservoirs for the characterization of interwell stress and pressure evolutions and for the modeling of interwell hydraulic fracture geometry. An in-house finite element model coupling fluid flow and geomechanics is first introduced and used to characterize production-induced stress and pressure changes in the reservoir. Then, an in-house complex fracture propagation model coupling fracture mechanics and wellbore/fracture fluid flow is used for the simulation of hydraulic fractures along infill wells. A parallel solver is also implemented in a reservoir geomechanics simulator in a separate study to investigate the potential of improving computational efficiency. Results show that differential stress (DS), parent well fracture geometry, legacy production time, bottomhole pressure (BHP) for legacy production, and perforation cluster location are key parameters affecting interwell fracture geometry and the occurrence of frac hits. In general, transverse infill well fractures are obtained in scenarios with large DS and small legacy producing time/BHP. Non-uniform parent well fracture geometry leads to frac hits in certain cases, while the assumption of uniform parent well fracture half-lengths in the numerical model could not capture the phenomenon of frac hits. Perforation cluster locations along infill wells do not play an important role in determining whether an infill well hydraulic fracture is transverse, while they are important for the occurrence of frac hits. In addition, the implementation of a parallel solver, PETSc, in a fortran-based simulator indicates that an overall speedup of 14 can be achieved for simulations with one million grid blocks. This result provides a reference for improving computational efficiency for geomechanical simulation involving large matrices using finite element methods (FEM)

Rotation-Based Mixed Formulations for an Elasticity-Poroelasticity Interface Problem

In this paper we introduce a new formulation for the stationary poroelasticity equations written using the rotation vector and the total fluid-solid pressure as additional unknowns, and we also write an extension to the elasticity-poroelasticity problem. The transmission conditions are imposed naturally in the weak formulation, and the analysis of the effective governing equations is conducted by an application of Fredholm's alternative. We also propose a monolithically coupled mixed finite element method for the numerical solution of the problem. Its convergence properties are rigorously derived and subsequently confirmed by a set of computational tests that include applications to subsurface flow in reservoirs as well as to dentistry-oriented problems.Fondo Nacional de Desarrollo Científico y Tecnológico/[11160706]/FONDECYT/ChilePrograma Concurso Apoyo a Centros Científicos y Tecnológicos de Excelencia/[AFB170001]/PIA/ChileUCR::Sedes Regionales::Sede de OccidenteUCR::Vicerrectoría de Docencia::Ciencias Básicas::Facultad de Ciencias::Escuela de Matemátic