Numerical Solution of Compositional Two-Phase Flow in Porous Media - - PowerPoint PPT Presentation
Numerical Solution of Compositional Two-Phase Flow in Porous Media P. Bastian Collaborators: M. Blatt, O. Ippisch, R. Neumann Universit at Heidelberg Interdisziplin ares Zentrum f ur Wissenschaftliches Rechnen Im Neuenheimer Feld 368,
Collaborators: M. Blatt, O. Ippisch, R. Neumann
Universit¨ at Heidelberg Interdisziplin¨ ares Zentrum f¨ ur Wissenschaftliches Rechnen Im Neuenheimer Feld 368, D-69120 Heidelberg email: Peter.Bastian@iwr.uni-heidelberg.de
Linz, October 6, 2011
Compositional Two-Phase Flow Linz, October 6, 2011 1 / 34
Overview
1 Discontinuous Galerkin method for incompressible pure
two-phase flow
◮ Implicit, fully-coupled ◮ Media discontinuities, internal Signorini-type interface
conditions
2 Compositional two-phase flow with equilibrium phase exchange ◮ Phase disappearance problem, obstacle-type problem ◮ No switching, no additional variables 3 Software framework DUNE ◮ Rapid prototyping ◮ Many different meshes, adaptive refinement,
dimension-independence
◮ Parallelization and efficient solvers
Compositional Two-Phase Flow Linz, October 6, 2011 2 / 34
DG for Two-Phase Flow
Setting
◮ Both phases incompressible ◮ No dissolution, pure two-phase flow ◮ 1d, 2d, 3d, including gravity ◮ Different entry pressures in subdomains (media heterogeneities)
Why DG ?
◮ Potentially higher order of convergence ◮ Full tensors ◮ Local conservation ◮ Handles elliptic, parabolic and hyperbolic equations ◮ nonconforming, unstructured grids / refinement
Compositional Two-Phase Flow Linz, October 6, 2011 3 / 34
DG for Two-Phase Flow
Bastian & Rivi` ere ’2004
◮ Global pressure / saturation formulation, splitting ◮ Implicit/explicit saturation(+limiters), H(div)-projection
Eslinger ’2005
◮ Splitting: Implicit Pressure, Implicit/explicit saturation ◮ Kirchhoff transformation, media disc., gravity, compressible
Epshteyn & Rivi` ere ’2007
◮ Fully implicit, Fully-coupled approach, pw, sn formulations ◮ no media discontinuities, no gravity, very coarse grids
Ern, Mozolevski, Schuh ’2010
◮ Splitting, global pressure, implicit saturation, H(div)-projection ◮ Media discontinuities, 1D, no gravity
This work
◮ Fully- coupled, higher-order in time, pw, pc formulation, ◮ Media discontinuities, 1-3d
Compositional Two-Phase Flow Linz, October 6, 2011 4 / 34
DG for Two-Phase Flow
P = {w, n}. Two coupled equations for pw, pc: −∇ · {(λw + λn)K∇pw + λnK∇pc − (λwρw + λnρn)Kg} = f ∂t{φ(1 − ψ(pc))} − ∇ · {λnK(∇pw + ∇pc − ρng)} = fn Nonlinearities: sw = ψ(pc), λw(pc) = krw(ψ(pc)) µw , λn(pc) = krn(1 − ψ(pc)) µn . pc sw 1 pe ψ(pc) sw 1 1 krw krn pc < pe: sn = 0, ∂t(φsn) = 0, λn = 0 → singular! Allow ψ(pc) > 1.
Compositional Two-Phase Flow Linz, October 6, 2011 5 / 34
DG for Two-Phase Flow
Two regions with different ψ(pc) curves: Ω ψ(l) ψ(h) nF pc sw 1 p(l)
e
ψ(l) p(h)
e
ψ(h) s∗
w
Interface conditions (van Duijn et al., ’1995): p(h)
c
=
c
p(l)
c
≥ p(h)
e
p(h)
e
else At the media interface let nF point into the (l) region and set J(pc) =
c
− p(l)
c
p(l)
c
≥ p(h)
e
p(h)
c
− p(h)
e
else
Compositional Two-Phase Flow Linz, October 6, 2011 6 / 34
DG for Two-Phase Flow
Setting u = −(λw + λn)K∇pw − λnK∇pc + (λwρw + λnρn)Kg, then
−
u · ∇v +
h
nF · {u}ωv + (γFv − θnF · {λK∇v}ω)pw =
fv −
h
h
(γFv − θnF · (λK∇v))) pD
w
Weighted averages (Di Pietro, Ern, Guermond ’2008): δ±
Kn = nt FK ±nF,
ω− = δ+
Kn
δ−
Kn + δ+ Kn
, ω+ = δ−
Kn
δ−
Kn + δ+ Kn
. Penalty term (· denotes harmonic average): γF = λδKnM, M = αk(k + n − 1) |F| min(|T −(F)|, |T +(F)|)
Compositional Two-Phase Flow Linz, October 6, 2011 7 / 34
DG for Two-Phase Flow
Setting un = λnqn, qn = −K(∇pw + ∇pc − ρng), then
∂t
φ(1 − ψ(pc))w −
un · ∇w +
h
λ↑
nnF · {qn}ωw + (γn,Fw − θλ↑ nnF · {K∇w}ω)J(pc)
=
fnw −
h
h
(γn,Fw − θλnnF · (K∇w))) pD
c
Upwinding of mobility:
p↑
c =
c
nF · {qn}ω ≥ 0 p+
c
else , s±
n = 1 − ψ±(p↑ c ),
λ↑
n =
2k−
rn (s− n )k+ rn(s+ n )
(k−
rn (s− n ) + k+ rn(s+ n ))µn
.
Penalty parameter:
γn,F = λ↑
nδKnM,
F ∈ Γ, γn,F = λ−
n + λ+ n
2 δKnM, else.
Compositional Two-Phase Flow Linz, October 6, 2011 8 / 34
DG for Two-Phase Flow
Compositional Two-Phase Flow Linz, October 6, 2011 9 / 34
DG for Two-Phase Flow
Parameters
◮ Brooks Corey capillary pressure ◮ Quadratic relative permeability ◮ scalar absolute permeability ◮ P1 on quadrilaterals, 2nd order Alexander scheme in time ◮ P2 on quadrilaterals, 3rd order Alexander scheme in time
Reference scheme
◮ Cell-centered finite volume scheme, two-point flux
approximation, implicit Euler
◮ pw, pc formulation ◮ Upwinding of mobility, as described above
Meshes and timesteps
◮ From 40 × 24, ∆t = 60s to 1280x768, ∆t = 1.875,
T = 1800s, 3600s
◮ 2 million DOF, 2000 time steps for CCFV reference solution
Compositional Two-Phase Flow Linz, October 6, 2011 10 / 34
DG for Two-Phase Flow
Compositional Two-Phase Flow Linz, October 6, 2011 11 / 34
DG for Two-Phase Flow
Compositional Two-Phase Flow Linz, October 6, 2011 12 / 34
DG for Two-Phase Flow
Compositional Two-Phase Flow Linz, October 6, 2011 13 / 34
DG for Two-Phase Flow
Compositional Two-Phase Flow Linz, October 6, 2011 14 / 34
DG for Two-Phase Flow
Compositional Two-Phase Flow Linz, October 6, 2011 15 / 34
DG for Two-Phase Flow
Number of time steps (Newton convergence) in 2d: T = 1800s T = 3600s h−1 ∆t CCFV DG CCFV DG 40 60.000 30 30 60 60 80 30.000 60 60 120 120 160 15.000 120 120 240 240 320 7.500 240 260 557 566 640 3.750 480 960 1280 1.875 960 1925 Conclusion: Fully coupled DG scheme runs robustly and agrees with CCFV scheme Saves about two levels of refinement in space and time
Compositional Two-Phase Flow Linz, October 6, 2011 16 / 34
Compositional Two-Phase Flow
Consider two phases P = {l, g}, two components C = {w, g} and equilibrium phase exchange. Mole conservation for each component: φ∂t(νlxw
l sl + νgxw g sg) + ∇ · (νlxw l ul + νgxw g ug + jw l + jw g ) = f w
φ∂t(νlxg
l sl + νgxg g sg) + ∇ · (νlxg l ul + νgxg g ug + jg l + jg g ) = f g
Extended Darcy’s law and diffusion terms: uα = −λαK(∇pα − ραg), jκ
α = −(τφsαDκ α)να∇xκ α,
jw
α + jg α = 0.
Capillary pressure: pc = pg − pl, sl = ψ(pc), sl + sg = 1. Solubility (given in parametrized form): xg
l = ξ(pg, T, . . .),
xw
l = 1 − xg l ,
xw
g = η(pg, T, . . .),
xg
g = 1 − xw g .
Compositional Two-Phase Flow Linz, October 6, 2011 17 / 34
Compositional Two-Phase Flow
If both phases are present, 0 < sl, sg < 1, two variables describe the state of the system, e.g (pβ, sn) or (pβ, pc). If sg = 0 (gas phase disappears) then the balance equations reduce to φ∂tνl − ∇ · (νlµ−1
l K(∇pl − ρlg)) = f w + f g
φ∂t(νlxg
l ) + ∇ · (νlxg l ul + jg l ) = f g
→ Single phase flow + transport problem. Variables (pl, xg
l ) describe
the state of the system. Note that xg
l is now an independent variable.
Several ways to overcome this: Variable switching (Class, ’2000) Complementarity conditions (Jaffr´ e, Sboui, ’2010) Extended formulations (Abadbour & Panfilov, ’2009, Bourgeat et al., preprint)
Compositional Two-Phase Flow Linz, October 6, 2011 18 / 34
Compositional Two-Phase Flow
Introduced by Ippisch (2003) in his Ph.D. thesis:
1 In the two-phase region, pg, pc is a valid set of primary variables. 2 In the one-phase region (no gas phase), pl, xg
l are natural
primary variables. Consider a transformation of variables: pl = pg − pc, xg
l = ξ(pg).
With invertible ξ ∈ C 1, (pg, pc) can be used as alternative set of variables.
3 At the interface between one-phase and two-phase region pl, xg
l
are continuous and therefore also pg, pc. This choice is not unique. pg is preferred over pl because then ξ(pg) and ψ(pc) depend on just one variable.
Compositional Two-Phase Flow Linz, October 6, 2011 19 / 34
Compositional Two-Phase Flow
CCFV, two-point flux pc upwinding, harmonic averaging for corresponding mobility Fully implicit, fully coupled Newton method, numerical Jacobian Agglomeration based algebraic multigrid applied to full system Dimension-independent implementation Parallel Work by Rebecca Neumann
Compositional Two-Phase Flow Linz, October 6, 2011 20 / 34
Compositional Two-Phase Flow
P = {l, g}, C = {water, hydrogen}, no water in gas phase.
Γout noflux noflux Γin
sg and pressures at the inflow boundary over time:
0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0e+00 2e+05 4e+05 6e+05 8e+05 1e+06 Sn Time [years] Sn 700000 800000 900000 1e+06 1.1e+06 1.2e+06 1.3e+06 1.4e+06 1.5e+06 0e+00 2e+05 4e+05 6e+05 8e+05 1e+06 p[Pa] Time [years] pl pn
Compositional Two-Phase Flow Linz, October 6, 2011 21 / 34
Compositional Two-Phase Flow
2d-setup
Γout noflux noflux Γin
Van Genuchten ψ, krα Solubility after Spycher, Pruess, Ennis-King (2005) Liquid phase density after Garcia (2001) CO2 phase density after Duan et al. (1992) Liquid phase viscosity after Atkins (1990) CO2 phase viscosity after Fenghour et al. (1998) Tortuosity after Millington, Quirk (1961)
Compositional Two-Phase Flow Linz, October 6, 2011 22 / 34
Compositional Two-Phase Flow
7 days, max(sg) = 0.82, max(xg
l ) = 0.022
20 days, max(sg) = 1, max(xg
l ) = 0.02
80 days, max(sg) = 1, max(xg
l ) = 0.02
Compositional Two-Phase Flow Linz, October 6, 2011 23 / 34
Compositional Two-Phase Flow
Compositional Two-Phase Flow Linz, October 6, 2011 24 / 34
DUNE Software Framework
Distributed and Unified Numerics Environment Software for the numerical solution of PDEs with grid based methods. Goals: Flexibility: Meshes, discretizations, adaptivity, solvers. Efficiency: Pay only for functionality you need. Parallelization. Reuse of existing code. Enable team work through standardized interfaces.
Compositional Two-Phase Flow Linz, October 6, 2011 25 / 34
DUNE Software Framework
Developed since 2002 by groups at
◮ Free University of Berlin: O. Sander and R. Kornhuber ◮ Freiburg University: R. Kl¨
◮ Warwick: A. Dedner ◮ M¨
unster University: C. Engwer, M. Ohlberger
◮ Heidelberg University: M. Blatt, P. Bastian
Available under GNU LGPL license with linking exception. Platform for “Open Reservoir Simulator” (U Stuttgart, U Bergen, SINTEF, StatOil, . . . ) DUNE courses given every spring (at least).
Compositional Two-Phase Flow Linz, October 6, 2011 26 / 34
DUNE Software Framework
Separation of data structures and algorithms
Mesh Interface E.g. FE discretization Algorithm Structured grid Unstructured simplicial grid Unstructured multi−element grid
Realization with generic programming (templates) in C++. Static polymorphism:
◮ Inlining of “small” methods. ◮ Allows global optimizations. ◮ Interface code is removed at compile-time.
Template Meta Programs: compile-time algorithms. Standard Template Library (STL) is a prominent example.
Compositional Two-Phase Flow Linz, October 6, 2011 27 / 34
DUNE Software Framework
Major DUNE modules are:
ALU UG dune−grid−howto dune−fem dune−common Alberta NeuronGrid dune−pdelab−howto dune−pdelab dune−localfunctions VTK Gmsh SuperLU Metis dune−grid dune−istl
Compositional Two-Phase Flow Linz, October 6, 2011 28 / 34
DUNE Software Framework
Alberta 2d, 3d, ALU3dGrid, simplices, cubes UG 2d, simplices, cubes, 3d, simplices, cubes (Viz.: ParaView/VTK)
Compositional Two-Phase Flow Linz, October 6, 2011 29 / 34
DUNE Software Framework
BlueGene/P at J¨ ulich Supercomputing Center P · 803 degrees of freedom (51203 finest mesh), CCFV Clipped random permeability, σ2 = 8, λ = 4h, 10−8 reduction AMG used as preconditioner in BiCGStab (2 V-Cycles!) procs 1/h lev. TB TS It TIt TT 1 80 5 19.93 49.39 12 4.116 69.32 8 160 6 28.1 73.7 15 4.91 102 64 320 7 75.1 105 20 5.26 180 512 640 8 80.11 134 25 5.362 214.1 4096 1280 10 84.71 171.7 33 5.203 256.4 32768 2560 11 93.24 189.5 36 5.264 282.7 262144 5120 12 195.9 386.5 72 5.368 582.5
Compositional Two-Phase Flow Linz, October 6, 2011 30 / 34
DUNE Software Framework
4×12 AMD Magny Cours, 2.1 GHz, 12×0.5MB L2, 12MB L3. SIPG(k = 1) DG discretization of Poisson Problem 4 × 4 blocks, 160 × 160 × 90, 9.2 Mio degrees of freedom BiCGStab + inexact block Jacobi prec. (No AMG!) P #IT(max) Tit S Tass S 1 66 11.59
100 1.50 7.7 103.3 6.4 12 104 1.87 6.2 70.4 9.5 24 110 0.61 18.9 35.7 18.7 48 110 0.40 29.1 18.3 36.5
Compositional Two-Phase Flow Linz, October 6, 2011 31 / 34
DUNE Software Framework
Rapid prototyping: Substantially reduce time to implement discretizations and solvers for systems of PDEs based on DUNE. Simple things should be simple — suitable for teaching. Discrete function spaces:
◮ Conforming and non-conforming, ◮ hp-refinement, ◮ general approach to constraints, ◮ simple construction of product spaces for systems.
Operators based on weighted residual formulation:
◮ Linear and nonlinear, ◮ stationary and transient, ◮ FE and FV schemes requiring at most face-neighbors.
Exchangeable linear algebra backend. User only involved with “local” view on (reference) element.
Compositional Two-Phase Flow Linz, October 6, 2011 32 / 34
DUNE Software Framework
Problem Scheme LOC ∇ · {vu − k∇u} = f Pk, Qk 643 DG 1034 −∇ · {k∇u} = f CCFV 223 RT0 289 Mimetic 278 Two-phase flow in porous media pl, pc, CCFV 700 pl, pc, DG 931 CCFV, 2c 986 Stokes P2/P1, Q2/Q1 539 DG 1122 Linear acoustics DG 700 Maxwell in time domain DG 857 Nedelec 488
This does not include finite element spaces and problem setup.
Compositional Two-Phase Flow Linz, October 6, 2011 33 / 34
DUNE Software Framework
DG for two-phase flow
◮ Fully-coupled pl, pc formulation including media heterogeneities
in 3d.
◮ Saves about two levels of mesh refinment compared to low
◮ More efficient linear solvers are required.
Compositional two-phase flow
◮ pg, pc formulation handling phase disappearance ◮ First results on nuclear waste and CO2 problems
DUNE software framework
◮ Allows rapid implementation for wide variety of discretization
schemes for complex problems
◮ Use adaptive refinement, fast solvers and parallelization with
relative ease
Compositional Two-Phase Flow Linz, October 6, 2011 34 / 34