DG approximation of two-component miscible liquid-gas porous media - - PowerPoint PPT Presentation

Setting Numerical method Results

DG approximation of two-component miscible liquid-gas porous media flows

Alexandre Ern and Igor Mozolevski Universit´ e Paris-Est, CERMICS, Ecole des Ponts, France RICAM Workshop, Linz, October 2011

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Introduction

◮ Multicomponent multiphase porous media flows are encountered in

several applications

◮ petroleum engineering (oil-water) ◮ agricultural engineering, groundwater remediation (air-water)

◮ Such flows have received enhanced attention recently

◮ gas sequestration ◮ underground repositories of radioactive waste Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

The issue of radioactive waste I

◮ Production in France estimated at 2 kg/year/citizen

◮ various sources: electro-nuclear plants, health treatments, industrial

processes, research, army

◮ Total conditioned radioactive waste amounts to 1.2 Mm3 ◮ 0.2% of this volume contains 99% of radioactivity ◮ HAVL waste High Activity, Long Life 1 Myears

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

The issue of radioactive waste II

◮ HAVL waste could be stored in underground repository

◮ ANDRA (French Agency for Radioactive Waste Management) ◮ preliminary study established feasibility in 2005 ◮ clay host rock located at ≈ 500 m depth ◮ decision to be taken in 2014, operating could start in 2025 Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

The issue of radioactive waste III

◮ Underground research facility currently operating ◮ Various academic research programs have been launched ◮ GNR MOMAS Mathematical modeling and numerical

simulation www.gdrmomas.org

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

The issue of radioactive waste IV

◮ Multiple barriers to contain radionuclides (RN)

◮ conditioning of waste, steel containers ◮ manufactured barriers (bentonite, steel) ◮ host rock (quite favorable properties)

◮ Main time scales

◮ 102 years operating and observing the facility, reversibility ◮ 104 years degradation of manufactured barriers ◮ 106 years migration of RN through geosphere up to biosphere Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Hydrogen production and migration

◮ Corrosion of metallic components (and marginally water radiolysis)

Fe + 2H20 → Fe(OH)2 + H2

◮ Understand hydrogen migration through host rock ◮ Two-phase (liquid, gas), two-component (water, hydrogen) flow, gas

phase is compressible

◮ Gas phase (dis)appearance gas saturation at inflow position of gas phase front

2 4 6 8 10 x 10

5

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 t (year) Gas saturation at x=0 2 4 6 8 10 x 10

5

20 40 60 80 100 120 140 160 180 200 t (year) Saturation front position

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Outline

◮ Setting ◮ Numerical method ◮ Results

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Setting

◮ Governing equations ◮ Choice of main unknowns ◮ Mathematical model

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Governing equations I

◮ Basic notation

◮ subscript α ∈ {l, g} for phase ◮ superscript β ∈ {w, h} for component ◮ ̺β

α density of component β in phase α

◮ sα saturation of phase α, sl + sg = 1 ◮ pα pressure of phase α

◮ Mass conservation equation for each component

Φ

• α∈{l,g}

∂t(sα̺β

α) +

• α∈{l,g}

∇ · (̺β

αqα + jβ α) = F β

with Φ porosity, qα volumetric flow rate of phase α, jβ

α mass

diffusion flux of component β in phase α

◮ [Bear ’78, Chavent & Jaffr´

e ’78, Helmig ’97]

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Governing equations II

◮ Darcy–Muskat law for volumetric flow rates (neglecting gravity)

qα = −Kλα(sα)∇pα ∀α ∈ {l, g} with K absolute permeability, λα mobility of phase α (λα = 0 if phase α is absent)

◮ Capillary pressure π : [0, 1) → [0, +∞)

pg = pl + π(sg)

◮ Assume incompressibility in liquid phase and neglect water

vaporization ̺w

l = ̺std l

̺w

g = 0

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Governing equations III

◮ Ideal gas law in gas phase and hydrogen phase changes in

thermodynamic equilibrium (Henry’s law) ̺h

g = Cgpg

̺h

l = Chpg ◮ Fick’s law for dissolved hydrogen diffusion flux (dilute

approximation) jh

l = −ΦslDh l ∇̺h l

jw

l = −jh l ◮ Governing equations

Φ̺std

l

∂tsl + ∇ · (̺std

l

ql − jh

l ) = F w

Φ∂t(̺h

l sl + Cgpgsg) + ∇ · (̺h l ql + Cgpgqg + jh l ) = F h

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Choice of main unknowns I

◮ Possible absence of gas phase in a priori unknown parts of domain ◮ Choosing one of the saturations as one of the main unknowns is

inappropriate if gas phase disappears

◮ sg is identically 0 ◮ sl is identically 1

◮ Unified formulation of governing equations highly desirable to avoid

intricate numerical solvers

◮ Artificially enforcing sg ≥ ǫ > 0 is inappropriate (can lead to both

dissolved hydrogen and gas pressure overestimation)

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Choice of main unknowns II

◮ Following [Bourgeat, Jurak & Sma¨

ı ’09], we choose as main unknowns y = (y1, y2), y1 := pl, y2 := ̺h

l

allowing for a unified formulation including gas phase disappearance

pl ̺h

l

sg = ǫ ̺l

h = Chpl

sg = 0, no gas sg > 0, gas present ◮ See also

◮ [Jaffr´

e & Sboui ’10] for a reformulation based on complementary constraints

◮ [Abadpour & Panfilov ’09] for method with negative saturations Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Choice of main unknowns III

◮ Gas saturation recovered from capillarity and thermodynamic

equilibrium sg(pl, ̺h

l ) = sg(y) = π−1

y2 Ch − y1

• π−1 : R → [0, 1) inverse capillary pressure function extended by zero

1 ◮ sg is a continuous function of y

◮ continuously differentiable for van Genuchten capillary pressure

model

◮ not differentiable at entry pressure for Brooks–Corey model Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Mathematical model I

◮ Nondimensional form

◮ reference pressure p0 (1 MPa), reference density Chp0 (15g/m3) ◮ mass conservation equations scaled by ̺std

l

and Chp0

◮ Governing equations

∂tb1(y) − ∇ · (A11(y)∇y1 + A12(y)∇y2) = F1 ∂tb2(y) − ∇ · (A21(y)∇y1 + A22(y)∇y2) = F2 with bi(y) nondimensional and Aij(y) in m2/s

◮ IC on y, BC either Dirichlet on y or Neumann on total fluxes

σi(y) =

• j∈{1,2}

Aij(y)∇yj

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Mathematical model II

◮ System coefficients

b1(y) = −Φsg(y) b2(y) = Φa(sg(y))y2 A11(y) = p0Kλl(1 − sg(y)) A12(y) = −(Chp0/̺std

l

)Φ(1 − sg(y))Dh

l

A21(y) = y2p0Kλl(1 − sg(y)) A22(y) = y2p0Kλg(sg(y))ω + Φ(1 − sg(y))Dh

l

with a(s) = 1 + (ω − 1)s and ω = Cg/Ch (≈ 50)

◮ First equation is parabolic in y1, degenerating into elliptic if gas

phase disappears

◮ Second equation is parabolic in y2, degenerating into elliptic if

dissolved hydrogen disappears

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Mathematical model III

◮ Structure of space differential operator ◮ Ellipticity iff

(A12 + A21)2 < 4A11A22

◮ Typical orders of magnitude (in µm2/s)

A11 ≈ 50, A12 ≈ 0.7, A21 ≈ 50, A22 ≈ 450

◮ Under assumption A12 ≈ 0, ellipticity iff p0Kλl(1) < 4ΦDh l ◮ Alternative assumption for ellipticity is smallness condition on

hydrogen [Mikeli´ c ’09]

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Mathematical model IV

◮ Under assumption A12 ≈ 0, change of variables

y1 = u1 + ω−1eωu2 y2 = eωu2 yields coupled elliptic-parabolic system [Sma¨ ı ’09] ∂tb∗(u) − ∇ · (A∗(u)∇u) = F

◮ b∗ is the gradient of a convex potential ◮ A∗(u) is symmetric positive definite ◮ Fits Alt–Luckhaus theory for existence of weak solutions ◮ See also [Amaziane, Jurak & ˇ

Zgali´ c-Keko ’11; Khalil & Saad ’11] for mathematical analysis of compressible immiscible case

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Numerical method

◮ Time discretization: backward Euler ◮ Linearization: inexact Newton ◮ Space discretization: discontinuous Galerkin (dG)

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Time discretization

◮ Recall governing equations (i ∈ {1, 2})

∂tbi(y) −

• j∈{1,2}

∇ · (Aij(y)∇yj) = Fi

◮ Time discretization by backward Euler scheme (i ∈ {1, 2})

1 τ m (bi(y m) − bi(y m−1)) −

• j∈{1,2}

∇ ·

• Aij(y m)∇y m

j

• = F m

i

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Linearization

◮ Incomplete Newton method with fully coupled approach

◮ fixed-point for nonlinearities in diffusion operator ◮ linearization in time derivative

◮ For all m = 1, . . . , M (time loop) and for all l ≥ 0 (linearization

loop), find y m

l+1 s.t. (i ∈ {1, 2})

• j∈{1,2}

1 τ m ∂jbi(y m

l )y m j,l+1 −

• j∈{1,2}

∇ ·

• Aij(y m

l )∇y m j,l+1

• = G m

i,l

yielding a system of linear coupled PDEs in space

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Space discretization I

◮ dG methods can be viewed as

◮ FE-based methods using discrete functions with jumps ◮ FV-based high-order methods using numerical fluxes

◮ Vigorous development since the late 90s

◮ [Cockburn, Karniadakis & Shu ’00; Hesthaven & Warburton ’08] ◮ unified analysis Poisson problem [Arnold, Brezzi, Cockburn & Marini

’02], Friedrichs systems [AE & Guermond ’06]

◮ DG methods for two-phase immiscible porous media flows

◮ [Bastian ’99; Bastian & Rivi`

ere ’03; Eslinger ’05; Klieber & Rivi` ere ’06; Epshteyn & Rivi` ere ’07; AE, Mozolevski & Schuh ’09]

◮ Upcoming book Mathematical aspects of discontinuous

Galerkin methods, D. Di Pietro & AE, Springer, 2011

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Space discretization II

◮ {Tδ}δ>0 family of shape-regular meshes (possibly with hanging

nodes)

◮ Fi δ collects interfaces, Fb δ boundary faces, Fδ := Fb δ ∪ Fi δ ◮ Average operator {·}, jump operator [[·]], face normal nF ◮ For polynomial degree k ≥ 1,

V k

δ := {vδ ∈ L2(Ω); ∀T ∈ Tδ, vδ|T ∈ Pk(T)} ◮ Discrete pressure and hydrogen density both sought in V k δ

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Space discretization III

◮ Interior penalty dG bilinear form (i, j ∈ {1, 2})

aij

δ(yδ; uδ, vδ) =

• T∈Tδ
• T

Aij(yδ)∇uδ · ∇vδ −

• F∈Fi

δ∪FD δ

• F

nF·{Aij(yδ)∇uδ}[[vδ]] − θij

• F∈Fi

δ∪FD δ

• F

nF·{Aij(yδ)∇vδ}[[uδ]] +

• F∈Fi

δ∪FD δ

ηij

F

σk2 δF

• F

[[uδ]][[vδ]]

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Space discretization IV

◮ Discrete problem: For all vi,δ ∈ V k δ , i ∈ {1, 2},

• j∈{1,2}

1 τ m

• Ω

∂jbi(y m

δ,l)y m j,δ,l+1vi,δ +

• j∈{1,2}

aij

δ(y m δ,l; y m j,δ,l+1, vi,δ)

=

• Ω

G m

i,lvi,δ + bc’s ◮ Penalty and symmetry terms only on diagonal blocks (i = j) ◮ Penalty strategy reasonable if ellipticity can be asserted

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Space discretization V

◮ DG method weakly enforces

◮ zero PDE residual on each element T ∈ Tδ ◮ Dirichlet/Neumann BC’s on all boundary faces ◮ on all mesh interfaces F ∈ F i

δ, the transmission conditions

[[yi]] = 0, nF · [[σi(y)]] = 0, ∀i ∈ {1, 2}

◮ For heterogeneous media with different rock types

◮ weighted-averages can be considered ([AE, Stephansen & Zunino

’09] for linear convection-diffusion)

◮ for immiscible flows where gas saturation is one of the main

unknowns, nonlinear interface conditions must be enforced [AE, Mozolevski & Schuh ’09]

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Results

◮ Three test cases (1D)

◮ Gas-phase (dis)appearance [MOMAS benchmark] ◮ Ill-prepared initial condition [MOMAS benchmark] ◮ Heterogeneous medium [Bourgeat, Jurak & Sma¨

ı ’11]

◮ In all cases, we use

◮ van Genuchten model for capillary pressure and Mualem model for

relative permeability

◮ piecewise linear polynomials in dG method (k = 1) ◮ 10−8 convergence criterion in L2-norm for inexact Newton Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Gas-phase (dis)appearance I

◮ Domain Ω = (0, 200) (m) initially saturated by water ◮ Hydrogen injection at x = 0 during 105 years ◮ Injection rate qinj = 5.57·10−6 kg /m2/year ◮ Simulation time 106 years ◮ Initial and BC’s

−n · σ1|x=0 = 0 −n · σ2|x=0 = qinjχ[0,Tinj](t) pl|x=200 = 106 (Pa) ̺h

l |x=200 = 0

pl|t=0 = 106 (Pa) ̺h

l |t=0 = 0 ◮ Uniform mesh with 200 elements, time steps within [125;5000] years

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Gas-phase (dis)appearance II

Problem data Porous medium Fluid characteristics Param. Value Param. Value Φ 0.15 (-) Dh

l

3 10−9 m2/s K 5 10−20 m2 µl 1 10−3 Pa·s Pr 2 106 Pa µg 9 10−6 Pa·s n 1.49 (-) H(303K) 7.65 10−6 mol/Pa/m3 slr 0.4 (-) Mh 2 10−3 kg/mol sgr 0 (-) ̺std

l

103 kg/m3

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Gas-phase (dis)appearance III

◮ Snapshots for times in [0, 105] years liquid pressure (MPa)

• diss. hyd. molar density (mol/m3)

gas saturation (%)

50 100 150 200 1 1.05 1.1 1.15 x (m) Liquid pressure 2 103 1.4 104 105 1.4 105 5 105 50 100 150 200 2 4 6 8 10 12 x (m) Dissolved hydrogen molar density 5.1 105 5.6 105 6.1 105 6.6 105 106 50 100 150 200 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x (m) Gas saturation 2 103 1.4 104 105 1.4 105 5 105

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Gas-phase (dis)appearance IV

◮ Snapshots for times in [105, 106] years liquid pressure (MPa)

• diss. hyd. molar density (mol/m3)

gas saturation (%)

50 100 150 200 0.75 0.8 0.85 0.9 0.95 1 1.05 x (m) Liquid pressure 5.1 105 5.6 105 6.1 105 6.6 105 106 50 100 150 200 2 4 6 8 10 12 x (m) Dissolved hydrogen molar density 5.1 105 5.6 105 6.1 105 6.6 105 106 50 100 150 200 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x (m) Gas saturation 5.1 105 5.6 105 6.1 105 6.6 105 106

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Ill-prepared initial condition I

◮ Domain Ω = (0, 1) (m) with zero flux BC’s and no external sources ◮ Initial conditions

pl(x, 0) = 106 (Pa) x ∈ (0, 1) pg(x, 0) =

• 1.5·106 (Pa)

if x ∈ (0, 0.5) 2.5·106 (Pa) if x ∈ (0.5, 1) so that system is initially out of mechanical equilibrium

◮ Similar medium and fluid parameters (except K = 10−16 m2) as

previous case

◮ Simulation time 106 s ◮ Uniform mesh with 512 elements, time steps within [1;104] s

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Ill-prepared initial condition II

◮ Snapshots for times [0,103] s liquid pressure (MPa)

• diss. hyd. molar density (mol/m3)

gas saturation (%)

0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 x (m) Liquid pressure 10 100 500 1000 0.2 0.4 0.6 0.8 1 10 11 12 13 14 15 16 17 18 19 20 x (m) Dissolved hydrogen molar density 10 100 500 1000 0.2 0.4 0.6 0.8 1 2 4 6 8 10 12 14 16 18 x (m) Gas saturation 10 100 500 1000

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Ill-prepared initial condition III

◮ Snapshots for times [5·103,105] s liquid pressure (MPa)

• diss. hyd. molar density (mol/m3)

gas saturation (%)

0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 x (m) Liquid pressure 5 103 104 5 104 105 0.2 0.4 0.6 0.8 1 10 11 12 13 14 15 16 17 18 19 20 x (m) Dissolved hydrogen molar density 5 103 104 5 104 105 0.2 0.4 0.6 0.8 1 2 4 6 8 10 12 14 16 18 x (m) Gas saturation 5 103 104 5 104 105

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Ill-prepared initial condition IV

◮ Snapshots for times [2·105,106] s liquid pressure (MPa)

• diss. hyd. molar density (mol/m3)

gas saturation (%)

0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 x (m) Liquid pressure 2 105 5 105 106 0.2 0.4 0.6 0.8 1 10 11 12 13 14 15 16 17 18 19 20 x (m) Dissolved hydrogen molar density 2 105 5 105 106 0.2 0.4 0.6 0.8 1 2 4 6 8 10 12 14 16 18 x (m) Gas saturation 2 105 5 105 106

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Heterogeneous medium I

◮ Similar to test case 1 except that medium consists of two rock types

• ccupying Ω1 = (0, 20) and Ω2 = (20, 200) (m) and that hydrogen

injection is not stopped

◮ Rock occupying Ω2 has a finer texture ⇒ capillary barrier at

interface x = 20

◮ Uniform mesh with 160 elements (fitted to rock interface), time

steps within [200;20000] years

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Heterogeneous medium II

Snapshots at times {3·104, 4.2·104, 1.3·105, 106} years

50 100 150 200 1 1.05 1.1 1.15 1.2 1.25 x (m) Liquid pressure 3 104 4.2 104 1.3 105 106 50 100 150 200 2 4 6 8 10 12 14 16 x (m) Dissolved hydrogen molar density 3 104 4.2 104 1.3 105 106 50 100 150 200 10

−3

10

−2

10

−1

10 10

1

x (m) Gas saturation 3 104 4.2 104 1.3 105 106 50 100 150 200 0.2 0.4 0.6 0.8 1 x (m) Capillary pressure 3 104 4.2 104 1.3 105 106

Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows

Setting Numerical method Results

Conclusions

◮ Two-component miscible liquid-gas porous media flow model

motivated by hydrogen migration in underground radioactive waste repositories

◮ Choice of main unknowns to handle absence of gas phase ◮ DG space discretization combined with backward Euler scheme and

incomplete Newton linearization

◮ Further work

◮ deeper analysis of mathematical model ◮ multidimensional test cases with heterogeneities Alexandre Ern Universit´ e Paris-Est, CERMICS Two-component two-phase flows