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

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)上式右辺第 1 項の密度変化に関し、密度流を考慮する場合には密度を次のように取り扱うこととする。

 w   wf 1  c (3)ここに、c は間隙水中の物質の濃度、ρwf は物質の濃度が 0 のときの間隙水の密度、χ は濃度の増加に対する密度の増加を定める係数である。
このとき、式(2)の右辺第 1 項は次式となる。

nS w

 w
c
 c
 nS w  wf  nS w  w
tt
1  c 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  とすることでこれらが正の物性値となる。
したがって、式(2)は次のように書き換えることができる。

 n w S w S 
 c
n 
 nS w  w
 wn w
 Swwt
1  c t
 t
 t
 S
 cn   nS w  w
  w  n w  S w
1  c t  t (7)

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

1

こととする。

1   w v i  c  S w
n  
 nS w  n Sw
 w xi
1  c t    t(8)右辺第 2 項の係数を 1 K b と書けば、最終的に次式が得られる。

1   w vi  c 1  nS w
 w xi
1  c t K b t(9)S1n n w  S w
Kb(10)

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

S
1n n w  Kb(11)

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

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

1.2

物質移行の支配方程式

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

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



(13)

ここに、vc は分散と拡散による物質の移動速度（流束）であり、vw は間隙中の水の速度である。なお、間隙水の平均流速 v と vw には、次の関係がある。

v  nS w v w(14)

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

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

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

  (15) nS w  w vic  nS w  w viw c  nS w  w vic
 nS w  w viwcc nS w  w viw
(16)
xi
xi
xi
xi
xi

2

上式の右辺第 2 項は、式(1)の間隙水の質量保存則、   w vi nS w  w xitと式(14)を考慮して、

c

(17)

 n w S w viw
nS w  w  cxi
t(18)であることから、式(16)の右辺は次のように書き換えられる。

   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(19)

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

 nS w  w

c  nS w  w vic
c nS w  w viwtxi
xi(20)あるいは、次のように書くことができる。

c1  nS w  w vic
c

 viw
tnS w  w
xi
xi(21)

他方、分散と拡散による物質の流束 vc は、Fick の法則に従うものとする。

vic   Dij

cxi(22)

D は拡散係数であり、次のような形で記述される。

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

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

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

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

2.

離散化手法

2.1

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

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

vcc
Dx(24)

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

  E c tT N T  v dV  N   c t dVe  0e e1e 1 D E(25)

ここに、[N]は内挿関数からなるマトリクス、[∂]は空間での偏導関数マトリクスである。また、添え字の t は現在の時刻を表しており、後に示す時間差分では既知量として取り扱うことを示している。
物質流束と分散係数が節点で定義され、濃度は要素内で一定とすれば、

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

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

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

 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(28)となり、勾配マトリクス[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 項を対角化すると、

4



(29)

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

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

1 0 0
Ve W   0 1 03
0 0 1(31)である。
要素単位での合算が終わった状態では、式(20)は次のような形となる。

W  v    Cet  Cst   0c t



(32)


 D 

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

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

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

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

  me 1  Ve Dij Dij   e1 m
Ve  e 11(35)

ここに、[ ]-1 は逆行列である。このとき、それぞれの節点における物質の流束は次のようにして求める。

5v   D
c ti2.2

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

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

i cxi(37)

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

   N   dV     N  c dV   0EET

Ttee 1tee 1(38)

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

t 2.3

Cet  C stW(39)

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

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

vkx(40)

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

   N   k dV     N   dV   0Ee 1T vt eEe 1

Tte(41)式(26)から式(31)と同様な手続きを踏み、要素単位での合算が終わった状態では、式(41)は次のような形となる。

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

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

 Fet  Fstv kWt(43)

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

k m1(44)mVee 1 k eVe 1e

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

2.4

要素濃度と間隙水圧の算定式(9)の方程式では、要素内の流速の勾配は内挿関数を用いて節点値から求め、濃度や間隙水圧の時間微分は差分で近似する。

1N  c t  t  c t  1  t  t   t v   nS w
tKbt1  c tttw, nod(45) wt   f 1  c t (46)

 wt

あるいは、

1 wt

B c t t  c t  1  t t   t v  nS w
tKbt1  c tttw,nod(47)

ここに、ρw,nod は節点での密度であり、節点を囲む要素の値を要素体積の重みを付けて平均した値とする。
m

 w,nod  Ve  wee 1mVe 1(48)e

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

ct  t       T ct1N  nS w  w tnod v c t  v w t  tttnS w  w あるいは、7(49)ct  t       T ct
1tc tw tBnSv
v
tw w nodt
tnS w  w (50)

ここに、(nSwρw)nod は節点での値であり、節点を囲む要素の値を要素体積の重みを付けて平均した値とする。
mme1e1

nS w  w nod  Ve nS w  w  Ve(51)また、間隙内の流速は次式で求められる。

v w tvt
nS w

(52)

間隙水の平均流速 vt は、式(43)に示したとおり既知の間隙水圧φt より求めることができる。したがって、式(50)により次の時間の濃度を現在の間隙水圧と濃度より求めることができる

c t t  1

BnS w  w tnod v c t t  v w t   t t  c tTnS w  w t(53)濃度が求まれば、式(47)より次の時間の間隙水圧も求めることができる。

 1

 t  t  K b B wt ,nod v t  nS wt w c t  t  c t  t   ttt1  c


8(54)

3.

計算フロー
密度流の陽解法における計算フローは、次のとおりである。
要素・節点、物性、境界条件入力

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

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

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

濃度の更新
時間増分

間隙水圧の更新

物性の変更（必要があれば）

結果の出力

図-1

密度流の陽解法における計算フロー9