物質移行に関する支配方程式と陽解法
地層科学研究所

1.

支配方程式

1.1

地下水流れの支配方程式

地下水流れの支配方程式は、次のとおりである。

  w vi n w S w xi
t(1)ここに、ρw、v 、Sw、n はそれぞれ間隙水の密度、平均流速、飽和度、地盤の空隙率である。
式(1)の右辺は、次のように分解することができる。

n w S w Sn nS w w  n w w  S w  w
t
t
t
t

(2)

間隙水では間隙水圧φ、体積を V とすると、質量が一定であることから、d w  w
VdV    w d w    w

1dKw(3)が成り立つ。ここに、Kw は間隙水の体積弾性定数、εw は間隙水のひずみである。したがって、式(2)は次のように表される。

n w S w nS Sn w w
 n w w  S w  w
tK w tt
t(4)上式の右辺第 2 項と第 3 項を、次のように書き換える。

n w

S w
S 
 n w wt t(5)nn  Sw wt t(6)Sw w

ここに、 S w  は比水分容量、 n  は比貯留係数と呼ばれる物性値である。なお、ここでは間隙水圧は応力と同様に圧縮を負とするため、比水分容量を  S w  、比貯留係数を  n  とすることでこれらが正の物性値となる。
式(4)は、次のように書き換えることができる。

n w S w nS S n  w w wn w
 Sww
tK w t
 t
 t
 nSSn  
  w  w  n w  S w
  t
 Kw

(7)

最終的に、式(1)は次のようになる。なお、両辺は密度で除して体積変化に変換しておく。

1S
1  w vi   nS w
n  
 
 n w  Sw
 w xi  t
 Kw

(8)

右辺の係数を 1 K b と書けば、次式が得られる。

1  w vi  1  w xiK b t(9)nSS1n w  n w  SwKb K w
(10)

ここに、Kb は空隙弾性定数であり、水の入出量差によって生じた体積変化に対する間隙水圧変化を定める。また、不飽和領域では間隙水圧変化が空隙率に影響を及ぼさないと考えられることから、飽和領域で 1、不飽和領域で 0 となる係数 α を用いて式(10)を書き直す。

nS
S1n w  n w KbKw
(11)

なお、ここでは間隙水の密度変化が充分小さいことを前提とし、次式を採用する。

vi
1 
xi K b t(12)S1n n w  Kb(13)

他方、間隙水の平均流速に関しては、Darcy の法則が成り立つものとする。

 x v i  k ij  w g 3  x
x j  j(14)ここに、kij は地盤の透水係数である。

1.2

物質移行の支配方程式

地下水の流れとともに物質が移行する際には、物質の質量保存則が満たされると考える。
これは次のように表すことができる。

 nS w  w c   nS w  w vic  nS w  w viw ctxixi

(15)

ここに、c は間隙水中の物質の濃度である。また、vc は分散と拡散による物質の移動速度（流束）であり、vw は間隙中の水の速度である。
上式で、右辺第 1 項は分散と拡散により単位時間に単位体積中で生じた物質の流入出量の差を表しており、第 2 項は間隙水により運ばれることで同様に生じた物質の流入出量の差を表している。また、左辺はこの物質量の変化が濃度の変化をもたらすことを表してい2る。なお、間隙水の平均流速 v と vw には、次の関係がある。

v  nS w v w(16)

式(15)で、左辺は次のように分解することができる。

nS w  w c nS w  w c c nS w  wttt

また、右辺は次のように分解できる。

  (17) nS w  w vic  nS w  w viw c  nS w  w vic
 nS w  w viw
c(18)
c
 nS w  w viw
xi
xi
xi
xi
xi
上式の右辺第 2 項は、式(1)の間隙水の質量保存則、   w vi nS w  w xitと式(16)を考慮して、c(19) n w S w viw
nS w  w  cxi
t(20)であることから、式(18)の右辺は次のように書き換えられる。

   nS w  w vic  nS w  w viw c  nS w  w vic
nS w  w cc
 nS w  w viw
xi
xi
xi
t
xi

(21)

改めて式(15)を書き直すと、式(17)の右辺第 1 項と式(21)の右辺第 2 項が消去されることから、最終的に次式が得られる。

c  nS w  w vic
c nS w  w
 nS w  w viw
t
xi
xi

(22)

あるいは、次のように書くことができる。

c
1  nS w  w vicc viw
tnS w  wxixi(23)さらに、間隙水の密度変化が十分小さいと仮定すれば、質量保存則は最終的に次式となる。

viccc
 viw
txixi(24)他方、分散と拡散による物質の流束 vc は、Fick の法則に従うものとする。

vic   Dijcxi(25)D は拡散係数であり、次のような形で記述される。

Dij   T v  ij   L   T w

viw v wjvw

 Dm ij3(26)ここに、αT とαL は横分散長（間隙内流速方向と垂直方向の分散長）と縦分散長（間隙内流速方向の分散長）、Dm は分子拡散係数、τ は屈曲率、δij はクロネッカのデルタである。また、 vw

は間隙内水流速ベクトルの長さであり、2 次元であれば

4v   v  である。
w 2xw 2y2.

離散化手法

2.1

Fick 則の離散化と節点物質流束の算定

要素分割した領域に対して、内挿関数を重みとして式(25)にガラーキン法を適用する。最初に、式(25)を次のように変形する。わかりやすくするため添え字は省略した。

vcc
Dx(27)

上式に内挿関数を乗じ、対象領域で積分すると次式が得られる。ただし、積分は要素単位で行い節点で合算することとする。

  E
 c t
T
 N T  v dV  N   c t dVe  0
e D  e1e 1E(28)ここに、[N]は内挿関数からなるマトリクス、[∂]は空間での偏導関数マトリクスである。また、添え字の t は現在の時刻を表しており、後に示す時間差分では既知量として取り扱うことを示している。
物質流束と分散係数が節点で定義され、濃度は要素内で一定とすれば、

 
 E
 vc t Tt
 N T N dV e     N   c dVe  0e 1
 D  e1E(29)

ここに、{ }は節点ベクトルである。変形すると次式となる。

 E

 E
 vc t 
T tT N T N dV
NcdV N  c t dVe  0ee e1
e 1e 1 D E  (30)また、発散定理より上式は、

 E E
 vc t 
T t N T N N T c t dS e  0dVNcdV e  eDe 1e 1 e1E(31)となり、勾配マトリクス[B]を用いると次式が得られる。

 E

 E
 vc t 
T t N T N 
N T c t dS e  0dV

BcdV e  ee 1e 1 D  e1E

左辺第 3 項は、境界面の濃度と等価な物質の節点流束である。
左辺第 1 項を対角化すると、

5

(32)

 tE vc   ET t W N T c t dS e  0
BcdV  ee 1e 1 D   e1 E

ここに、[W]は





(33)

 N  N dVe を対角化したマトリクスである。例えば、2 次元の三角形要素
Tで内挿関数が 1 次の、いわゆる定ひずみ要素では、

1 0 0
Ve W   0 1 0
3
0 0 1

(34)

である。
要素単位での合算が終わった状態では、式(33)は次のような形となる。

W  v    Cet  Cst   0c t(35) D 

ここに、{Cet}は各節点における等価濃度勾配、{Cst}は境界面の濃度と等価な節点濃度勾配である。
左辺第 1 項の係数は対角マトリクスであることから、式(35)は節点ごとの並列式となり、各節点での割り算によりそれぞれの節点における物質の流束を求めることができる。

v   D  CW Ctec tts(36)なお、拡散係数は要素の物性として定める。節点での拡散係数は、節点の周囲の要素について、拡散係数の逆数の要素体積の重み付き平均を求め、この逆数をもって節点値とする。

D mVee 1 De1(37)mVe 1e

ここに、m は節点の周りの要素の数、De は要素の拡散係数である。
拡散係数は一般にはテンソルであることから、次のように節点値を求める。

  me 1  Ve Dij Dij   e1 mVe  e 11(38)ここに、[ ]-1 は逆行列である。このとき、それぞれの節点における物質の流束は次のようにして求める。

6v   Dc ti2.2

    Cet j  C stij(39)W

節点濃度勾配の算定式(24)に示す保存則を解くためには、濃度勾配が必要となる。ω を濃度勾配として、次式をガラーキン法により離散化する。

i cxi(40)

上式に内挿関数を乗じ、対象領域で積分すると次式が得られる。ただし、積分は要素単位で行い節点で合算することとする。

   N   dV     N  c dV   0EETTtee 1tee 1(41)

式(29)から式(35)と同様の手続きを踏むことで、それぞれの節点での値が求まる。

t 2.3Cet  C stW(42)

Darcy 則の離散化と節点水流速の算定

要素分割した領域に対して、内挿関数を重みとして式(14)にガラーキン法を適用する。なお、説明を簡単にするため右辺第 2 項は省くこととする。また、わかりやすくするため添え字は省略した。
最初に、式(14)を次のように変形する。

vkx(43)

上式に内挿関数を乗じ、対象領域で積分すると次式が得られる。ただし、積分は要素単位で行い節点で合算することとする。

   N   k dV     N   dV   0Ee 1T vt eEe 1Tte(44)式(29)から式(35)と同様な手続きを踏み、要素単位での合算が終わった状態では、式(44)は次のような形となる。

W  v   Fet  Fst   0t(45)k 

ここに、{Fet}は各節点における等価水圧勾配、{Fst}は境界面の水圧と等価な節点水圧勾配7である。
左辺第 1 項の係数は対角マトリクスであることから、式(45)は節点ごとの並列式となり、各節点での割り算によりそれぞれの節点における流速を求めることができる。

 Fet  Fstv kWt(46)

なお、透水係数は要素の物性として定めることとする。このため、節点での透水係数は、節点の周囲の要素について透水係数の逆数の要素体積の重み付き平均を求め、この逆数をもって節点値とする。

k m1(47)mVee 1 k eVee 1

ここに、m は節点の周りの要素の数、ke は要素の透水係数である。

2.4

要素濃度の算定式(24)の保存則において、要素内の物質流束の勾配は内挿関数を用いて節点値から求め、時間微分は差分で近似する。

ct  t              tt T ct
  N  vc  v w
t
t(48)

あるいは、

ct  ttt T ct B v c  v w
tt(49)

なお、間隙内の水流速は次式で求められる。

v   nSvw tt(50)w

したがって、既知の値から次の時間の要素の濃度を求めることができる。

  t  v    t  cc t t  B v ctTw tt8t(51)

3.

計算フロー
物質移行の陽解法における計算フローは、次のとおりである。

要素・節点、物性、境界条件入力

Δt の設定濃度からの節点物質流束の算定

移流項のための節点濃度勾配の算定

濃度の更新

時間増分

物性の変更（必要があれば）結果の出力

図-1 物質移行の陽解法における計算フロー9