Shuyu Sun

Earth Science and Engineering program KAUST Presented at the 2009 annual UTAM meeting, 2:05-2:40pm January 7, 2010 at the Sutton Building, University of Utah, Salt Lake City, Utah

Energy and Environment Problems

Single Phase Flow in Porous Media • Continuity equation – from mass conservation  

t

     

q m

, (

x

,

t

)   (0,

T

] • Thermodynamic model     (

T

,

P

) • For impressible fluid (constant density):   

u



q

(

x

,

t

)    (0,

T

] • Still need one more equation 

Relate velocity with pressure: Darcy's law

• Experiment by Henry Darcy (1855–1856)

Darcy's law

• Can be derived from the Navier-Stokes equations via homogenization. • It is analogous to – Fourier's law in the field of heat conduction, – Ohm's law in the field of electrical networks, – Fick's law in diffusion theory.

• In 3D:

 Incompressible Single Phase Flow • Continuity equation , ( • Darcy’s law ( ) , ] • Boundary conditions:

p



p B

u



n



u B

( (

x

,

t

)

x

,

t

)  

D

 

N

 (0,

T

 (0,

T

] ]

Transport in Porous Media • Transport equation 

t

  ( ) ( • Boundary conditions  uc 

D



c

   

D



c

 

n n

 

c B

u



n

0

t t

 (0,

T

],  (0,

T

],

x

 

in

(

t

)

x

 

out

(

t

)  • Initial condition (  0  

t

Numerical Methods for Flow & Transport • Challenge #1: Require the numerical method to be: – Locally conservative for the volume/mass of fluid (flow equation) – Locally conservative for the mass of species (transport equation) – Provides fluxes that is continuous in the normal direction across the entire domain.

• Methods that are not locally conservative without post-processing – Point-Centered Finite Difference Methods – Continuous Galerkin Finite Element Methods – Collocation methods – ……

Example: importance of local conservation

Example: importance of local conservation

Numerical Methods for Flow & Transport • Challenge #2: Fractured Porous Media – Different spatial scale: fracture much smaller – Different temporal scale: flow in fracture much faster • Solutions: – Mesh adaptation for spatial scale difference – Time step adaptation for temporal scale difference

Example: flow/transport in fractured media

Example: flow/transport in fractured media Locally refined mesh: FEM and FVM are better than FD for adaptive meshes and complex geometry

Example: flow/transport in fractured media CFL condition requires much smaller time step in fractures than in matrix: adaptive time stepping.

Numerical Methods for Flow & Transport • Challenge #3: Sharp fronts or shocks – Require a numerical method with little numerical diffusion – Especially important for nonlinearly coupled system, with sharp gradients or shocks easily being formed • Solutions: – Characteristic finite element methods – Discontinuous Galerkin methods

Example: Comparison of DG and FVM Upwind-FVM on 40 elements Linear DG on 40 elements Advection of an injected species from the left boundary under constant Darcy velocity. Plots show concentration profile at 0.5 PVI.

Example: Comparison of DG and FVM Flow in a medium with high permeability region (red) and low permeability region (blue) with flow rate specified on left boundary. Contaminated fluid flood into clean media.

Example: Comparison of DG and FVM FVM Linear DG Advection of an injected species from the left. Plots show concentration profiles at 3 years (0.6 PVI).

Numerical Method for Flow & Transport • Challenge #4: Time dependent local phenomena – For example: moving contaminant plume • Solutions: – Dynamic mesh adaptation • Based on conforming mesh adaptation • Based on non-conforming mesh adaptation

Adaptive DG methods – an example • Sorption occurs only in the lower half sub domain, • SIPG is used.

Adaptive DG example (cont.)

Adaptive DG example (cont.)

Adaptive DG example (cont.)

A Posteriori

Error Estimators • Residual based – L2(L2) – L2(H1) • Implicit – Solve a dual problem, can give estimates on a target functional – Disadvantages: computational costly and not flexible – Advantages: More accurate estimates • Hierarchical bases – Brute-force: difference between solutions of two discretizations (most expensive) – Local problems-based – Advantage: can guide anisotropic

hp

-adaptivity • Superconvergence points-based – Difficult for unstructured and non-conforming meshes

A posteriori

error estimates • Residuals – Interior residuals

I



t

 – (Element-)boundary residuals

R R

  

C



B

   

u D u n D



n

         

h h

, ,

Diri

,

out



A posteriori

error estimate in L2(L2) for SIPG

DG C



L

)     

E

4

h

4

r I

2

L

( )) 1 2 

r

  2

R

2  1 2 3 

r

2

R

(     

E

3

h

3

r

2

R

(  Proof Sketch: Compare with L2 projection; Cauchy Schwarz; Properties of cut-off operator; Approximation results; Inverse and Gronwell’s inequalities; Relation of residue and error

Dynamic mesh adaptation with DG • Nonconforming meshes – Effective implementation of mesh adaptation, – Elements will not degenerate unless using anisotropic refinement on purpose. • Dynamic mesh adaptation – Time slices = a number of time steps; only change mesh for time slices. – Refinement + coarsening  elements remain constant. number of

Concentration projections during dynamic mesh modification • Standard L2 projection used – Computation involved only in elements being coarsened • L2 projection is a local computation for discontinuous spaces – This results in computational efficiency for DG – L2 projection is a global computation for CG • L2 projection is locally mass conservative – This maintains solution accuracy for DG – Interpolation or interpolation-based projection used in CG is NOT locally conservative

Adaptive DG example (quads)

Adaptive DG example (quads)

Adaptive DG example (quads)

Adaptive DG (with triangles) Initial mesh

Adaptive DG (with triangles) T=0.5

T=1.0

Adaptive DG (with triangles) T=1.5

T=2.0

Adaptive DG example in 3D T=1.5

T=0.1

T=0.5

T=2.0

T=1.0

ANDRA-Couplex1 case Background – ANDRA: the French National Radioactive Waste Management Agency – Couplex1 Test Case • Nuclear waste management: Simplified 2D Far Field model • Flow, Advection, Diffusion-dispersion, Adsorption Challenges – Parameters are highly varying • permeability; retardation factor; effective porosity; effective diffusivity – Very concentrated nature of source • concentrated in space • concentrated in time – Long time simulation • 10 million years – Multiple space scales • Around source / Far from source – Multiple time scales • Short time behavior (Diffusion dominated) • Long time behavior (Advection dominated)

200k years ANDRA-Couplex1 case (cont.) 2m years

Compositional Three-Phase Flow

• Mass Conservation (without molecular diffusion)

U

i

   

w

,

o

• Darcy’s Law ,

c



g x i

, 

u



u

  

k r

  

K

 

P

   

g

 ,  

w

,

o

,

g

Numerical Modeling for Flow & Transport • Challenge #5: Importance of capillarity – Capillary pressure usually ignored in compositional flow modeling – Even the immiscible two-phase flow or the black oil model usually assumes only a single capillary function (i.e. assuming a single uniform rock)

Example: Reservoir Description • Two-dimensional 400x200m^2 domain • Contains a less-permeable (K=1md) rock in the center of the domain while the rest has K=100md. Isotropic permeability tensor used. • Porosity = 0.2 • Densities: 1000 kg/m^3 (W) and 660 kg/m^3 (O) • Viscosities: 1 cp (W) and 0.45 cp (O) • Inject on the left edge, and produce on the right edge • Injection rate: 0.1 PV/year • Initial water saturation: 0.0; Injected saturation: 1.0

Reservoir Description (cont.) • Relative permeabilities (assuming zero residual saturations):

k rw



S m we

,

k rn

  1 

S we



m

,

S we



S w

,

m

 2 • Capillary pressure

p c

(

S we

)  

B c

log

S we

,

S we



S w

,

B c

 5 and 50 bars K=100md K=1md

Discretization • DG-MFE-Iterative • Pressure time step: 10years / 1000 timeSteps • Saturation time step = 1/100 pressure time step • Mesh: 32x64 uniform rectangular grid:

Comparison: if ignore capillary pressure … With nonzero capPres With zero capPres Saturation at 10 years: Iter-DG-MFE

Numerical Modeling for Flow & Transport • Challenge #6: Discontinuous saturation distribution – Saturation usually is discontinuous across different rock type, which is ignored in many works in literature – When permeability changes, the capillary function usually also changes! • Solutions: – Discontinuous Galerkin methods

Saturation at 3 years Notice that Sw is continuous within each rock,

but Sw is discontinuous across the two rocks

Iter-DG-MFE Simulation

Saturation at 5 years Notice that Sw is continuous within each rock,

but Sw is discontinuous across the two rocks

Iter-DG-MFE Simulation

Saturation at 10 years Notice that Sw is continuous within each rock,

but Sw is discontinuous across the two rocks

Iter-DG-MFE Simulation

Water pressure at 10 years Notice that Pw is continuous within the entire domain.

Iter-DG-MFE Simulation (pressure unit: Pa)

Capillary pressure at 10 years Notice that Pc is continuous within the entire domain.

Iter-DG-MFE Simulation (pressure unit: Pa)

Numerical Modeling for Flow & Transport • Challenge #7:

Multiscale heterogeneous permeability

– Fine scale permeability has pronounced influence on coarse scale flow behaviors – Direct simulation on fine scale is intractable with available computational power • Solutions: – Upscaling schemes – Multiscale finite element methods

Recall: DG scheme for flow equation • Bilinear form    

e

   

e K

form  

e

    [

n

e w s

form   1    0 1 • Linear functional 

e e e

  

e

  0 0 SIPG IIPG NIPG NIPG

a h l

• Fine mesh DG on two meshes ( ( • Coarse mesh  ) (

Space decomposition • Introduce

V f

(   • Solution

p p p h H f

.



H

,

RH f



f

:

a H

 

f f

 

pv H

( )

RH vV f

Closure Assumption • Introduce

V f

0  :

f



H

• Two-scale solution 

p D H



H

   

H

.

( )

RH

Implementation • Multiscale basis functions: For each 

H T RH

(  0

a v v

,

H

• Multiscale approximation space:  :  • Two-scale DG solution  )

V

Other Closure Options • Local problems for solving multiscale basis functions need a closure assumption. • In previous derivation, we strongly impose zero Dirichlet boundary condition on local problems. • Other options: – Weakly impose zero Dirichlet boundary condition on local problems. – Strongly impose zero Neumann boundary condition on local problems.

– Weakly impose zero Neumann boundary condition on local problems.

– Combination of zero Neumann and zero Dirichlet.

Comparison with direct DG • Memory requirement – Direct DG solution in fine mesh: – Multiscale DG solution

r d O

(

h d

)

d R

( )

d H

(

d r d

(/ ) )

hH

• Computational time – Direct DG solution in fine mesh: 1 s f

d r d h

o – Multiscale DG solution

H O h

o

H

t

O H

• Conductivity: Example • Boundary conditions: – Left: p=1; Right: p=0; top & bottom: u=0. • Discretization: – R=r=1; – Coarse mesh 16x16; Fine mesh 256x256

Example (cont.) Coarse DG solution Brute-force Fine DG solution

Example (cont.) Multiscale DG solution Brute-force Fine DG solution

Future work • Multiscale DG methods for compositional multiple-phase flow in heterogeneous media, • Stochastic PDE simulations, • Multigrid solver for DG (including p-multigrid), • Other future works: – Automatically adaptive time stepping, – Implicit

a posteriori

error estimators, – Fully automatically

hp-

adaptivity for DG, –

A posteriori

and flow. estimators for coupled reactive transport