[PPT] - Fully Conservative Characteristic Methods for Flow and Transport: PowerPoint Presentation

SLIDE 1

Fully Conservative Characteristic Methods for Flow and Transport: Part II, Theoretical Considerations

Todd Arbogast Department of Mathematics

and

Center for Subsurface Modeling, Institute for Computational Engineering and Sciences (ICES) The University of Texas at Austin Wenhao Wang, The University of Texas at Austin

This work was supported by

U.S. National Science Foundation
KAUST through the Academic Excellence Alliance

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 1 The University of Texas at Austin, USA

SLIDE 2

Outline

1. The Volume Corrected Characteristics-Mixed Method
2. Monotonicity and Stability of VCCMM
3. Convergence of VCCMM
4. Numerical Convergence Tests
5. Some Numerical Examples
6. Extension to Compressible Flows
7. Summary and Conclusions

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 2 The University of Texas at Austin, USA

SLIDE 3

The Volume Corrected Characteristics-Mixed Method

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 3 The University of Texas at Austin, USA

SLIDE 4

Characteristic Trace-backs

Points. The characteristic trace-back of x is ˇ

x = ˇ x(x; t):

dˇ

x

dt = u(ˇ

x, t)

φ(ˇ

x)

, tn ≤ t < tn+1 ˇ

x(tn+1) = x

Assume u · ν = 0 so no particles trace to the domain boundary.

✲

tn tn+1

✻

x

ˇ

x

Space Time

✲

tn tn+1

✻

E ˇ E EE

Elements. Let Th be the collection of elements.

Particles in E ∈ Th trace back in space-time to EE = {(ˇ

x, t) ∈ Ω×[tn, tn+1] : ˇ x = ˇ x(x; t), x ∈ E}.

The “bottom” is ˇ E = {ˇ

x ∈ Ω : ˇ x = ˇ x(x; tn), x ∈ E}.

Let the pore volume of R be |R|φ =

R φ dx.

The local volume constraint: |E|φ = | ˇ E|φ +

EE

q dx dt

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 4 The University of Texas at Austin, USA

SLIDE 5

The VCCMM Approximate local volume constraint: Perturb The polygonal approximation of ˇ E is perturbed to ˜ E so that |E|φ = | ˜ E|φ +

˜

EE

q dx dt

② ② ② ② ② ② ② ②

ˇ E ≈ ˜ E Vol( ˇ E) = Vol( ˜ E)

② ② ② ② ② ② ② ②

Volume Corrected Characteristics Mixed Method (VCCMM): c(x, tn) ≈ cn

h(x) = cn E ∈ R

when x ∈ E For each E ∈ Th, define cn+1

E

by cn+1

E

|E|φ =

˜

E φcn h dx +

˜

EE

qcn

h dx dt Center for Subsurface Modeling Institute for Computational Engineering and Sciences 5 The University of Texas at Austin, USA

SLIDE 6

A Note on Updating Concentrations It is convenient to store the trace-back elements ˜ E as a polyline, i.e., a doubly linked list of vertices.

V1 V2 V3 V4 V5 V6 V7 V8 V1 V2 V3 V4 V5 V6 V7 V8 V9=V1 Polygon Polyline class

For each E ∈ Th, in the absence of sources, cn+1

E

|E|φ =

˜

E φcn h dx

=

F∈Th
˜

E∩F φcn F dx

=

F∈Th

cn

F| ˜

E ∩ F|φ ˜ E ∩ F is calculated efficiently by a 2-D Sutherland-Hodgman clipping algorithm (Foley, Van Dam, Feiner, & Hughes, 1995)

? E F ˜ E F ˜ E

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 6 The University of Texas at Austin, USA

SLIDE 7

An Algebraic View of VCCMM—1 VCCMM reduces to cn+1

E

= |E|−1

φ ˜ En cn h dx +

˜

En

E

qch dx dt

= |E|−1

φ F∈Th

cn

F

| ˜

En ∩ F|φ +

˜

En

E∩In F

q− dx dt

+
˜

En

E

cIq+ dx dt

=
F

An

E,F cn F + bn E,

where In

F is the space-time cylinder F × (tn, tn+1) over F.

Consider the piecewise constant function cn

h as a vector indexed by Th

cn = {cn

E}E∈Th

Thus

cn+1 = Ancn + bn

where bn

E = |E|−1 φ

˜

En

E

cIq+ dx dt.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 7 The University of Texas at Austin, USA

SLIDE 8

An Algebraic View of VCCMM—2 Case 1: E / ∈ Th,P (not at any production well, q− = 0 in ˜ En

E)

An

E,F = |E|−1 φ | ˜

En ∩ F|φ Case 2: E ∈ Th,P (a production well) Note that production wells expand during trace-back, so In

E ⊂ ˜

En

E, and q− = 0 for F

near E, F = E. Thus |E|φAn

E,F = | ˜

En ∩ F|φ +

˜

En

E∩In F

q− dx dt = (1 − δE,F)| ˜ En ∩ F|φ + δE,F

|E|φ +
In

E

q− dx dt

V n

E

where δE,F = 1 if E = F and 0 if E = F.

✲

tn tn+1

✻

E IE ˜ EE q < 0 q = 0 q = 0 E F Note: V n

E ≥ 0 is the volume of fluid remaining in E after production.

Problem: However, V n

E might be negative if production is large!

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 8 The University of Texas at Austin, USA

SLIDE 9

An Algebraic View of VCCMM—3 Modification when V n

E < 0:

1. Set V n

E = 0 in above formula.

2. Reduce the volume of fluid in each nearby F proportionately by

| ˜ En ∩ F| | ˜ En \ E|V n

E

Case 2 (Modified): E ∈ Th,P (a production well) |E|φAn

E,F =(1 − δE,F)

| ˜

En ∩ F|φ + | ˜ En ∩ F|φ | ˜ En \ E|φ (V n

E)−

if negative
+ δE,F (V n

E)+

if positive

Remarks:

This is consistent with the unmodified case when V n

E ≥ 0.

Note that the local volume constraint still holds.
In practice, we might ignore the production wells, since there fluid

leaves the domain (and we know how much by mass conservation).

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 9 The University of Texas at Austin, USA

SLIDE 10

An Algebraic View of VCCMM—4 VCCMM (Vector Form/Algebraic View):

cn+1 = Ancn + bn

Matrix An = (An

E,F)E,F∈Th and vector bn = (bn E)E∈Th are

An

E,F =

                          

| ˜ En ∩ F|φ |E|φ , E / ∈ Th,P, | ˜ En ∩ F|φ |E|φ

1 +

(V n

E)−

| ˜ En \ E|φ

,

E ∈ Th,P, F = E, (V n

E)+

|E|φ , E = F ∈ Th,P, bn

E =

1 |E|φ

˜

En

E

cIq+ dx dt, V n

E = |E|φ +

In

E

q− dx dt.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 10 The University of Texas at Austin, USA

SLIDE 11

Monotonicity and Stability of VCCMM

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 11 The University of Texas at Austin, USA

SLIDE 12

Properties of Matrix An—1 Lemma 1: Each An

E,F ≥ 0.

Proof: We need only consider E ∈ Th,P and F = E. Then An

E,F ∝ | ˜

En \ E|φ + (V n

E)− ≥ | ˜

En \ E|φ +

In

E

q− dx dt = | ˜ En|φ − |E|φ +

˜

En

E

q dx dt = 0

local volume constraint

.

Center for Subsurface Modeling

Institute for Computational Engineering and Sciences 12 The University of Texas at Austin, USA

SLIDE 13

Properties of Matrix An—2 Lemma 2: Each row sum

F∈Th

An

E,F ≤ 1.

Proof: We need only consider E ∈ Th,P and F = E. Then, since production wells expand, when E / ∈ Th,P,

F∈Th

An

E,F = | ˜

En|φ |E|φ = 1 |E|φ

|E|φ −
˜

En

E

q dx dt

≥0
≤ 1,

again using the local volume constraint. When E ∈ Th,P,

F∈Th

An

E,F = | ˜

En \ E|φ |E|φ

1 +

(V n

E)−

| ˜ En \ E|φ

+ (V n

E)+

|E|φ = 1 |E|φ (| ˜ En \ E|φ + V n

E) =

1 |E|φ

| ˜

En|φ +

In

E

q− dx dt

=

1 |E|φ

| ˜

En|φ +

˜

En

E

q dx dt

= 1
local volume constraint

.

Center for Subsurface Modeling

Institute for Computational Engineering and Sciences 13 The University of Texas at Austin, USA

SLIDE 14

Properties of Matrix An—3 Corollary: Matrix An does not increase the l∞-norm of a vector, i.e., |Anc|∞ ≤ |c|∞, where | · |∞ is the maximum component in absolute value: |x|∞ = max

E∈Th

|xE|. Proof: For E ∈ Th,

(Anc)E
=
F∈Th

An

E,FcF

≤
F∈Th

An

E,F|cF| ≤

F∈Th

An

E,F|c|∞ ≤ |c|∞.

Center for Subsurface Modeling

Institute for Computational Engineering and Sciences 14 The University of Texas at Austin, USA

SLIDE 15

Monotonicity of VCCMM VCCMM satisfies a maximum and minimum principle. Theorem: Up to the effect of injection wells, VCCMM can produce neither overshoots nor undershoots. Proof: Without bn (i.e., injection wells),

cn+1 = Ancn,

and so each cn+1

E

is a positive average of previous values cn

F.

Center for Subsurface Modeling

Institute for Computational Engineering and Sciences 15 The University of Texas at Austin, USA

SLIDE 16

Stability of VCCMM Theorem: VCCMM is stable; that is, if cn and ˜

cn solve cn+1 = Ancn + bn

(scheme), ˜

cn+1 = An˜ cn + bn + ∆tnδn

(perturbed scheme), then the error over time [0, T] is bounded by the initial error and the perturbation, i.e., max

0≤n≤N |˜

cn − cn|∞ ≤ |˜ c0 − c0|∞ + T

max

0≤n≤N |δn|∞

Proof: Taking differences and l∞-norms, |˜

cn+1 − cn+1|∞ ≤ |An(˜ cn − cn)|∞ + ∆tn|δn|∞

≤ |˜

cn − cn|∞ + ∆tn max

0≤n≤N |δn|∞.

Interating on n completes the proof.

Center for Subsurface Modeling

Institute for Computational Engineering and Sciences 16 The University of Texas at Austin, USA

SLIDE 17

Convergence of VCCMM

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 17 The University of Texas at Austin, USA

SLIDE 18

The Convergence Result Definition: The L1 and L∞ norms of f(x) over the domain Ω are fL1 =

Ω |f(x)| dx

and fL∞ = ess max

x∈Ω

|f(x)| Theorem: (T. Arbogast and W. Wang, 2009) For time step ∆t and maximal element diameter h, Eh = max

0≤n≤N cn − cn hL1 ≤ C

h

√ ∆t + h + (∆t)r

where r ≥ 1 is the order of approximately solving for ˇ

x.

Remarks:

1. The L1-norm is better than the L2-norm, since fluid volume and

tracer mass are more like L1 quantities.

2. There is no CFL constraint on ∆t !
3. Optimally, ∆t = O(h2/(2r+1)), and Eh = O(h2r/(2r+1)).
4. In practice, ∆t = O(h), and Eh = O(

√ h).

5. For Godunov: ∆t ≤

h 2uL∞ ≡ ∆tCFL and Eh = O( √ h).

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 18 The University of Texas at Austin, USA

SLIDE 19

Perturbed Velocity Key idea 1: There is a perturbed velocity ˜

u ≈ u such that E traces back

to ˜

E. That is, the transport is exact with the perturbed velocity! (cf.

Arbogast and Wheeler, 1995) Reasons for the perturbation:

1. Approximately solve for the characteristics ˇ

x.

2. Polygonal approximation of ˇ

E.

3. Point adjustment.

Lemma: u − ˜

uL∞ + ∇ · (u − ˜ u)L∞ ≤ C

h + (∆t)r

, where r ≥ 1 is the order of approximately solving for ˇ

x.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 19 The University of Texas at Austin, USA

SLIDE 20

Sketch of Proof of Perturbed Velocity There are three main steps to the local construction of ˜

u.

1. Perturbation in the characteristic direction (layer adjustment):

˜ x(x, t) = ˇ

x

x, t + α(tn+1 − t)
,

= ⇒

˜ u(x, t) = (1 − α) u

x, t + α(tn+1 − t)
.

Need to show α = O(h).

2. Perturbation in normal direction (midpoint adjustment):

˜ x(x, t) = ˇ

x(x, t) + s

tn+1 − t

∆tn

ν,

= ⇒

˜ u(x, t) = u

x − s
tn+1 − t

∆tn

ν, t
−

s ∆tn ν. Need to show s = O(h∆t).

3. Interpolation of the local definition of ˜

u gives the perturbed velocity

field in all of Ω × [0, T].

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 20 The University of Texas at Austin, USA

SLIDE 21

VCCMM as Exact Tracing Interpretation of VCCMM: Let ξ solve the perturbed system (φξ)t + ∇ · (ξ˜

u) = qξ

in Ω × [tn, tn+1], ξ(·, tn) = cn

h

at Ω × {tn}. Then cn+1

h

= Pξ(·, tn+1), where P is a projection into constants over the elements: Pξ(·, tn+1)

E = 1

|E|

E ξ(x) dx.

Remark: There are only two sources of error:

1. Using ˜

u instead of u;

2. Projection onto constants each time step.

The second error limits accuracy to O( √ h) when ∆t ∼ h. It can be reduced by post-processing onto linears, combined with a slope limiting procedure to maintain monotonicity properties.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 21 The University of Texas at Austin, USA

SLIDE 22

Approximation of L1-errors Key idea 2: For ǫ > 0, approximate the L1-error (cf. Kuznetsov, 1976) ρn

ǫ,h =

Ω×Ω Kǫ(x − y) |cn(x) − cn

h(y)| dx dy

where Kǫ(x) = ǫ−dK0(x/ǫ) is an approximation to the identity in Ω. Properties of ρn

ǫ,h:

1. First order approximation of L1-errors:
ρn

ǫ,h − cn − cn hL1

≤ Cǫ.
2. Change of ρn

ǫ,h in time: (derived from the entropy inequality)

ρn+1−

ǫ,h

− ρn

ǫ,h ≤ C

ǫ + h + (∆t)r

∆tn.

3. Projection error: (requires analysis in the space BV)

ρn

ǫ,h − ρn− ǫ,h ≤ Ch2

ǫ |cn−

h |BV (Ω).

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 22 The University of Texas at Austin, USA

SLIDE 23

Sketch of Convergence Proof Proof of Convergence: Summing changes of ρn

ǫ,h in each time steps gives

ρn

ǫ,h ≤ ρ0 ǫ,h + En ǫ,h + O(ǫ + h + (∆t)r),

where the total projection error is En

ǫ,h := n

k=1

(ρk

ǫ,h − ρk− ǫ,h) = O

h2

ǫ∆t

by summing the projection error in each time step.

Replacing ρn

ǫ,h and ρ0 ǫ,h with errors cn h − cnL1 + O(ǫ) and

c0

h − c0L1 + O(ǫ) gives

cn

h − cnL1 ≤ c0 h − c0L1 + C

ǫ + h2

ǫ∆t + h + (∆t)r

,

where the optimal choice is to take ǫ = h/ √ ∆t.

Center for Subsurface Modeling

Institute for Computational Engineering and Sciences 23 The University of Texas at Austin, USA

SLIDE 24

Numerical Convergence Tests

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 24 The University of Texas at Austin, USA

SLIDE 25

Convergence Test 1 A quarter 5-spot problem

1. Field Parameters
Ω = 15m×20m
k = 10 md
u · ν = 0 on ∂Ω
Rate q = 1.2 m2/min
2. Simulation Parameters
Initial c0 = 0
Inject cI = 1
3. Discretization
Uniform grid n × n

5 10 15 5 10 15 20 Center for Subsurface Modeling Institute for Computational Engineering and Sciences 25 The University of Texas at Austin, USA

SLIDE 26

Errors for Convergence Test 1 “Exact” solution: Higher Order Godunov’s Method on 256 × 256 grid with ∆tCFL,256 ≈ 0.2289 sec in T = 1 hour. Numerical solution: VCCMM on n × n grid with optimal ∆t = Ch2/3 (Using Euler to solve for characteristics, so r = 1). The error is Eh = 1 |Ω| max

0≤n≤N cn − cn hL1 ?

≤ Ch2/3 Define the ratio Ch = Eh h2/3 n h (m) ∆t (sec) Eh Ch 8 3.1250 115.34577 0.468699 0.219278 16 1.5625 72.66328 0.191371 0.142123 32 0.7813 45.77500 0.078365 0.092384 64 0.3907 28.83644 0.035072 0.065632 128 0.1953 18.16582 0.017668 0.052486 256 0.0977 11.44375 0.008821 0.041596 Conclusion: Ch is indeed bounded!

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 26 The University of Texas at Austin, USA

SLIDE 27

Convergence Test 2 A constant velocity field

1. Field Parameters
Ω = 15m×20m
u = (0.03, 0.04) m/sec
2. Simulation Parameters
Initial c0 = 0
c = 1 on Γin × [0, T]
No source q = 0
T = 500 sec
3. Discretization
Uniform grid n × n
4. No point adjustment is

needed in this test.

u c=0 c=1 Γin

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 27 The University of Texas at Austin, USA

SLIDE 28

Errors for Convergence Test 2 Exact solution: ct + u · ∇c = 0 in Ω × [0, T], c(x, t) = 1

n Γin × [0, T],

c0(x) = 0 in Ω. VCCMM: Expect that Eh = 1 |Ω| max

0≤n≤N cn − cn hL1 ?

≤ Ch2/3, so let Ch = Eh h2/3 n h (m) ∆t (sec) E′

h

Ch 8 3.1250 166.67 0.072524 0.033930 16 1.5625 105.00 0.045082 0.033480 32 0.7813 66.14 0.025061 0.029544 64 0.3906 41.67 0.016387 0.030666 128 0.1953 26.25 0.009717 0.028866 256 0.0977 16.54 0.006700 0.031596 512 0.0488 10.42 0.003963 0.029664 1024 0.0244 6.56 0.002621 0.031141 Conclusion: Ch is indeed constant!

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 28 The University of Texas at Austin, USA

SLIDE 29

Some Numerical Examples

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 29 The University of Texas at Austin, USA

SLIDE 30

Rotating Pollutant Problem—1

1. Field Parameters
u(x) = (−x2, x1) in R2
q = 0 in R2
2. Simulation Parameters
Initial c0(x) = 1

in [9, 10] × [0, 1], and 0 elsewhere

T = 2π
3. Discretization
Uniform grid on

integers in R2

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 30 The University of Texas at Austin, USA

SLIDE 31

Rotating Pollutant Problem—2 Pollutant Concentrations at T = 2π 5 Time Steps 10 Time Steps

30
20
10

10 20 30

30
20
10

10 20 30

C

1.0E+00 1.0E-05 1.0E-10 1.0E-15 1.0E-20 1.0E-25 1.0E-30 1.0E-35

30
20
10

10 20 30

30
20
10

10 20 30

C

1.0E+00 1.0E-05 1.0E-10 1.0E-15 1.0E-20 1.0E-25 1.0E-30 1.0E-35

20 Time Steps 40 Time Steps

30
20
10

10 20 30

30
20
10

10 20 30

C

1.0E+00 1.0E-05 1.0E-10 1.0E-15 1.0E-20 1.0E-25 1.0E-30 1.0E-35

30
20
10

10 20 30

30
20
10

10 20 30

C

1.0E+00 1.0E-05 1.0E-10 1.0E-15 1.0E-20 1.0E-25 1.0E-30 1.0E-35

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 31 The University of Texas at Austin, USA

SLIDE 32

Quarter Five-Spot Problem—1

1. Field Parameters
Ω = 15m×20m
Heterogeneous k:

Mk = 10 md Cv = 2.5

u · ν = 0 on ∂Ω
Rate q = 1.2 m2/min
2. Simulation Parameters
Initial clean c0 = 0
Inject cI = 1
3. Discretization
Uniform grid n × n

5 10 15 5 10 15 20

k

500 100 50 10 5 1 0.5 0.1 0.05 0.01

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 32 The University of Texas at Austin, USA

SLIDE 33

Quarter Five-Spot Problem—2 Comparison of CMM and VCCMM Concentrations at 100 min (Grid 100 × 100, ∆t = 30 sec) CMM VCCMM

5 10 15 5 10 15 20

C

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

5 10 15 5 10 15 20

C

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

(uncorrected) (corrected)

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 33 The University of Texas at Austin, USA

SLIDE 34

Quarter Five-Spot Problem—3 Comparison of Godunov and VCCMM Concentrations at 100 min (Grid 100 × 100) Godunov (∆tCFL = 1.5 sec) VCCMM (∆t = 30 sec = 20∆tCFL)

5 10 15 5 10 15 20

C

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

5 10 15 5 10 15 20

C

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 34 The University of Texas at Austin, USA

SLIDE 35

Nuclear Waste Disposal—1 Problem: Concentration of Iodine-129 leaking from a repository satisfies φR(ct + λc) + ∇ · (cu − D∇c) = 0 in Ω × [0, T] with initial leak concentration c0 = 0.133 mol/m3.

5000 10000 15000 20000 25000 100 200 300 400 500 600

uy = 0 ux = 0

Marl K = 3.2e-5

ux = 0

p = 200 Limestone K = 6.3 p = 310

ux = 0

Clay K = 3.2e-6 Repository

ux = 0

p = 286 Dogger K = 25.2 p = 289

uy = 0

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 35 The University of Texas at Austin, USA

SLIDE 36

Nuclear Waste Disposal—2 Grid mesh (non-uniform 108 × 70)

5000 10000 15000 20000 25000 100 200 300 400 500 600

Characteristic trace-back mesh

5000 10000 15000 20000 25000 100 200 300 400 500 600

Remark: Only a little work is required for volume adjustment.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 36 The University of Texas at Austin, USA

SLIDE 37

Nuclear Waste Disposal—3 Flow Approximation Pressure

5000 10000 15000 20000 25000 100 200 300 400 500 600

pressure

300 290 280 270 260 250 240 230 220 210

Speed (magnitude of velocity |u|)

5000 10000 15000 20000 25000 100 200 300 400 500 600

Velocity Magnitu

1.0E-02 1.0E-03 1.0E-04 1.0E-05 1.0E-06 1.0E-07 1.0E-08 1.0E-09

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 37 The University of Texas at Austin, USA

SLIDE 38