有限体積法による差分方程式の導出

変換方程式(4.16)と式(4.17)の離散化には､数値的な安定性を考慮して､有限体積法を用い た差分方程式を導出した。これは､地下構造モデルには、隣接する節点で導電率が激しく変 化する場合もあることから､節点での導電率およびその空間的微分の直接的な定義を避け､

連続性を有するfluxを扱う有限体積法を用いた方が良いと考えられるからである。

いま､ある1つの節点pに注目し､その節点pを中心とした

∆Γ

_pを境界面とする微小体積

Vp

∆

について､式(4.16)の両辺で体積積分を行う。ここで、Gaussの定理を用いれば、体積積 分を面積積分に置き換えることができ、次式が得られる。

 

V

Px

∆ ∆ V

^P_ξ

1 Px

∆Γ+

1 Px

∆Γ− 2 Px

∆Γ− ₂

Px

∆Γ+

3 Px

∆Γ+ 3 Px

∆Γ− ₃

Pξ−

∆Γ

3Pξ+

∆Γ

2 Pξ+

∆Γ

2 Pξ−

∆Γ 1

P−ξ

∆Γ

1 Pξ+

∆Γ

1 2 3

ξ ξ ξ

1 2 3

x x x

Physical Space Intermediate Space 

Fig.4.3 Illustration of volume element and its out surface in phys ical and intermediate space: Coordinate transformation to the intermediate space facilitates the surface integration in the finite volume method, because the boundary surface of

∆

V_pis piecewise coincident with certain coordinate surface elements in the intermediate space.

( ) 0

x x

p p

x x i x i

V

s dv s d

σ σ

∆ ∆Γ

∇ ⋅ ∇ = ∇ ⋅ =

∫ Ñ ∫ ^Γ

^(4.25)

一般に､物理空間および計算空間では､境界面は座標面と平行関係にはなく、複雑な曲面 を形成するため､その境界面積分は容易ではない。しかし､媒介空間では､その定義から､

∆Γ

p^の6枚の側面がそれぞれ座標平面の面要素と重なる(Fig.4.3)ことから、計算は容易とな る。

次の恒等式、

3

1 i m x i

i m

s s

_ξ

=

ξ

∇ = ∂

∑ ∂ ^a

および

3 3

1 1

( , , 1,2,3)

j

k k l l k l

j j

d d d g d d

j k l circulatingwith

ξ ξ

ξ ξ

ξ ξ ξ ξ

= =

= × =

⇒

∑ ∑

Γ a a a

を式(4.25)に代入すると、媒介空間で対応する境界面、

( )

3

1

j j

p p p

j

ξ − +

=

∆Γ = ∑ ∆Γ +∆Γ

での面積積分へと変換でき、次式が得られる。

3 3 3

1 1 1

( )

0

( 1,2,3; , , 1,2,3)

p

j j

p p

i V

jm i jm i

k l k l

j m m m m

s dv

s s

g g d d g g d d

i j k l circulatingwith

ξ ξ ξ ξ

σ

σ ξ ξ σ ξ ξ

ξ ξ

+ −

∆

= ∆Γ = ∆Γ =

∇⋅ ∇ =

  ∂   ∂  

   −    =

  ∂   ∂  

 

= ⇒

∫

∑ ∫ ∑ ∫ ∑

^(4.26)

式(4.26)の導出過程では、基本計量テンソル反変成分の定義式、

jm j m

x j x m

g

_ξ

= a

_ξ

⋅ a

_ξ

= ∇ ξ ⋅ ∇ ξ

を用いている。

同じような手順で、楕円型変換方程式(4.17)を微小体積

∆ V

_p^xについて体積積分すると、次 式が得られる。

3 3 3

1 1 1

( )

0

( 1,2,3; , , 1,2,3)

p

j j

p p

i V

jm i jm i

k l k l

j m m m m

x dv

x x

g g d d g g d d

i j k l circulating with

ξ ξ

ξ ξ

ξ ξ

ξ ξ

ξ ξ

+ −

∆

= ∆Γ = ∆Γ =

∇⋅ ∇ =

  ∂   ∂  

   −    =

  ∂   ∂  

 

= ⇒

∫

∑ ∫ ∑ ∫ ∑

^(4.27)

さらに、式(4.4)と式(4.5)から、

2

jm jm

g

_ξ

= σ G

_ξ

3

g_ξ 1 G_ξ

=

σ

がわかるので、この2つの式を式(4.27)に代入すると、次式が得られる。

3 3 3

1 1 1

( )

1 1

0

( 1,2,3; , , 1,2,3)

p

j j

p p

i V

jm i jm i

k l k l

j m m m m

x dv

x x

G G d d G G d d

i j k l circulatingwith

ξ ξ

ξ ξ

ξ ξ

ξ ξ

σ ξ σ ξ

+ −

∆

= ∆Γ = ∆Γ =

∇⋅ ∇ =

  ∂   ∂  

   −    =

  ∂   ∂  

 

= ⇒

∫

∑ ∫ ∑ ∫ ∑

^(4.28)

ただし､

∆Γ

^j_p⁺

, ∆Γ

^j_p^{−は、それぞれ境界面}

∆Γ

p^のξj正方向と負方向の座標面要素を示す。ま た､j k l, , は､jが1のときk,l は2,3､jが2のときk,lは3,1､jが3のときk,l は1,2である。

jm

x j x m

g

_ξ

= ∇ ξ g ∇ ξ

、

J = g

_ξ

= ( g

^ξ_jm

)

は､それぞれ物理空間に対する媒介空間の基本計 量テンソルの反変成分、変換ヤコビアンであり、

G

_ξ^jm

= ∇

_s

ξ

_j

g ∇

_{s m}

ξ

､

J = G

_ξ

= ( G

^ξ_jm

)

^は､

それぞれ計算空間に対する媒介空間の基本計量テンソルの反変成分、変換ヤコビアンであ る。

式(4.26)と式(4.28)に、それぞれ

0

i i i

m m m

s s s

ξ ξ ξ

∂ = ∂ + ∂∆

∂ ∂ ∂

^と

0

i i i

m m m

x x x

ξ ξ ξ

∂ = ∂ + ∂∆

∂ ∂ ∂

^{を代入すると、}

座標修正量を未知数とする、ピカード反復計算法への適用可能な変換方程式が得られる。

3 3 3

1 1 1

( 1,2,3; , , 1,2,3)

j j

p p

jm i jm i

k l k l p

j m m m m

s s

g g d d g g d d Q

i j k l circulating with

ξ ξ ξ ξ

σ ξ ξ σ ξ ξ

ξ ξ

+ −

= ∆Γ = ∆Γ =

  ∂∆   ∂∆  

   −    =

  ∂   ∂  

 

= ⇒

∑ ∫ ∑ ∫ ∑

  (4.29)

3 3 3

1 1 1

1 1

( 1,2,3; , , 1,2,3)

j j

p p

jm i jm i

k l k l p

j m m m m

x x

G G d d G G d d R

i j k l circulating with

ξ ξ ξ ξ ξ ξ ξ ξ

σ ξ σ ξ

+ −

= ∆Γ = ∆Γ =

  ∂∆   ∂∆  

   −    =

  ∂   ∂  

 

= ⇒

∑ ∫ ∑ ∫ ∑

  (4.30)

ただし、

0 0

3 3 3

1 j 1 j 1

p p

jm i jm i

p k l k l

j m m m m

s s

Q σ g_ξ g_ξ dξ ξd σ g_ξ g_ξ d dξ ξ

ξ ξ

+ −

= ∆Γ = ∆Γ =

  ∂   ∂  

 

= − ∑   ∫   ∑ ∂   − ∫   ∑ ∂    

0 0

3 3 3

1 1 1

1 1

j j

p p

jm i jm i

p k l k l

j m m m m

x x

R G_ξ G_ξ d dξ ξ G_ξ G_ξ d dξ ξ

σ ξ σ ξ

+ −

= ∆Γ = ∆Γ =

  ∂   ∂  

 

= − ∑   ∫   ∑ ∂   − ∫   ∑ ∂    

さらに、式(4.29)の

∂∆ ∂ s

_i

ξ

_m^{を中心差分近似し、}σ g g_{ξ ξ}^jmを各面要素での平均値で近似 し、整理すると、未知変数

s

_iの修正量に関する次の19点差分方程式が得られる。

19

1 n

n i p

n p

C s Q

=

 ∆  =

 

 ∑ 

^(4.31)

ただし、

23 32

1

(

3 3 3 5 5 5

) 4

C = σ g g + σ g g

13 31

2

(

1 1 1 5 5 5

) 4

C = σ g g + σ g g

33 13 13 23 23

3 5 5 5

(

1 1 1 2 2 2 3 3 3 4 4 4

) 4

C = σ g g + σ g g − σ g g + σ g g − σ g g

13 31

4

(

2 2 2 5 5 5

) 4

C = − σ g g + σ g g

32 23

5

(

5 5 5 4 4 4

) 4

C = − σ g g + σ g g

12 21

6

(

1 1 1 3 3 3

) 4

C = σ g g + σ g g

22 12 12 32 32

7 3 3 3

(

1 1 1 2 2 2 5 5 5 6 6 6

) 4

C = σ g g + σ g g − σ g g + σ g g − σ g g

12 21

8

(

2 2 2 3 3 3

) 4

C = − σ g g + σ g g

11 21 21 31 31

9 1 1 1

(

3 3 3 4 4 4 5 5 5 6 6 6

) 4

C = σ g g + σ g g − σ g g + σ g g − σ g g

11 11 22 22 33 33

10

(

1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6

)

C = − σ g g + σ g g + σ g g + σ g g + σ g g + σ g g

 

11 21 21 31 31

11 2 2 2

(

4 4 4 3 3 3 6 6 6 5 5 5

) 4

C = σ g g + σ g g − σ g g + σ g g − σ g g

12 21

12

(

1 1 1 4 4 4

) 4

C = − σ g g + σ g g

22 12 12 32 32

13 4 4 4

(

2 2 2 1 1 1 6 6 6 5 5 5

) 4

C = σ g g + σ g g − σ g g + σ g g − σ g g

12 21

14

(

2 2 2 4 4 4

) 4

C = σ g g + σ g g

23 32

15

(

3 3 3 6 6 6

) 4

C = − σ g g + σ g g

13 31

16

(

1 1 1 6 6 6

) 4

C = − σ g g + σ g g

33 13 13 23 23

17 6 6 6

(

2 2 2 1 1 1 4 4 4 4 4 4

) 4

C = σ g g + σ g g − σ g g + σ g g − σ g g

13 31

18

(

2 2 2 6 6 6

) 4

C = σ g g + σ g g

23 32

19

(

4 4 4 6 6 6

) 4

C = σ g g + σ g g

ここで、

σ

_k

g g

_{k k}^jm

( k = 1 : 6)

は、k 番目の面要素でσ g g_{ξ ξ}^{jmの平均値を意味する。}

( ) ( 1)

n n n

i i i

s s it s it

∆ = − −

^{であり､このとき}s it_iⁿ( )^は､節点p とかかわる19 の節点の中､n 番節 点の計算空間座標変数

s

_iに関する第it回反復計算結果である。Fig.4.4に､相関する各節点の

位置関係を示す。

式(4.30) において､偏導関数

∂ ∂ x

_i

ξ

_mを中心差分で近似すると､未知変数

x

_iの修正量に関 する19点差分方程式が得られる。

1

s

i

∆

2

s

i

∆ ∆ s

_i³

∆ s

ⁱ⁴

5

s

i

∆

6

s

i

∆ ∆ s

_i⁷

8

s

i

∆

9

s

i

∆ ∆ s

¹⁰i

11

s

i

∆

12

s

i

∆ ∆ s

¹³i

14

s

i

∆

15

s

i

∆

16

s

i

∆ ∆ s

i¹⁷

∆ s

¹⁸ⁱ

19

s

i

∆

Fig.4.4 Illustration of the finite difference scheme: Central difference approximation of equation (4.29) leads to a 19-point finite difference scheme.

For equation (4.30), a similar scheme also obtained

19

1 n

n i p

n p

D x R

=

 ∆  =

 

 ∑ 

^(4.32)

ただし、

23 32

1

(

3 3 3 5 5 5

) 4

D = ρ G G + ρ G G

13 31

2

(

1 1 1 5 5 5

) 4

D = ρ G G + ρ G G

33 13 13 23 23

3 5 5 5

(

1 1 1 2 2 2 3 3 3 4 4 4

) 4

D = ρ G G + ρ G G − ρ G G + ρ G G − ρ G G

13 31

4

(

2 2 2 5 5 5

) 4

D = − ρ G G + ρ G G

32 23

5

(

5 5 5 4 4 4

) 4

D = − ρ G G + ρ G G

12 21

6

(

1 1 1 3 3 3

) 4

D = ρ G G + ρ G G

22 12 12 32 32

7 3 3 3

(

1 1 1 2 2 2 5 5 5 6 6 6

) 4

D = ρ G G + ρ G G − ρ G G + ρ G G − ρ G G

12 21

8

(

2 2 2 3 3 3

) 4

D = − ρ G G + ρ G G

11 21 21 31 31

9 1 1 1

(

3 3 3 4 4 4 5 5 5 6 6 6

) 4

D = ρ G G + ρ G G − ρ G G + ρ G G − ρ G G

11 11 22 22 33 33

10

(

1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6

)

D = − ρ G G + ρ G G + ρ G G + ρ G G + ρ G G + ρ G G

 

11 21 21 31 31

11 2 2 2

(

4 4 4 3 3 3 6 6 6 5 5 5

) 4

D = ρ G G + ρ G G − ρ G G + ρ G G − ρ G G

12 21

12

(

1 1 1 4 4 4

) 4

D = − ρ G G + ρ G G

22 12 12 32 32

13 4 4 4

(

2 2 2 1 1 1 6 6 6 5 5 5

) 4

D = ρ G G + ρ G G − ρ G G + ρ G G − ρ G G

12 21

14

(

2 2 2 4 4 4

) 4

D = ρ G G + ρ G G

23 32

15

(

3 3 3 6 6 6

) 4

D = − ρ G G + ρ G G

13 31

16

(

1 1 1 6 6 6

) 4

D = − ρ G G + ρ G G

33 13 13 23 23

17 6 6 6

(

2 2 2 1 1 1 4 4 4 4 4 4

) 4

D = ρ G G + ρ G G − ρ G G + ρ G G − ρ G G

13 31

18

(

2 2 2 6 6 6

) 4

D = ρ G G + ρ G G

23 32

19

(

4 4 4 6 6 6

) 4

D = ρ G G + ρ G G

ここで、

ρ

_k

G G

_{k k}^jm

( k = 1 : 6)

は、k 番目の面要素でρ G G_{ξ ξ}^{jmの平均値を意味する。}

( ) ( 1)

n n n

i i i

x x it x it

∆ = − −

^{であり､このとき}x it_iⁿ( )^は､節点pとかかわる19の節点の中､n番節

点の物理空間座標変数

x

_iに関する第it回反復計算結果である。

式(4.31)(あるいは式(4.32))は、グリッド系のすべての内点に対して適用することができる。

境界点の場合には､境界条件式(4.18)と式(4.19) (あるいは式(4.20)と式(4.21))を用いて､境界 外の仮想点での未知量を、境界点あるいは内点での未知量に置き換えればよい。これより､

グリッド系のすべての節点について、それぞれ一つの差分方程式が与えられることから､構 成される全体方程式を解くことにより､計算空間(あるいは物理空間)グリッド系の各節点座 標の修正量が求められる。

ドキュメント内 数値座標変換による 比抵抗法 3 次元データ解析に関する研究 (ページ 55-63)