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

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-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

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)

Domain Decomposition And Time-Splitting Methods For The Biot System Of Poroelasticity

In this thesis, we develop efficient mixed finite element methods to solve the Biot system of poroelasticity, which models the flow of a viscous fluid through a porous medium along with the deformation of the medium. We study non-overlapping domain decomposition techniques and sequential splitting methods to reduce the computational complexity of the problem. The solid deformation is modeled with a mixed three-field formulation with weak stress symmetry. The fluid flow is modeled with a mixed Darcy formulation. We introduce displacement and pressure Lagrange multipliers on the subdomain interfaces to impose weakly the continuity of normal stress and normal velocity, respectively. The global problem is reduced to an interface problem for the Lagrange multipliers, which is solved by a Krylov space iterative method. We study both monolithic and split methods. For the monolithic method, the cases of matching and non-matching subdomain grid interfaces are analyzed separately. For both cases, a coupled displacement-pressure interface problem is solved, with each iteration requiring the solution of local Biot problems. For the case of matching subdomain grids, we show that the resulting interface operator is positive definite and analyze the convergence of the iteration. For the non-matching subdomain grid case, we use a multiscale mortar mixed finite element (MMMFE) approach. We further study drained split and fixed stress Biot splittings, in which case we solve separate interface problems requiring elasticity and Darcy solves. We analyze the stability of the split formulations. We also use numerical experiments to illustrate the convergence of the domain decomposition methods and compare their accuracy and efficiency in the monolithic and time-splitting settings. Finally, we present a novel space-time domain decomposition technique for the mixed finite element formulation of a parabolic equation. This method is motivated by the MMMFE method, where we split the space-time domain into multiple subdomains with space-time grids of different sizes. Scalar Lagrange multiplier (mortar) functions are introduced to enforce weakly the continuity of the normal component of the mixed finite element flux variable over the space-time interfaces. We analyze the new method and numerical experiments are developed to illustrate and confirm the theoretical results

Anisotropic dual-continuum representations for multiscale poroelastic materials:Development and numerical modelling

Dual-continuum (DC) models can be tractable alternatives to explicit approaches for the numerical modelling of multiscale materials with multiphysics behaviours. This work concerns the conceptual and numerical modelling of poroelastically coupled dual-scale materials such as naturally fractured rock. Apart from a few exceptions, previous poroelastic DC models have assumed isotropy of the constituents and the dual-material. Additionally, it is common to assume that only one continuum has intrinsic stiffness properties. Finally, little has been done into validating whether the DC paradigm can capture the global poroelastic behaviours of explicit numerical representations at the DC modelling scale. We address the aforementioned knowledge gaps in two steps. First, we utilise a homogenisation approach based on Levin's theorem to develop a previously derived anisotropic poroelastic constitutive model. Our development incorporates anisotropic intrinsic stiffness properties of both continua. This addition is in analogy to anisotropic fractured rock masses with stiff fractures. Second, we perform numerical modelling to test the dual-continuum model against fine-scale explicit equivalents. In doing, we present our hybrid numerical framework, as well as the conditions required for interpretation of the numerical results. The tests themselves progress from materials with isotropic to anisotropic mechanical and flow properties. The fine-scale simulations show anisotropy can have noticeable effects on deformation and flow behaviour. However, our numerical experiments show the DC approach can capture the global poroelastic behaviours of both isotropic and anisotropic fine-scale representations