変形・地下水・熱の連成問題に関する支配方程式と陽解法
地層科学研究所

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

 1
dV    w d w    w d  dT V Kw
(3)が成り立つ。ここに、Kw は間隙水の体積弾性定数、εw は間隙水のひずみ、ξ は間隙水の熱（体積）膨張係数である。
したがって、式(2)は次のように表される。

n w S w nS S
Tn w w
  w nS w
 n w w  S w  w
tK w t
t
t
t(4)

さらに、地盤の構成粒子の体積弾性定数が極めて大きいと仮定すると、空隙率の変化は地盤の体積ひずみ変化 εv として計測される。したがって、

n w S w nS ST
 w w
  w nS w
 n w w  S w  w v
t
K w tttt(5)上式の右辺第 2 項を、次のように書き換える。

n w

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

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

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

(7)

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

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

(8)

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

1  w vi  1 T
 nS w
 Sw v
 w xiK b ttt(9)nSS1 w n w
Kb K w
(10)

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

1  w vi  1 T n
 v
 w xiK b ttt(11)

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

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

ここに、kij は間隙水の透水係数である。

2

1.2

熱移動の支配方程式

地盤中を熱が移行する際には、熱量の保存則が満たされると考える。これは次のように表すことができる。

 Cv



 viwTT qi nS wCvw
t xi
xi



(13)

ここに、T は地盤の温度、q は熱伝導による熱流束、Cv は地盤の単位体積あたりの熱容量（体積熱容量）、Cvw は間隙水の単位体積あたりの熱容量、vw は間隙中の水の速度である。ただし、地盤と間隙水の温度は常に同じであると仮定する。
地盤の熱容量は次式でモデル化する。

Cv  Cvd  nS wCvw

(14)

ここに、Cvd は乾燥した地盤の単位体積あたりの熱容量である。
式(13)で、右辺第 1 項は熱伝導により単位時間に単位体積中で生じた熱量の流入出量の差を表しており、第 2 項は間隙水により運ばれることで同様に生じた熱量の流入出量の差を表している。また、左辺はこの熱量の変化が温度濃度の変化をもたらすことを表している。
なお、間隙水の平均流速 v と vw には、次の関係がある。

v  nS w v w(15)

したがって、式(13)は間隙水の平均流速を用いて書くことができ、

 Cv

viT T qi Cvw
t xixi(16)

さらに、次のように変形できる。

 Cvv
T qiT Cvw vi
 CvwT it xi
xi
xi(17)

上式で、右辺第 2 項は間隙水の流れにより運搬された熱量であり、第 3 項は間隙水の流入出量の差によって生じた熱量の変化を表している。
熱伝導による熱流束は、Fourier の法則に従うものとする。

qi  ij

Txi(18)ここに、λ は地盤の熱伝導率である。

1.3

構成方程式

地盤では、いわゆる有効応力則にしたがって変形が生ずると考える。

 m  K d  v   vi   (19)3 viT  fttt(20)ここに、σm は平均応力、Kd は排水条件で計測された地盤の体積弾性定数、εvi は非弾性体積ひずみ、η は地盤の熱膨張率、εf は凍結膨張ひずみである。

1.4Kb の求め方

支配方程式に現れる物性値のうち、空隙弾性定数 Kb は直接実験で求めることができないが、次のようにして試験結果に基づき値を定めることができる。なお、以下の説明では温度は一定に保たれているものとする。
非排水、飽和条件下で等方三軸圧縮試験を行い、このときの間隙水圧が計測され、平均応力増加に対する間隙水圧の増加比率 β が求められたとする。


 m
t
t(21)

他方、式(11)からは飽和で非排水条件であることを考慮して次式が得られ、1  w vi  1   v0
 w xiK b t
t


(22)


 Kb vtt(23)

また、構成方程式より次式が得られる。

 m
 Kd v ttt(24)式(21)を式(24)に代入すると次式となり、

1 
 Kd v  ttt(25)



Kd v
t 1  t(26)

式(23)と式(26)の係数を比較すると、次の関係があることがわかる。

Kb Kd1 (27)このように、空隙弾性定数 Kb は、排水条件での三軸圧縮試験で求められた体積弾性定数
Kd と、非排水条件下での等方三軸圧縮試験で求められた、平均応力増加に対する間隙水圧の増加比率 β より求めることができる。

41.5

凍結の取り扱い

いま、地盤の温度が 0℃となり、さらに冷却によって単位時間に単位体積あたり ΔQ の熱量が奪われていく場合について考える。ΔQ のうち μ の割合が潜熱 Ql により消費され、残りが温度変化に寄与すると仮定する。これは、地盤の構造骨格を介した熱移動と間隙水の凍結時の潜熱をモデル化したものである。このときの熱移動の支配方程式は次式となる。

 Cv

 qv TT 1    i  Cvw vi
 CvwT i txi
xi  xi(28)あるいは、次のように書くことができる。

Cv T  1   Q(29)

この変化が t 時間継続し、凍結終了温度 Tf まで温度が変化し、その後は潜熱の影響は無くなると考える。この間の温度変化は、次式で求められる。

T f  Q  1   t Cv(30)

他方、潜熱分が全て消費される時間は、次のとおりである。

t  Ql Q(31)

これらの式より、潜熱に消費される割合 μ は次のように求められる。

Ql
T f Cv  Ql(32)したがって、実験により Tf が求めることができれば、μ を定めることができる。
なお、地盤の潜熱は水分量の変化を考慮し、次のように求められる。

Ql  nS wQls(33)ここに、Qls は単位体積あたりの水の潜熱である。
また、凍結膨張ひずみ εf は、飽和時におけるその最大値 εf0 に対して、潜熱に消費された熱量と潜熱の割合で生ずることと仮定する。

 f  S w f 0 Q   Ql(34)

書き直すと次式となる。

 ft

S w f 0   qi
v T Cvw vi
 CvwT i Ql  xi
xi
xi 

5(35)

2.

離散化手法

2.1

Fourier 則の離散化と節点熱流束の算定

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

qTx(36)

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

   N    dV     N  T dV   0Ee1T qt eEe1Tt(37)e

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


 ET
 qt Tt


NNdV
  e     N   T dVe  0
e 1   e1E(38)

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

  (39)(40)E ET qt T tN T T t dVe  0NNdVNTdVeee 1 e1  e1E

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

E ET qt T tN T T t dSe  0NNdVBTdV  e  ee 1 e 1  e1E

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

6(41)E  qt  ET tN T T t dSe  0WBTdV   ee 1 e 1    e1E

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

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

W  q   Qet  Qst   0t(44) 

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

qt   Qet  QstW(45)

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

 m1Vemee 1(46)  Ve 1e

ここに、m は節点の周りの要素の数、λe は要素の熱伝導率である。

2.2

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

i Txi(47)

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

7EE N T  t dV  N T T t dVe  0
e e1e1 (48)

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

 Qet  Qst Wt

2.3(49)Darcy 則の離散化と節点水流速の算定

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

v
kx(50)

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

t ETv TtNdV k  e     N    dVe  0e 1   e1
E(51)

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

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

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

vt  k Fet  FstW(53)

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

81k m(54)mVee 1 k eVee 1

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

2.4

要素温度と凍結ひずみの算定式(17)の方程式では、要素内の熱流速の勾配は内挿関数を用いて節点値から求め、時間微分は差分で近似する。

Tt  tあるいは、T  T t1N  q t  Cvw v t t  CvwT t N  v t
tCvt  t   (55)T t1
B q t  Cvw v t t  CvwT t B v t
tCv  (56)

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

T t  t      1
B q t  Cvw v t t  CvwT t B v t t  T t
Cv(57)

ただし、温度が 0℃から凍結終了温度までは潜熱の影響を考慮する。

T t t  

1    Bq t  C v t t  C T t Bv t t  T tvwCv(58)vw

式(35)に示す凍結ひずみも同様であり、次のようにして既知の値より求めることができる。

 tft S w f 0 Bq  C v   C T Bv t  QCtttvw

l2.5ttvw

tf(59)v

要素応力・ひずみ・間隙水圧の算定式(11)の方程式においても、要素内の流速の勾配は内挿関数を用いて節点値から求め、時間微分は差分で近似する。

N  wt ,nodv t   1 t

1w

Kbあるいは、

1tw

B





 t t   vt
tT t t  T t n vtttt  t vt t   vt
1  t t   tT t t  T tv  nKbt
t
tttw,nod

(60)

(61)

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

9mm

 w,nod  Ve  we

Ve 1(62)ee 1

なお、間隙水の密度は Kw と ξ が一定と近似して次のようにして求めることとする。

 wt   w0   w   w0 1 
1 t  0   T t  T0 Kw (63)ここに、ρw0、φ0、T0 は、それぞれ t=0 のときの間隙水の密度、間隙水圧、温度である。
式(61)を変形すると次式となる。

 1 w

 t t  K b B wt ,nodv t t  n T t t  T t     vt t   vt    tt

(64)

上式と構成方程式を併せて解くことによって、各時刻の値を得ることができる。

 mt t  K d  vt t   vit t    t t(65)

計算方法はつぎのとおりである。まず、式(64)を二つに分ける。

 t t  a  b   t

(66)

(67)

a  K b  vt t   vt

 1
b  K b  t B  wt ,nodv t t  n T t t  T t w(68)Δφa は地盤の変形によってもたらされた間隙水圧変化であり、Δφb は水の流れや温度変化によってもたらされる間隙水圧変化である。言い換えれば、水の流れや温度変化を発生させないで地盤を変形させた場合は Δφa の間隙水圧変化が生じ、変形を拘束して水の流れや温度変化を生じさせた場合は Δφb の間隙水圧変化が生ずることとなる。そこで、水の流れや温度変化を発生させないステップと、変形を拘束して水を流し温度を変化させるステップに分けて計算することとする。
水の流れを発生させないステップでは、式(65)と式(67)を動的緩和法で解く。水の流れを考慮しないため、時刻を t に固定して（時間に係わらず）計算することとなる。動的緩和法
における n 回目の繰り返し計算では、式(65)の応力を内力として計算を進め、変形からひずみが求まる。

 mt n   K d  vt n    vit    t n (69) vt n1  BU t n1

(70)10なお、式(57)や式(59)に示したとおり、温度や凍結ひずみなどは別途求められていることから、εvi は既知項として取り扱うことができる。

 vit t  T t t  T t    tft   tf    vit(71)

次に、式(67)からひずみ変化に伴う間隙水圧変化を求める。時間を止めているが、動的緩和法の計算過程で発生したひずみ変化に対しても式(67)は適用される。

a  Kb  vt n1   vt n (72)

間隙水圧を更新し、再度繰り返し計算を行い、

 t n1  a   t n (73) mt n1  K d  vt n1   vit    t n1(74)

応力の変化が十分小さくなった段階で時刻 t の解とする。
続いて、変形を拘束して水を流し温度を変化させるステップの計算を行う。この場合は、
時刻 t の間隙水圧をもとに、時刻が Δt 進んだ後の間隙水圧変化を求める。

T t  t      1
B q t  Cvw v t t  CvwT t B v t t  T tCv

 1b  K b  t B  wt ,nodv t t  n T t t  T t w t t  b   t

(75)

(76)

(77)

得られた間隙水圧変化を式(65)に加え動的緩和法を適用することにより、時刻 t+Δt における応力やひずみ、間隙水圧が求められる。計算のフローを以下に示す。

113.

計算フロー
陽解法による変形・地下水・熱の連成解析の計算フローは、次のとおりである。
要素・節点、物性、境界条件入力
Δt1（動的緩和法）とΔt2（連成時間増分）の設定
節点荷重項の設定
表面力＋物体力‐要素応力等価節点力
要素応力：
変位の算定

ひずみの算定動的緩和法
既知変位の更新
間隙水圧の算定

応力の算定

No
Yes
温度の更新

非弾性ひずみの更新
連成時間増分間隙水圧の更新

結果の出力

図-1

変形・地下水・熱の連成問題の陽解法における計算フロー12