一次元移流分散方程式の差分計算

第 1 章の緒論で述べたように，一次元移流分散方程式はダルシー則と併せて地下水流動や 物質移行挙動の評価に用いられてきた．本研究においては，式 3.1 に示すように，無次元化 した一次元移流分散方程式を用いた．また，式 3.2 に示すペクレ数（Pe [-]）は，移流によ る物質移行速度と分散による物質移行速度の比を示す無次元数である[1-4]．なお，C は無次 元濃度 ( = c/c^*) [-]，Xは無次元位置 ( = x/x1) [-]，Rdは遅延係数 [-]，Tは無次元時間 ( = t/t^*)

[-]，u はダルシー流速 [m/s]，x1は代表距離 [m]，Deは分散係数 [m²/s]，x は位置 [m]，c は

濃度 [mol/m³]，c^*は代表濃度 [mol/m³]，tは時間 [s]，t^*は平均滞在時間 [s]を示す[5-7]．

2 2 d e

1

C C C

X P X R T

  

  

   (3.1)

1 e

e

P ux

 D (3.2)

3.2.1 移流項の差分精度が計算結果へ及ぼす影響

一次元移流分散方程式は解析解が求められているものの，その導出に用いる物質濃度の初 期条件や境界条件が限られている[8]．それらの初期条件や境界条件は，実験条件のみならず，

モデル化する現象ごとに変化するため，その都度，解析解を導出する必要がある．しかしな がら，初期条件や境界条件が変化するごとに，解析解を求めることは現実的でないため，一 次元移流分散方程式を解く際には計算機による数値計算により求める．

偏微分方程式の数値計算においては，偏微分項を離散化することで，離散方程式として近 似値を計算する．本研究においては，偏微分方程式の離散化に差分法を用いた．差分法は，

以下の式 3.3 および式 3.4に示すTaylor 展開を利用して偏微分式を離散近似する方法である

[9, 11]．なお，O(Δx)^mはTaylor展開における打切り誤差を示し，mは打切り精度の次数を示

す．

2 2 3 3

4

2 3

( ) ( )

( ) ( ) ( )

1! 2! 3!

x f x f x f

f x x f x O x

x x x

     

       

   (3.3)

35

2 2 3 3

4

2 3

( ) ( )

( ) ( ) ( )

1! 2! 3!

x f x f x f

f x x f x O x

x x x

     

       

   (3.4)

式 3.1 に示した一次元移流分散方程式には，一階の偏微分項である移流項（左辺第一項）

と二階の偏微分項である分散項（同第二項）がある．この二階の偏微分項である分散項は，

式3.3および式3.4の両辺の和をとることにより，式3.5のように差分化される．この式は二 階の偏微分項における中心差分を意味し，打切り誤差はO(Δx)²であるために，この差分は2 次精度であることを意味する．しかし，式3.3および式3.4の和をとることは，3階の偏微分 項が相殺されていることから，見かけ上は 3 次精度の差分式であり，このために式 3.5 の差 分近似は良い計算精度を示す．

2

2 2

( ) 2 ( ) ( )

( )

f f x x f x f x x

x x

       

  (3.5)

一方で，移流項における一階の偏微分式は，式 3.3 または式 3.4 をそれぞれ式変換するこ とで以下のような，前進差分（式 3.6）と後進差分（式 3.7）が導かれる．また，式3.3 と式 3.4の差より，式3.8の中心差分が得られる．これらの差分は1次精度であり，差分の計算に おいて誤差が大きいことが知られる．

( ) ( )

f f x x f x

x x

    

  (3.6)

( ) ( )

f f x f x x

x x

    

  (3.7)

( ) ( )

2

f f x x f x x

x x

      

  (3.8)

ここで，風上差分の精度を確認するために，式3.9に示す移流方程式について，右辺を中 心差分により式3.10のように差分化する．なお，式3.10の左辺はExplicit Methodにより時 間ステップnとそのΔT後の時間ステップn+1において差分化しており，濃度（C）に付記 された添え字は位置を意味し，位置iを中心とした風上側i-1と風下側i+1の濃度を示す．

式3.10を変形すると，ΔT後の中心濃度は式3.11のようになる．なお，ωはクーラン数と 呼ばれる無次元数であり，差分計算が安定するためには1以下に設定する必要がある．この 無次元数は，単位時間メッシュ（ΔT）あたりにおける移流による物質移動（UΔT）が，空

36

間メッシュ（ΔX）を超えると数値的に安定しないことを意味する．

C C

T U X

   

  (3.9)

n+1 n n n

i i i+1 i-1

2

C C C C

T U X

 

    (3.10)

n+1 n n+1 n

i i+1 i i-1

1 1

2 2

C   C C  C (3.11)

U T

 ^X

 (3.12)

式3.11 より，U＞0の時に右辺の Ci+1nの係数が負となるために，左辺の値が数値的に振動 する可能性がある．この数値的な振動は，計算結果において物質濃度が負となることを意味 し，数値計算が成り立たない場合が生じる．このような場合，移流方程式や移流分散方程式 の差分計算において，一階の偏微分式である移流項における計算精度を向上する方法の一つ として，風上差分が用いられる．風上差分は，式3.13のように，移流項において常に風上側 に重みづけされた差分式を示し，式 3.11 における風下側の濃度（Ci+1n）を用いないために数 値計算の安定性が増す．そして，式 3.13 を U の正負によらない形式により表すと，式 3.14 のようになる．

i i-1

i+1 i

, ( 0)

, ( 0)

C C

C X U

U X C C

X U



 







  



 



 (3.13)

i+1 i-1 i+1 i-1

( ) ( 2 )

2 2

C U U

U C C C C C

X X X

     

   (3.14)

i+1 i-1 i+1 i-1

2

( ) ( 2 )

2 2 ( )

U X

C C C C C

U C U

X X X

   

  

   (3.15)

式3.14は，式3.15のように変形すると，左辺の中心差分に分散項を|U|ΔX/2倍した項を加

37

えた式となる．この，風上差分により加えられた見かけ上の分散項は，数値分散（数値拡散）

と呼ばれる．この数値分散を加えることにより，移流項の数値計算的な安定性は増すものの，

移流分散方程式においては，流速の増加に伴って分散が過大評価されることとなる．よって，

無次元化した一次元移流方程式を用いる際には，ペクレ数が非常に大きい場合や，移流が非 常に速いために局所収着平衡が成り立たない場合のような解析の際に，移流項における数値 拡散が問題となる．

この，移流項の計算精度が低くなる問題については多くの先行研究で議論されており，例 えば，空間メッシュ間の濃度を補完する QUICK 法[12]，特性曲線法[13, 14]，差分精度の打 切り誤差を高精度に行うKKスキームなどが対応策として示されつつある[15]．

本研究においては，3 次精度の風上差分を移流項に適用することで，この問題に対処する こととした．式 3.16 に示すように，3 次精度風上差分は着目する位置 i，風上側の 2 点，風 下の1 点を用いて一階の偏微分式を近似する．ここで，中心の濃度に対して風上側に2 点の 値を用いるために，境界条件における仮想点が1次精度の際より多く必要になる．

i-2 6 i-1 3 i 2 i+2

1 6

C C C C

C

X X X

  

 

   (3.16)

ここで，一次元移流分散方程式に 3 次精度の風上差分を適用し，1 次精度の風上差分と解 析解による計算結果と比較することとした．なお，式3.17に示す解析解はOgataらによるも のであり，境界条件は式3.18のClosed Vessel B.C.，初期条件は実験条件と同じステップ入力 である[8, 16]．また，Closed Vessel B.C.は1次精度と3次精度のいずれにおいても，中心差分 により与えた．この際，三次精度の差分計算を行うためには，一時精度の場合よりも計算領 域外における仮想点を多く定義する必要がある．この仮想点の定義については Appendix に 示す．

e e

e

( ) 1 (1 ) exp( ) erfc (1 )

2 4 4

P P

C T erfc T P T

T T

    

       

   

    

 

(3.17)

x 0 x 0

x=0+

e

1 C

C C

P P

   

  

 (3.18)

38

図3.1 移流項における差分に1次精度または3次精度を適用した際の

一次元移流分散方程式の数値計算結果と解析解との比較

図3.1に，移流項における差分の精度が計算結果に及ぼす影響を示した．なお，ペクレ数

は10，25，50および100の四通りに設定し，一次精度風上差分と三次精度風上差分を適用

した数値計算結果を解析解による計算結果と比較した．ペクレ数が10および25の場合に は，移流項における差分の精度に依存した計算結果への誤差は生じなかった．一方で，ペク レ数が50を超えると，一次精度の風上差分による計算結果は解析解から乖離した．しか し，3次精度の風上差分を適用することで，ペクレ数が100を超えるような，移流による物 質移動が卓越した条件においても数値拡散による誤差が制限され，解析解と変わらない計算 結果が得られることが示された．なお，本論文における数値計算はFortran90/95により行 い，そのコンパイルはNumerical Algorithms Group（NAG）のFortran Builder 6.0 for Windows により行った．

39

ドキュメント内 放射性廃棄物の浅地中処分を考慮した不飽和層の物質移動の評価に関する基礎的研究 (ページ 46-51)