地下水流れに関する支配方程式と陽解法
地層科学研究所

1.

支配方程式
地盤中を地下水（間隙水）が移動する際には、水の質量保存則が満たされると考える。

これは次のように表すことができる。

   w vi n w S w xit(1)

ここに、ρw、Sw はそれぞれ間隙水の密度と飽和度、n は地盤の空隙率である。また、v は間隙水の平均流速であり、ある断面より単位時間に流出（入）した間隙水の体積を、その断面積で除して求められる。
上式で、左辺は単位時間に単位体積中で生じた間隙水の流入出質量の差を表しており、右辺は、この流入出質量の差が地盤の空隙率の変化と間隙水の密度変化、及び飽和度の変化によって補われることを示している。
式(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
t
K 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)は、次のように書き換えることができる。

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

(7)

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

S
1  w vi   nS w
n   
 n w  Sw w xi  t
 Kw(8)右辺の係数を 1 K b と書けば、次式が得られる。

1  w vi  1  w xi
K 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 の法則が成り立つものとする。

v i  k r

 ij  x  w g 3   x j
x j (14)上式で kr は相対透過係数、κij は地盤の固有透過度テンソルである。また、μ は間隙水の粘性係数である。地盤の透水係数を kij とすれば、次式となる。

kij  k r

 ij
(15)

 x v i  k ij  w g 3  x
x j 
 j(16)

22.

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

vkx(17)

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

t ETv TtNdV k  e     N    dVe  0e 1   e1E(18)ここに、[N]は内挿関数からなるマトリクス、[∂]は空間での偏導関数マトリクスである。また、t は現在の時刻を表しており、後に示す時間差分では既知量として取り扱うことを示している。
流速と透水係数が節点で定義され、間隙水圧は要素内で一定とすれば、 ETvt TtNNdVe     N    dVe  0
ke 1   e1E(19)ここに、{ }は節点ベクトルである。変形すると次式となる。

  (20)





(21)

E
 ET
vt T t

N T  t dVe  0NNdVNdV  e  ee 1 e 1k  e1E

また、発散定理より上式は、
E
 ETvt T t
N T  t dS e  0NNdVNdV  e  ee 1 e 1k  e1Eとなり、勾配マトリクス[B]を用いると次式が得られる。

E ETvt T tN T  t dS e  0NNdV
BdVeee 1 e 1k  e1E

左辺第 3 項は、境界面の水圧と等価な節点流速である。
左辺第 1 項を対角化すると、

3(22)





E
 vt  ET tN T  t dS e  0WBdV   ee 1 e 1 k   e1E

ここに、[W]は(23) N  N dVe を対角化したマトリクスである。例えば、2 次元の三角形要素Tで内挿関数が 1 次の要素では次のようになる。

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

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

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

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

vt  k Fet  FstW(26)

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

k m1Vee 1 k e(27)mVe 1e

ここに、m は節点の周りの要素の数、ke は要素の透水係数である。
式(12)の方程式では、要素内の流速の勾配は内挿関数を用いて節点値から求め、時間微分は差分で近似する。

N v t   1 あるいは、

Kbttt  t

1  t  t   t
B v Kbt

t

(28)

(29)

したがって、既知の値から次の時間の間隙水圧を求めることができる。

 t t  K b Bv t t   t(30)

43.

計算フロー
地下水流れの陽解法における計算フローは、次のとおりである。

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

Δt の設定間隙水圧からの節点流速の算定

間隙水圧の更新
時間増分

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

図-1 地下水流れの陽解法における計算フロー

5