Stochastic multiscale modeling of subsurface and surface flows. Part - - PowerPoint PPT Presentation

Stochastic multiscale modeling of subsurface and surface flows. Part III: Multiscale mortar finite elements for coupled Stokes-Darcy flows

Ivan Yotov

Department of Mathematics, University of Pittsburgh

KAUST WEP Workshop January 30–February 1, 2010

Joint work with Vivette Girault, Paris VI, and Danail Vassilev, University of Pittsburgh

Department of Mathematics, University of Pittsburgh 1

Outline

• Mathematical model for the coupled flow problem

– Interface conditions – Existence and uniqueness for a global weak formulation – Equivalence to a domain decomposition weak formulation

• Finite element discretizations

– Conforming Stokes elements and mixed finite elements for Darcy – Mortar finite elements on all subdomain interfaces – Discrete inf-sup condition – Existence and uniqueness of a discrete solution – Convergence analysis

• Non-overlapping domain decomposition - reduction to a mortar interface problem
• Coupling of Stokes-Darcy flow with transport

Department of Mathematics, University of Pittsburgh 2

Applications of coupled Stokes and Darcy flows

• Interactions between surface water and groundwater flows
• Flows in fractured porous media
• Flows through vuggy rocks
• Flows through industrial filters
• Fuel cells
• Blood flows
• Glaucoma

Department of Mathematics, University of Pittsburgh 3

Model problem

Γ1 Γ1 Γ1 Γ2 Γ2 Γ2 Ω Ω

1 2

= fluid region (Stokes flow) = saturated porous medium (Darcy’s law) n n

2 1

ΓI

uj : Ωj → Rd, fluid velocity in Ωj, pj : Ωj → R , fluid pressure in Ωj. Deformation tensor D and stress tensor T in Ωs: D(us) := 1 2(∇us + (∇us)T), T(us, ps) := −psI + 2µD(us).

Department of Mathematics, University of Pittsburgh 4

Flow equations

Stokes flow on Ωs:        −div T(us, ps) ≡ 2µ div D(us) + ∇ps = fs in Ω1 (conservation of momentum), div us = 0 in Ωs (conservation of mass), us = 0

• n Γs

(no slip). Darcy flow on Ωd:        µK−1ud + ∇pd = fd in Ωd (Darcy’s law), div ud = qd in Ωd (conservation of mass), ud · nd = 0

• n Γd

(no flow), Solvability condition:

• Ωd

qd dx = 0

Department of Mathematics, University of Pittsburgh 5

Interface conditions

Mass conservation across Γsd: us · ns + ud · nd = 0 on Γsd. Continuity of normal stress on Γd: −ns · T · ns ≡ ps − 2µns · D(us) · ns = pd on Γsd. Slip with friction interface condition: (Beavers-Joseph (1967),Saffman (1971), Jones (1973), J¨ ager and Mikeli´ c (2000)) −ns · T · τ j ≡ −2µns · D(us) · τ j = µα

• Kj

u1 · τ j, j = 1, d − 1, on Γsd, where Kj = τ j · K · τ j.

Department of Mathematics, University of Pittsburgh 6

Global variational formulation

For ud, vd in L2(Ωd)d and us, vs in H1(Ωs)d, a(u, v) = µ

• Ωd

K−1ud·vd+2 µ

• Ωs

D(us) : D(vs)+

d−1

• j=1
• Γsd

µα

• Kj

(us·τ j)(vs·τ j). For vd ∈ H(div; Ωd), vs ∈ H(div; Ωs) and q ∈ L2(Ω), b(v, q) = −

• Ωd

q div vd −

• Ωs

q div vs. X = {v ∈ H(div; Ω) ; vs ∈ H1(Ωd)d, v|Γs = 0, (v · n)|Γd = 0}, ∀v ∈ X , vX =

• v2

H(div;Ω) + |vs|2 H1(Ωs)

1/2. Find (u, p) ∈ X × L2

0(Ω) such that

∀v ∈ X , a(u, v) + b(v, p) =

• Ω

f · v, ∀r ∈ L2

0(Ω) , b(u, r) = −

• Ω

r qd.

Department of Mathematics, University of Pittsburgh 7

Existence and uniqueness of a weak solution

Lemma: ∀q ∈ L2

0(Ω) , sup v∈X

b(v, q) vX ≥ βqL2(Ω) Lemma: a(v, v) ≥ γv2

X

∀ v ∈ V = {v ∈ X : div = 0} Proof: By Korn’s inequality, a(v, v) ≥ C(|vs|2

H1(Ωs) + vd2 L2(Ωd)).

Coercivity now follows from ∀v ∈ V , vsL2(Ωs) ≤ C

• |vs|2

H1(Ωs) + vd2 L2(Ωd)

1/2. Lemma: The variational problem has a unique solution.

Department of Mathematics, University of Pittsburgh 8

Domain decomposition variational formulation

Let Ωs = ∪Ωs,i, Ωd = ∪Ωd,i. Interface conditions Stokes-Stokes interfaces: [vs] = 0, [T · n] = 0 on Γss Darcy-Darcy interfaces: [ud · n] = 0, [pd] = 0 on Γdd ˜ X = {v ∈ L2(Ω)d : v|Ωs,i ∈ H1(Ωs,i)d, v|Ωd,i ∈ H(div, Ωd,i), v · n ∈ Lr(∂Ωd,i)(r > 1), v|Γs = 0, (v · n)|Γd = 0}, W = L2

0(Ω),

Λ = {ξ : ξ|Γij ∈ H−1/2(Γij) ∀ Γij ⊂ Γss, ξ|Γij ∈ H1/2(Γij) ∀ Γij ⊂ Γdd ∪ Γsd}.

Department of Mathematics, University of Pittsburgh 9

Domain decomposition variational formulation, cont.

as,i(us,i, vs,i) = 2 µ

• Ωs,i

D(us,i) : D(vs,i)+

d−1

• j=1
• ∂Ωs,i∩Γsd

µα

• Kj

(us,i·τ j)(vs,i·τ j) ad,i(ud,i, vd,i) = µ

• Ωd,i

K−1ud,i · vd,i, bi(vi, wi) = −

• Ωi

wi div vi a(·, ·) =

• as,i(·, ·) +
• ad,i(·, ·),

b(·, ·) =

• bi(·, ·)

bΓ(v, λ) =

• Γss

[v] λ +

• Γdd

[v · n] λ +

• Γsd

[v · n] λ Find u ∈ ˜ X, p ∈ W, λ ∈ Λ: ∀v ∈ ˜ X, a(u, v) + b(v, p) + bΓ(v, λ) =

• Ω

f · v, ∀w ∈ W, b(u, w) = −

• Ω

w qd, ∀ξ ∈ Λ, bΓ(u, ξ) = 0.

Department of Mathematics, University of Pittsburgh 10

Equivalence

Lemma: The two variational formulations are equivalent.

Department of Mathematics, University of Pittsburgh 11

Finite element discretization

Partition T h

i on Ωi; T h i and T h j need not match at Γij.

Stokes elements Xh

s,i × W h 1

in Ωs,i: MINI (Arnold-Brezzi-Fortin), Taylor- Hood, Bernardi-Raugel; contain at least polynomials of degree rs and rs − 1 resp.

Pressure Velocity

Mixed finite elements Xh

d,i × W h d,i in Ωd,i:

RT, BDM, BDFM, BDDF; contain at least polynomials of degree rd

velocity Lagrange pressure multplier

Xh :=

• Xh

i , W h :=

• W h

i ∩ L2 0(Ω)

T H

ij - partition of Γij, possibly different from the traces of T h i and T h j

ΛH

ij: continuous or discontinuous piecewise polynomials of degree at least ms on

Γss or md on Γdd and Γsd ΛH :=

• ΛH

ij

Nonconforming approximation: Λh ⊂ Λ

Department of Mathematics, University of Pittsburgh 12

Mortar finite element method

Find uh ∈ Xh, ph ∈ W h, λH ∈ ΛH: ∀vh ∈ Xh, a(uh, vh) + b(vh, ph) + bΓ(vh, λH) =

• Ω

f · vh, ∀wh ∈ W h, b(uh, wh) = −

• Ω

wh qd, ∀ξH ∈ ΛH, bΓ(uh, ξH) = 0. Equivalently, letting V h = {vh ∈ Xh : bΓ(vh, ξH) = 0 ∀ ξH ∈ ΛH}, Find uh ∈ V h, ph ∈ W h: ∀vh ∈ V h, a(uh, vh) + b(vh, ph) =

• Ω

f · vh, ∀wh ∈ W h, b(uh, wh) = −

• Ω

wh qd.

Department of Mathematics, University of Pittsburgh 13

Mortar compatibility conditions

On Γij ⊂ Γdd ∪ Γsd, i < j (assume that Ωi is a Darcy domain), ∀ ξH ∈ ΛH, sup

ϕh

i ∈V h i ·n

ϕh

i , ξHΓij

ϕh

i L2(Γij)

≥ βdξHL2(Γij). On Γij ⊂ Γss, i < j, ∀ ξH ∈ ΛH, sup ϕh

i ∈V h i |Γij∩H1/2 00 (Γij)d

ϕh

i , ξHΓij

ϕh

i H1/2(Γij)

≥ βsξHH−1/2(Γij). These are much more general than the mortar conditions used in [Bernardi-Maday- Patera], [Ben Belgacem] for Laplace and Stokes, and in [Layton-Schieweck-Y.], [Riviere-Y.], [Galvis-Sarkis] for Stokes-Darcy. The above condition on Γdd is similar to the one used in [Arbogast-Cowsar-Wheeler- Y.] and [Arbogast-Pencheva-Wheeler-Y.].

Department of Mathematics, University of Pittsburgh 14

Weakly continuous interpolants

V h

s = {vh ∈ Xh s : [vh], ξHΓss = 0 ∀ ξH ∈ ΛH}

V h

d = {vh ∈ Xh d : [vh · n], ξHΓdd = 0 ∀ ξH ∈ ΛH}

There exisit Πh

s : H1(Ωs) → V h d such that

(div(v − Πh

sv), wh)Ωs,i = 0 ∀wh ∈ W h, ∀Ωs,i ⊂ Ωs,

• Πh

svH1(Ωs,i) ≤ C

• vH1(Ωs,i)

There exist Πh

d : H1(Ωd) → V h s such that

(div(v − Πh

dv), wh)Ωd,i = 0 ∀wh ∈ W h, ∀Ωd,i ⊂ Ωd,

• Πh

dvH(div;Ωs,i) ≤ C

• vH1(Ωs,i)

Department of Mathematics, University of Pittsburgh 15

Discrete inf-sup condition

Lemma: ∀wh ∈ W h, sup

vh∈Vh

b(vh, wh) vX ≥ βwW Proof: Given wh ∈ W h, construct v ∈ H1(Ω): div v = −wh, vH1(Ω) ≤ CwhL2(Ω). Then construct an interpolant Πh : H1(Ω) → V h such that b(Πhv − v, qh) = 0 ∀ qh ∈ W h, ΠhvX ≤ CvH1(Ω).

Department of Mathematics, University of Pittsburgh 16

Existence and uniqueness of a discrete solution

Lemma: There exists a unique solution to: Find uh ∈ V h, ph ∈ W h: ∀vh ∈ V h, a(uh, vh) + b(vh, ph) = (f, vh), ∀wh ∈ W h, b(uh, wh) = −(qd, wh). Proof: Coercivity: The discrete Korn’s inequality gives as,i(vh, vh) ≥ αs,i|vh|2

H1(Ωs,i),

which, combined with ∀vh ∈ V h,

• vh2

L2(Ωs,i) ≤

• |vh|2

H1(Ωs,i),

implies

• as,i(vh, vh) ≥ αs
• vh2

H1(Ωs,i).

On Ωd,

• ad,i(vh, vh) ≥ αd
• vh2

H(div;Ωd,i) ∀vh : (div vh, wh) = 0 ∀wh ∈ W h.

Coercivity and inf-sup condition imply existence and uniqueness of a solution.

Department of Mathematics, University of Pittsburgh 17

Approximation properties of V h

Construct an interpolant Ih : H1(Ω) → V h, Ihv = (Ih

s v, Ih d v),

• v − Ih

s vH1(Ωs,i) ≤ C

• hrsvHrs+1(Ωs,i),
• v − Ih

d vH(div;Ωs,i) ≤ C

• (hrd+1vHrd+1(Ωd,i) + hrsvHrs+1(Ωs,i)).

Department of Mathematics, University of Pittsburgh 18

Convergence analysis

Theorem: u − uhX + p − phW ≤ C(hrs + hrd+1 + Hms+1/2 + Hmd+1/2) Proof: u − uhX + p − phW ≤≤ C( inf

vh∈V h u − vhX +

inf

wh∈W h p − whW) + Rh,

where Rh = sup

vh∈V h

|a(u, vh) + b(vh, p) − (f, vh)| vhX Bound on the non-conforming error: Θ(vh) := a(u, vh) + b(vh, p) − (f, vh) = −bΓ(vh, λ) |[vh · n], λΓsd| ≤ CHmd+1vhX |[vh · n], λΓdd| ≤ CHmd+1/2vhX |[vh], λΓss| ≤ CHms+1/2vhX

• Department of Mathematics, University of Pittsburgh

19

Example 1: Smooth solution

Ω = Ω1 ∪ Ω2, where Ω1 = [0, 1] × [1

2, 1] and Ω2 = [0, 1] × [0, 1 2]

Taylor-Hood on triangles in Ω1; lowest order Raviart-Thomas on rectangules in Ω2. u1 = u2 =

• sin( x

G + ω)ey/G

−cos( x

G + ω)ey/G

• ,

p1 = (G K − µ G)cos( x G + ω)e1/(2G) + y − 0.5, p2 = G Kcos( x G + ω)ey/G, where µ = 0.1, K = 1, α = 0.5, G = √µK α , and ω = 1.05.

Department of Mathematics, University of Pittsburgh 20

Convergence rates for Example 1

mesh u1 − u1,h1,Ω1 rate p − p1,h0,Ω1 rate 4x4 9.74e-02 3.73e-03 8x8 2.41e-02 2.01 9.62e-04 1.96 16x16 6.15e-03 1.97 2.58e-04 1.90 32x32 1.59e-03 1.95 6.90e-05 1.90 64x64 4.15e-04 1.94 1.84e-05 1.91 Table 1: Numerical errors and convergence rates in Ω1. mesh u2 − u2,h0,Ω2 rate p − p2,h0,Ω2 rate 4x4 3.18e-02 8.10e-03 8x8 9.20e-03 1.79 2.47e-03 1.71 16x16 2.48e-03 1.89 6.56e-04 1.91 32x32 6.56e-04 1.92 1.67e-04 1.97 64x64 1.72e-04 1.93 4.19e-05 1.99 Table 2: Numerical errors and convergence rates in Ω2.

Department of Mathematics, University of Pittsburgh 21

Example 2: discontinuous tangential velocity

u1 =

• (2 − x)(1.5 − y)(y − ξ)

−y3

3 + y2 2 (ξ + 1.5) − 1.5ξy − 0.5 + sin(ωx)

• ,

u2 =

• ω cos(ωx)y

χ(y + 0.5) + sin(ωx)

• ,

p1 = −sin(ωx) + χ 2K + µ(0.5 − ξ) + cos(πy), p2 = − χ K (y + 0.5)2 2 − sin(ωx)y K , where ω = 6 and the other parameters are defined as in Example 1.

Department of Mathematics, University of Pittsburgh 22

Computed pressure and velocity for Example 2

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

U 2.8 2.4 2 1.6 1.2 0.8 0.4

• 0.4
• 0.8
• 1.2
• 1.6
• 2
• 2.4
• 2.8

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

V 0.7 0.5 0.3 0.1

• 0.1
• 0.3
• 0.5
• 0.7
• 0.9
• 1.1
• 1.3

Left: horizontal velocity; Right: vertical velocity

Department of Mathematics, University of Pittsburgh 23

Convergence rates for Example 2

mesh u1 − u1,h1,Ω1 rate p − p1,h0,Ω1 rate 4x4 3.54e-01 3.00e-02 8x8 8.60e-02 2.04 7.09e-03 2.08 16x16 2.15e-02 2.00 1.76e-03 2.01 32x32 5.47e-03 1.97 4.44e-04 1.99 64x64 1.40e-03 1.97 1.12e-04 1.99 Table 3: Numerical errors and convergence rates in Ω1. mesh u2 − u2,h0,Ω2 rate p − p2,h0,Ω2 rate 4x4 2.16e-01 1.18e-01 8x8 5.79e-02 1.90 2.87e-02 2.04 16x16 1.47e-02 1.98 7.13e-03 2.01 32x32 3.70e-03 1.99 1.78e-03 2.00 64x64 9.27e-04 2.00 4.45e-04 2.00 Table 4: Numerical errors and convergence rates in Ω2.

Department of Mathematics, University of Pittsburgh 24

Solution of the algebraic flow system

  A B L BT LT     u p λ   =   f1 f2   ⇔ R L LT x λ

• =
• f
• x - subdomain unknowns; λ - interface unknowns

Matrix is symmetric but indefinite. Form the Shur complement system: LTR−1Lλ = LTR−1f (1) Iterative method for (1) (e.g. CG). Each iteration requires evaluating R−1 =   R−1

1

. . . R−1

n

  , i.e., solving subdomain problems. Advantages: subdomain solves in parallel; reuse existing Stokes and Darcy codes

Department of Mathematics, University of Pittsburgh 25

Domain decomposition

Subdomain problems with specified interface normal stress (R−1Lλ): find (u∗

h,i(λ), p∗ h,i(λ)) ∈ Vi h × W i h such that

ai(u∗

h,i(λ), vi) + bi(vi, p∗ h,i(λ))

= −

• ΓI

λ vi · ni, vi ∈ Vi

h,

bi(u∗

h,i(λ), wi)

= 0, wi ∈ W i

h.

λ

Complementary subdomain problems (R−1f): find (¯ uh, ¯ ph) ∈ Vi

h × W i h such

that ai(¯ uh, vi) + bi(vi, ¯ ph) =

• Ωi

fi · vi, vi ∈ Vi

h,

bi(¯ uh, wi) = −

• Ωi

qi wi, wi ∈ W i

h.

BC BC q 1 q2 BC BC BC BC

Interface problem: find λh ∈ Λh such that sh(λh, µ) ≡ −bI(u∗

h(λh), µ) = bI(¯

uh, µ), µ ∈ Λh, Recover global velocity and pressure: uh = u∗

h(λh) + ¯

uh, ph = p∗

h(λh) + ¯

ph.

Department of Mathematics, University of Pittsburgh 26

Subdomain problems

Both types of local problems are well posed due to the local discrete inf-sup conditions. Neumann-Robin BC on ΓI for the Stokes region: −n1 · T · n1 = λ, −n1 · T · τ j − µα

• Kj

u1 · τ j = 0. Dirichlet BC on ΓI for the Darcy region: p2 = λ. Interface (Steklov–Poincar´ e) operator: Sh : Λh → Λh, (Shλ, µ) = sh(λ, µ) ∀ λ, µ ∈ Λh; Sh : λ → [u∗

h(λ) · n]

Interface problem: Find λh ∈ Λh such that Shλh = gh.

Department of Mathematics, University of Pittsburgh 27

Condition number estimate

Lemma: For all λ ∈ M 0

h = {µ ∈ Mh :

• ΓI µ = 0},

C1,1hλ2

ΓI ≤ s1(λ, λ) ≤ C1,2λ2 ΓI

(S1 : H−1/2(ΓI) → H1/2(ΓI)) C2,1λ2

ΓI ≤ s2(λ, λ) ≤ C2,2h−1λ2 ΓI

(S2 : H1/2(ΓI) → H−1/2(ΓI)) Theorem: cond(Sh) = O(h−1) Proof: (C1,1h + C2,1)λ2

ΓI ≤ s(λ, λ) ≤ (C1,2 + C2,2h−1)λ2 Γ,

Department of Mathematics, University of Pittsburgh 28

Algorithm on non-matching grids

Apply the Conjugate Gradient method for Shλh = gh. Computing the action of the operator (needed at each CG iteration):

• Given mortar data λh ∈ Λh, project onto subdomain grids:

λh → Qhiλh

• Solve local Stokes and Darcy problems in parallel with boundary data Qhiλh
• Project fluxes onto the mortar space and compute the jump:

uh,i · ni → QT

hiuh,i · ni,

Shλh = QT

h1uh,1 · n1 + QT h2uh,2 · n2

Department of Mathematics, University of Pittsburgh 29

Multiple subdomains

Three types of interfaces: Stokes–Stokes, Stokes–Darcy, and Darcy–Darcy. The algorithm above works for Stokes–Darcy, and Darcy–Darcy interfaces: only u·n is continuous. On Stokes–Stokes interfaces u · τ also has to be continuous. Augment Sh with the tangential stress as primary variable: Sh :

• n1 · T · n1

n1 · T · τ j

• →
• [u1 · n1]

[u1 · τ j]

• r, equivalently, Sh : n1 · T → [u1],

where [·] is the jump operator. Sh is a Neumann-to-Dirichlet operator. Different number of primary variables on different interfaces.

Department of Mathematics, University of Pittsburgh 30

Condition number for multiple subdomains

ΓS: Stokes-Stokes interfaces ΓD: Stokes-Darcy and Darcy-Darcy interfaces C1(hλ2

ΓS + λ2 ΓD) ≤ s(λ, λ) ≤ C2(λ2 ΓS + h−1λ2 ΓD)

C1hλ2

Γ ≤ s(λ, λ) ≤ C2h−1λ2 Γ

Theorem: cond(Sh) = O(h−2)

Department of Mathematics, University of Pittsburgh 31

Condition number studies: two subdomains

1/h eigmin eigmax iter 5 2.41 4.505 9 10 1.977 7.392 14 20 1.980 14.259 21 40 1.981 28.284 28

Department of Mathematics, University of Pittsburgh 32

Condition number studies: multiple subdomains

4 subdomains 1/h eigmin eigmax iter 4 0.376 12.528 22 8 0.247 24.683 34 16 0.121 50.178 70 32 0.056 102.63 152 64 0.023 209.98 337 Varying number of subdomains 1/H eigmin eigmax iter 1 0.443 1.830 12 2 0.061 4.788 41 4 0.061 12.221 68 8 0.061 47.639 111 16 0.042 75.640 184

Department of Mathematics, University of Pittsburgh 33

Coupling Stokes-Darcy flow with transport

concentration of a chemical c(x, t) porosity of the medium φ(x) diffusion - dispersion tensor Dt(x, t) (symmetric and positive definite) source Q φct + ∇ · (uc − Dt∇c) = φQ , ∀(x, t) ∈ Ω × (0, T) IC: c(x, 0) = c0(x) , ∀x ∈ Ω BC:

• (uc − Dt∇c) · n = ucI · n
• n Γin

(Dt∇c) · n = 0

• n Γout

where Γin := {x ∈ ∂Ω : u · n < 0} and Γout := {x ∈ ∂Ω : u · n ≥ 0}.

Department of Mathematics, University of Pittsburgh 34

Stability and accuracy of a discontinuous Galerkin method

z = −Dt∇c; (ch, zh) : approximation of (c, z); Norm for (ch, zh): |||(ch, zh)|||2 := φ1/2ch(T)2

0,Ω +

T D−1/2

t

zh2

0,Ω dt

+ T

• |uh · n|, (c−

h )2Γ +

• l

|uh · nl|, [ch]2γl

• dt

where [ch] is the jump across the interior boundary edge γl. Theorem: (Vassilev and Y.) |||(ch, zh)||| ≤

• φ1/2c02

0,Ω +

T cI|u · n|1/22

0,Γin dt

1/2 + T φ1/2Q0,Ω dt Theorem: (Vassilev and Y.) |||(c − ch, z − zh)||| ≤ C(hk + u − uhX) ≤ C(hk + hr1 + hr2+1)

Department of Mathematics, University of Pittsburgh 35

Convergence test

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

U 2.8 2.4 2 1.6 1.2 0.8 0.4

• 0.4
• 0.8
• 1.2
• 1.6
• 2
• 2.4
• 2.8

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

V 0.7 0.5 0.3 0.1

• 0.1
• 0.3
• 0.5
• 0.7
• 0.9
• 1.1
• 1.3

D = 10−3 I D = 0 mesh c − chL∞(L2) rate z − zhL2(L2) rate c − chL∞(L2) rate 4x4 1.99e+00 8.95e-03 2.07e+00 8x8 3.27e-01 2.60 2.71e-03 1.72 3.39e-01 2.61 16x16 8.48e-02 1.95 1.20e-03 1.18 9.04e-02 1.91 32x32 2.23e-02 1.93 5.33e-04 1.17 2.59e-02 1.80 64x64 5.60e-03 2.00 1.77e-04 1.59 7.76e-03 1.74

Department of Mathematics, University of Pittsburgh 36

Contaminant transport example

X Y

0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5

permX 1500 1300 1100 900 700 500 300 100

Permeability

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

V

• 0.002
• 0.006
• 0.01
• 0.014
• 0.018
• 0.022
• 0.026
• 0.03
• 0.034
• 0.038
• 0.042
• 0.046
• 0.05

Computed velocity

Department of Mathematics, University of Pittsburgh 37

Transport: initial plume

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

C 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

Time = 0

X

0.2 0.4 0.6 0.8 1

Y

0.2 0.4 0.6 0.8 1

Time = 0

Department of Mathematics, University of Pittsburgh 38

Transport: initial plume

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

C 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

Time = 3

X

0.2 0.4 0.6 0.8 1

Y

0.2 0.4 0.6 0.8 1

Time = 3

Department of Mathematics, University of Pittsburgh 39

Transport: initial plume

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

C 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

Time = 5

X

0.2 0.4 0.6 0.8 1

Y

0.2 0.4 0.6 0.8 1

Time = 5

Department of Mathematics, University of Pittsburgh 40

Transport: initial plume

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

C 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

Time = 9

X

0.2 0.4 0.6 0.8 1

Y

0.2 0.4 0.6 0.8 1

Time = 9

Department of Mathematics, University of Pittsburgh 41

Transport: initial plume

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

C 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

Time = 16

X

0.2 0.4 0.6 0.8 1

Y

0.2 0.4 0.6 0.8 1

Time = 16

Department of Mathematics, University of Pittsburgh 42

Transport: inflow

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

C 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

Time = 2

X

0.2 0.4 0.6 0.8 1

Y

0.2 0.4 0.6 0.8 1

Time = 2

Department of Mathematics, University of Pittsburgh 43

Transport: inflow

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

C 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

Time = 11

X

0.2 0.4 0.6 0.8 1

Y

0.2 0.4 0.6 0.8 1

Time = 11

Department of Mathematics, University of Pittsburgh 44

Transport: inflow

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

C 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

Time = 17

X

0.2 0.4 0.6 0.8 1

Y

0.2 0.4 0.6 0.8 1

Time = 17

Department of Mathematics, University of Pittsburgh 45

Summary

• Well-posed mathematical flow model based on continuity of flux and normal

stress, and the Beavers-Joseph-Saffman condition

• Error analysis for mortar multiscale finite element approximations with very

general mortar conditions

• Non-overlapping domain decomposition: s.p.d. interface problem; local Stokes

and Darcy solves required at every CG iteration

• Condition number of interface problem: O(h−1) for two subdomains; O(h−2) for

multiple subdomains

• Conforming and DG Stokes elements and mixed finite elements for flow are

coupled with LDG for transport

Department of Mathematics, University of Pittsburgh 46

Current and future work

• Balancing preconditioner for the interface problem
• A posteriori error estimation and adaptive mesh refinement
• Extensions to Navier-Stokes and Brinkman
• Polygonal elements using MFD and DG

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8

K 100 50 20 10 5 2 1

Permeability

X Y

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

P 0.19 0.18 0.17 0.16 0.15 0.14 0.13 0.12 0.11 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01

Pressure and velocity

Department of Mathematics, University of Pittsburgh 47