Graziani F. (editor) Computational Methods in Transport

Подождите немного. Документ загружается.

Transport Approximations in Partially Diﬀusive Media 377

Deﬁne the harmonic decomposition

Fu(n)=ˆu

2π



2π

−inθ

u(ω)dθ(ω), F

−1

ˆu(θ)=u(θ)=



n∈Z

inθ

ˆu

(6)

In this basis, the operator G may be written as a (diagonal) Fourier multiplier:

G = F

−1

(σ − k

)F , (7)

where the coeﬃcients σ and k

are such that |k

| <k

<σfor |n|≥1.

Also, since G is a real-valued operator, we verify that k

−n

,wherethe

upper bar denotes complex conjugation. Because G is symmetric, we have

here k

−n

= k

real-valued. The above decomposition means that Ge

inθ

(σ − k

inθ

and implies that

−1

= F

−1

(σ − k

)

−1

F . (8)

We also verify the following convenient expression for the diﬀerential operator

ω ·∇= e

iθ

∂ + e

−iθ

∂, ∂ =



∂

∂x

− i

∂

∂y



∂ =



∂

∂x

+ i

∂

∂y



. (9)

This implies that



(ω ·∇v)

= ∂ˆv

n−1

∂ˆv

n+1

. We ﬁnally identify x =(x, y)

with z = x + iy in the complex plane.

With this notation, we verify that the transport equation is equivalent to



iθ

∂ + e

−iθ

∂



u + Gu = q (10)

with appropriate boundary conditions. In the Fourier domain this is nothing

but

∂ˆu

n−1

∂ˆu

n+1

+(σ − k

)ˆu

=ˆq

,Ω× Z .

(11)

Recall the Parseval relation



uvdµ =



n∈Z

ˆu

ˆv



n∈Z

ˆu

ˆv

−n

since u

and v are real-valued functions. The variational formulation (4) still holds

with a(u, v)andL(v) recast as:

a(u, v)=



Ω



n∈Z



a − k

(∂ˆu

n−1

∂ˆu

n+1

)(∂ˆv

−n−1

∂ˆv

−n+1

)

+(σ − k

)ˆu

ˆv

−n



dµ(z)+



uvω · ν dq ,

L(v)=



Ω



n∈Z



a − k

ˆq

(∂ˆv

−n−1

∂ˆv

−n+1

)+ˆq

ˆv

−n



dµ(z)



−

gv|ω · ν| dq .

(12)

The above formulae prove very useful in the analysis of the transport solution

in the diﬀusive regime.

378 G. Bal

2.2 Variational Formulation in the Diﬀusive Regime

We now consider the regime of high collisions and small absorption. This

regime is characterized by replacing G, q and g by [12]

u(x, ω)=



n−1

k(x, ω



− ω)



u(x, ω) − u(x, ω



)



dµ(ω



)+εσ

= εq, g

= εg .

(13)

Here g

is of order ε to avoid the presence of boundary layers at the leading

order [12] (we will see below that a term of order ε

1/2

rather than ε would have

been suﬃcient). We recast G

= ε

−1

Q + εσ

,whereσ

is uniformly bounded

from below by a positive constant so that G

−1

is bounded on L

(X) though

not uniformly in ε. The local ﬁrst-order transport equation (1) reads in this

regime:

ω ·∇u

+ G

= εq . (14)

We recast (4) in this regime as ﬁnding u ∈ W such that

(u, v)=L

(v), ∀v ∈ W, (15)

where

(u, v)=



((εG

)

−1

(ω ·∇u)ω ·∇v + ε

−1

(u)v)dp +



uvω · ν dq ,

(v)=





−1

(q)ω ·∇v + qv



dp +



−

gv|ω · ν| dq .

(16)

We now derive some properties of the above variational formulation and

the above bilinear form that make them attractive theoretically and compu-

tationally. We ﬁrst assume that G

is invertible (in L

(V )) with inverse given

−1

u =

εσ

¯u + εH

u, H

u =0, (17)

where H

u is a symmetric and bounded operator in L

(X) with norm

bounded by α

−1

< ∞,and¯u =



n−1

u(ω)dµ(ω) is the angular average

of u. This property holds when G

is decomposed over spherical harmonics;

see for instance (37) and (41) below. This allows us to recast (16) as

(u, v)=





ω ·∇u ω ·∇v + H

(ω ·∇u)ω ·∇v +





uvω · ν dq (18)

(v)=





εσ

ω ·∇v + εH

(ω ·∇v)+v



qdp+



−

gvω · ν dq .

Transport Approximations in Partially Diﬀusive Media 379

Let us assume that Q and H

arecoerciveonL

(X) with coercivity

constants α and σ

−1

, respectively, in the sense that:

(Qu, u) ≥ αu

, (H

u, u) ≥

u

. (19)

Here f is the L

(X)normoff . These are reasonable assumptions when

can be diagonalized in the basis of spherical harmonics; see (37) and (41)

below. Let us also assume that the absorption 0 <σ

≤ σ

(x) ≤ σ

on Ω.

We then verify that

ω ·∇u



ω ·∇u



u

− u



+ σ

u



u

ε|Γ



≤ a

) . (20)

Here f 

denotes the L

(Γ

; |ω ·ν|dq)normoff deﬁned on Γ

. This implies

that a

is coercive on W with constant of coercivity independent of ε.The

constant of coercivity also depends on the scattering kernel only through the

constants α and σ

−1

It is also straightforward to observe that a

(u, v) is continuous on W ×W ,

though with a continuity constant that depends on ε. In order to obtain a

constant of continuity independent of ε, the Hilbert space may be equipped

with a norm ·

deﬁned as the (square root of the) left-hand side of (20);

this is in essence the norm introduced in [24]. On W

(i.e., W equipped with

·

) the form a

is coercive and continuous with constants of coercivity

and continuity independent of ε (see also Sect. 2.7). The source term is also

bounded on W

with a constant independent of ε, and bounded on W with

a constant that depends on ε:

)| 

ω ·∇u

q +

εσ

¯qω ·∇u

 + qu

 + u

ε|Γ



g

(21)

The notation a  b stands for a ≤ Cb for some constant C independent of ε.

We deduce from the constraints on a

and L

that the unique solution u

the variational formulation (15) satisﬁes the (a priori) estimate

ω ·∇u

 +



ω ·∇u

 + u

  q+

√

εg

. (22)

This implies that

)|  q

+ εg

, (23)

so that from (20),

u

− ¯u



(X)



√

α(εq + ε

3/2

g

) . (24)

To summarize we have the following result:

380 G. Bal

Lemma 1. Provided that (19) holds, the unique solution u

of (15) veriﬁes

that

ω·∇u

+u

+



ω ·∇u

+

u

−¯u



(X)

√

u

ε|Γ



 q+

√

εg

(25)

provided that q ∈ L

(X) and g ∈ L

(Γ

−

; |ω · ν|dq). We recognize the above

left-hand side as u



Least-Square Formulation

Let us introduce the operator

=(εG

)

−1/2

(ω ·∇+ G

) . (26)

Using the divergence theorem on ∇·(uvω)=(ω ·∇u)v + uω ·∇v,weverify

that

(u, v)=(L

u, L

v)+

u, v , (27)

where (·, ·) is the usual inner product on L

(X)and·, · is the inner product

on L

(Γ

−

, |ω · ν|dq). We also have

(v)=(q

, L

v)+g, v,q

= ε

1/2

−1/2

(q) . (28)

Since G

, whence a

, is symmetric, the variational formulation (15) is equiv-

alent to minimizing the functional:

arg min

v∈W

(v, v) − L

(v) , (29)

which is also equivalent to minimizing the following least-square problem

arg min

u∈W

L

u − q



(X)

2ε

u − εg

(Γ

−

;|ω·ν|dq)

, (30)

as can easily be veriﬁed.

We have thus recast solving the transport solution as minimizing the least-

square problem associated to a variational form a

that is W -coercive with

constant of coercivity independent of ε. The boundary conditions are also

accounted for in a variational sense as the functions in the space W need not

satisfy the boundary conditions exactly. These very important properties are

a central element in the numerical methods developed in [24, 25], which are

based on Galerkin projections; see also Sect. 2.7. Note that the scaling oper-

ator S used in the above references is replaced in our analysis by (εG

)

−1/2

The method developed in this paper may thus be seen as a case of ﬁrst-order

system least-squares (FOSLS) [24, 25].

Because of (20), Galerkin methods, which are orthogonal projections

onto subspaces of W with respect to the bilinear form a

, provide lower-

dimensional approximations that are expected to be valid both in the trans-

port regime (ε ∼ 1) and the diﬀusive regime (ε  1) as in [24, 25]. We now

Transport Approximations in Partially Diﬀusive Media 381

consider two Galerkin methods, which are not discretizations of the transport

solution, but rather projections onto smaller subsets of W that are physi-

cally relevant in the diﬀusive regime. The ﬁrst method consists of projecting

the transport solution onto functions that only depend on space and not on

the angular variable. The resulting diﬀusion approximation is analyzed in the

rest of this section. In the next section we project the transport solution onto

functions that depend on the spatial variable only in parts of the domain.

This allows us to couple the diﬀusion approximation to the full transport

solution in regions where diﬀusion may not be valid.

2.3 Diﬀusion by Orthogonal Projection

We want to use the above variational formulation to deduce the limit of

the transport solution in the diﬀusion limit. Diﬀusion is characterized by

high scattering and small absorption. High scattering implies that the initial

directional content of the particles is quickly lost through interactions. It is

therefore reasonable to assume that u(x, ω) does not depend on ω in a ﬁrst

approximation. Such a condition is easy to implement in a variational setting:

we orthogonally (with respect to a

(·, ·)) project the solution u of (4) onto

functions that depend only on x. Let us deﬁne



f ∈ W, f = f(x)



≡ H

(Ω),

⊥

Dε



f ∈ W, a

(f,v)=0∀v ∈ W



. (31)

We verify that W = W

⊕ W

⊥

Dε

since a

(·, ·) is an inner product on W [5].

Then the orthogonal projection Π

of W to W

for the inner product a

(u, u)

allows us to deﬁne the “diﬀusion” solution

= Π

, (32)

where u

is the solution to (15). By deﬁnition, we have that U

is the solution

in W

,v)=L

(v),v∈ W

. (33)

We may recast the above equation as



Ω

∇U ·∇V + σ

UV)dx +



∂Ω

UVdσ(x)



Ω





n−1

−1

(q)ωdµ(ω)



·∇V + qV



dx +



−

gV |ω · ν| dq ,

for all V ∈ W

,where



n−1

(εG

)

−1

(ω) ⊗ωdµ(ω)



n−1

|ω · ν|dµ(ω) .

(34)

382 G. Bal

We verify that c

and that c

. After classical integrations by parts,

this means that U

is the solution of the following diﬀusion equation:

−∇ · D

∇U

+ σ

= q −∇·





n−1

−1

(q)ωdµ(ω)



ν · D

∇U

= J(x) ≡



ω·ν<0

g(x, ω)|ω · ν|dµ(ω) .

(35)

Here we assume that q vanishes at the domain boundary ∂Ω to simplify. This

is consistent with our choice of incoming conditions of size εg as O(1) source

terms at the domain boundary result in boundary layer analyses that we do

not want to dwell into here; see [12].

It remains to understand whether the solution U

, which is obviously

uniquely deﬁned, is uniformly bounded in H

(Ω) and to obtain an error

estimate for u

− U

and for similar diﬀusion approximations of u

Let us ﬁrst consider the simpler case where scattering is isotropic, which

implies that

u =

(u − ¯u)+εσ

u, (36)

where ¯u =



n−1

u(ω)dµ(ω)andwherek

(x)=k(x, ω −ω



). In this context

we verify that

−1

u =

¯u +

u, σ

= k

+ ε

. (37)

Assuming that q =¯q and that g ≡ 0 to simplify, this implies that the trans-

port equation takes the form





ω ·∇u +

ω ·∇u



ω ·∇v +

(u − ¯u)v + σ





uvω · ν dq =





ω ·∇v

εσ

+ v



dp .

(38)

Upon choosing u and v in W

in the above equation, we ﬁnd the variational

formulation for U



Ω



nσ

∇U ·∇V + σ



dx +



∂Ω

UVdσ =



Ω

qV dx . (39)

This is nothing but the variational formulation of (35) with the diﬀusion

coeﬃcient given by D

=(nσ

)

−1

as usual and the right-hand side given by

q(x). We thus obtain the classical diﬀusion equation [19].

2.4 Harmonic Decomposition and Diﬀusion Approximation

The generalization of the derivation of the above diﬀusion solution to arbi-

trary scattering kernels is relatively straightforward as long as G

−1

can be

Transport Approximations in Partially Diﬀusive Media 383

explicitly characterized. We restrict ourselves here to the two-dimensional

case with scattering kernel of convolution type and use explicitly the vari-

ational form written in the harmonic basis. Within this context and in the

diﬀusive regime, G

is given by

= F

−1



− k

+ εσ



F . (40)

We recall that the Fourier transform is considered in the angular variable

only. The inverse of G

is then

−1

= F

−1

− k

+ ε

F≡F

−1

− k

F (41)

where σ

= k

+ ε

, This yields (37) when k

=0for|n|≥1 and (17)

in general, where H

u is a bounded operator in L

(X). Let us introduce the

spectral gap of Q:

α =min

n≥1,x∈Ω

(x) − k

(x)) > 0 , (42)

and the minimum of σ

=min

x∈Ω

(x) > 0 . (43)

We then verify that the coercivity constraints (19) are veriﬁed.

With this explicit form for G

, the terms of the variational equation (15)

now take the form

(u, v)=



Ω



n∈Z



− k

+ε



(ω ·∇u)



(ω ·∇v)

−n



− k

+ σ



ˆu

ˆv

−n



dx + ε

−1



uvω · ν dq,

(v)=



Ω



n∈Z



εˆq



(ω ·∇v)

−n

− k

+ ε

+ˆq

ˆv

−n



dx +



−

gv|ω · ν| dq . (44)

The generalization of (39) to more general scattering terms may easily be

carried out thanks to the above formulae. Indeed, U and V belong to W

and only if all their Fourier coeﬃcients

and

vanish for |n|≥1. As a

result, U is the unique solution to the following variational problem



Ω



∂

− k

+ ε

+ σ



dx +



∂Ω

UVdσ



Ω



ˆq

∂

+ˆq

−1

∂

− k

+ ε

+ˆq





∂Ω





ω·ν<0

g(x, ω)|ω · ν|dµ(ω)



dσ . (45)

384 G. Bal

We have used here the symmetry of G implying that k

−1

= k

. We denote

by ˆq

−1

+ iq

) and identify it with the vector q =(q

) in Cartesian

coordinates. The above variational formulation shows that U is the (weak)

solution in W

of the following diﬀusion equation:

−∇ · D

∇U + σ

U = q − ε∇·D

q,Ω

∂U

∂n

U = J, ∂Ω .

(46)

The diﬀusion coeﬃcient is given by

2(σ

− k

)

2(k

− k

+ ε

)

. (47)

This is the usual expression for the diﬀusion coeﬃcient.

Remark 1. The property that the constant of coercivity of a

is independent

of the mean free path ε is very important to obtain the diﬀusion solution by

orthogonal projection. Consider for instance the transport problem

≡

ω ·∇u

= q(x)inX, u

=0 onΓ

−

, (48)

and the variational formulation: ﬁnd u

such that

˜a

,v)=(w

v)=(q,w

v), ∀v ∈ W. (49)

Here w

(x) is a weight function that could take the value 1 as in the introduc-

tion in [24] or σ

−1

as in [7]. The above problem is equivalent to (48) for w

uniformly bounded from above and below by positive constants as is shown

in [7]. This results from the fact that ˜a

,v)iscoerciveonW . However the

constant of coercivity is not independent of ε. The equivalence between (48)

and (49) thus somewhat degrades as ε → 0. The consequence is that (49)

becomes inaccurate in the diﬀusive regime with dire consequences when it

comes to discretizations as pointed out in [24]. The same consequence arises

when the natural orthogonal projection

of u

onto W

is considered. As-

suming that G

is isotropic as in (36) to simplify, we obtain that

solves the equation

−∇ ·

∇

+ w

= w

q, in Ω. (50)

Here n =2, 3 is the spatial dimension. Although the local equilibrium

≈ σ

−1

q is veriﬁed in the limit of very strong absorption and sources,

the diﬀusion tensor is not correct, independently of the choice of the weight

. This indicates that variational formulations of the form (49) should not

be used in the diﬀusive regime and should be replaced by (4) or by rescaled

formulations as in [24].

Transport Approximations in Partially Diﬀusive Media 385

2.5 Convergence Result and Error Estimates

Let u

(x, ω) be the solution of the transport equation (15) and U

= Π

the diﬀusion approximation solution of (33). Deﬁne δu

= u

− U

.Weﬁrst

show that δu

convergesto0inW as ε → 0. The proof is similar to that

in [5].

Thanks to the orthogonal projection, we ﬁrst obtain that δu

∈ W

⊥

Dε

Consequently, we have

(δu

,δu

)=a

,δu

)=L

(δu

) .

We then deduce from (20) and (21) that

δu



δu

− δu

  q+

√

εg

Since δu

is uniformly bounded in W , which is a Hilbert space with a unit ball

compact for the weak topology, we deduce the existence of a subsequence of

δu

converging to w. However the subsequence of δu

converges to the same

limit so that w =

w ∈ W

. This implies that

(δu

,w)=0.

Upon passing to the limit in the above expression, we ﬁnd that



Ω



− k





(ω ·∇w)



+ σ

|w|



dx =0

and that w

|∂Ω

= 0. This implies that w = 0, whence that δu

converges to

0. Note that all we have used to get this convergence result is that q and

g

are bounded (in the L

−sense).

In order to obtain error estimates for δu

, additional regularity of the

solutions and source terms is required. Let us recall (18):

(u, v)=





ω ·∇u ω ·∇v + H

(ω ·∇u)ω ·∇v +





uvω · ν dq ,

(v)=





q − ω ·∇

εσ

− εω ·∇H

(q)



vdp+



−

gvω · ν dq ,

(51)

assuming that q vanishes on ∂Ω to simplify. We split the source term as

= L

+ L

(v)=



qv dp −



εω ·∇H

(q)vdp+



−

gvω · ν dq ,

(v)=



−ω ·∇

εσ

vdp.

(52)

386 G. Bal

The term appearing in L

(v)isoforderε

−1

and cannot be estimated from

the variational formulation. In order to estimate it explicitly, we deﬁne W

the subspace of W of functions u(x, ω) such that ˆu

(x)=0for|n| =1.This

is thus the Hilbert space of functions linear in ω.LetΠ

1ε

be the orthogonal

projection of W onto W

with respect to a

(·, ·). We deﬁne U

1ε

= Π

1ε

, i.e.,

1ε

)=L

), ∀V

∈ W

. (53)

Thanks to the two terms of order ε

−2

in a

in (51), we verify that

U

1ε

 + ω ·∇U

1ε

  εq

+ ε

3/2

g

. (54)

In what follows, we need a similar estimate for the following solution:

(

1ε

)=L

), ∀V

∈ W

. (55)

We verify as above that



1ε

 + ω ·∇

1ε

  ε∇q+ ε

3/2

g

. (56)

In order to estimate the error coming from the source term L

(v), we

deﬁne ψ

2ε

((εG

)ψ

2ε

,v)=L

(v) − a

,v), ∀v ∈ W

2ε

=0.

(57)

Here (·, ·) is the usual inner product of L

(X). We verify that ψ

2ε

is well-posed

and is given by

2ε

(x, ω)=H



ω ·∇H

(ω) ·∇U

−∇·D

∇U

+ εω ·∇H

(q) − 2D

ε∇·q



. (58)

The term within brackets in the above expression is mean zero by construction

so that all equalities in (57) are veriﬁed. Let us now decompose the exact

solution as

= u

0ε

+ u

1ε

kε

,v)=L

(v), ∀v ∈ W, k =0, 1 . (59)

We introduce the following decompositions:

0ε

= U

+ ε

2ε

+ δ

0ε

1ε

+ δ

1ε

, (60)

and estimate both terms δ

kε

for k =0, 1. We start with k = 0 and obtain

from the transport equation a

0ε

,δ

0ε

) − L

(δ

0ε

)=0that