飽和・不飽和土の水・空気・熱連成解析手法の研究
張 鋒*・岩井裕正**・彭鵬***・熊勇林****1
1.

研究の目的

放射性廃棄物処分における坑道の人工バリア設計においては，空洞変形，廃棄体の熱と再冠水時の不飽和状態から飽和状態への移行による人工バリアに及ぼす影響を評価する必要がある．そこで、
本研究は、人工バリア材料の熱挙動および飽和/不飽和特性を評価できる解析技術を確立することを目的とし，名古屋工業大学地盤研究室で開発した飽和/不飽和土･水･空気･熱連成 FEM 解析（THMAcoupling analysis）プログラム Soft(Xiong et al., 2014)を，室内試験、模型実験および原位置試験のシミュレーションを実施することで，その適用性を検証する．
2.

研究の方法

2019 年度においては、名古屋工業大学で実施したトンネル模型実験に使用される人工軟岩の三軸圧縮･クリープ要素試験と模型実験の再現解析を実施し、温度の変化による影響を考慮した連成解析の有効性を検証した．2020 年度では、原位置試験の再現解析および実務で解析手法の適用を前提とした提案解析手法の検証を行う．また、同じ解析手法による掘削問題や不飽和問題の事例解析を行い、解析プログラムの信頼性と適用性を深める．さらに、実施工を考慮した解析や、軟岩を対象として温度を考慮した原位置試験への適応解析も行い、構成モデルおよび解析プログラムの汎用性を検証する．本研究で実施したすべての数値解析は著者が提案してきた精巧な熱･弾塑性飽和/不飽和土構成式（Zhang & Zhang, 2009; Zhang & Ikariya, 2011; Zhang et al., 2012）に基づいていることを付け加える．
2.1

不飽和地盤での加熱試験における THMA 連成解析

加熱実験は，図 1 に示すように，モンテリ地下岩石研究所（Mont Terri Underground RockLaboratory）で実施された大規模な原位置試験である（Munoz，2006）．加熱実験の目的は，THM 連成プロセス内の母岩とベントナイトバッファー（人工バリア）の挙動を調べるである．加熱試験場は，地下岩石研究所のニューギャラリーの西壁にある「HE ニッチ」の Opalinus 粘土層の頁岩相にありました．図 2 に，モンテリ地下岩盤研究所の鳥瞰図と HE ニッチの位置を示す．図 2 に示すように，ニッチフロアに直径 0.30m，深さ 7.0m の BHE-0 と呼ばれる垂直ボアホールを掘削した．ボアホール
BHE-0 の軸に外径 100mm の加熱管を配置した．加熱管に直径 75mm，長さ 2.0m のヒーターを設置し，試験装置の詳細を図 3 に示す．
加熱実験は 3 つの異なる段階で実施された．フェーズ１は，982 日間続いたベントナイトバッファーの給排水である．ベントナイト緩衝材の給排水は，ニッチの床と周囲の母岩の上にある 2.0m のピエゾメータヘッドを 4 つの異なる深さで与えた．フェーズ２は，ベントナイトバッファーが完全に飽和された後の加熱フェーズである．まず，ヒートバッファ接触部が 100℃の温度に達するまで，140W，150W，285W，580W の４ステップで加熱した．次に，ヒーターとベントナイトの接触部で 100℃の一定温度を 540 日間の加熱期間で維持した．最後に，加熱をやめ，冷却段階が始まり，最後まで続いた．実験では，選択された各地点で間隙水圧と温度を測定した．
数値解析においては対称であるため，対象物の半分のみを考慮する．図 4 に 2D FEM メッシュと境界条件を示す．解析では，図 4 にオレンジ色の線でマークされた各位置での THMA 連成解析の計算結果を実験結果と比較する．簡単にするために，材料が均質，等方であると想定し，岩盤母材で構成される緩衝材の THMA 連成変化に焦点を当てる．シミュレーションでは，Munoz（2006）の仮定と同じように，大きさ 5MPa の等方応力状態を初期応力状態とする．また，母岩の全水頭は 40m，初期温度は 15℃であり，試験条件と同じである．シミュレーションでは，実験と同様に空気圧は一定*

名古屋工業大学・教授，***同・助教，***同・博士後期課程 D2，****寧波大学・准教授

で，大気圧に等しいと設定した．計算で使用したベントナイトと岩石の材料パラメータを表 1 に示す．

図 1. モンテリ地下岩石研究所の鳥瞰図と HE ニッチの位置（Munoz, 2006）

図 2. ボアホール BHE-0 の垂直断面（Munoz, 2006）

図 3. 加熱試験の概観（Munoz, 2006）

図 5 に，シミュレーションで使用したベントナイトと母岩の水分特性曲線を示す．ベントナイトと母岩の計算された水分特性曲線が試験結果によく適合していることがわかります．ベントナイトと母岩の水分特性曲線のパラメータを表 2 に示す．シミュレーションでは，ベントナイトバッファ
ーの初期飽和度は 70％であり，これは現場で測定された s = 136MPa のサクションに相当する．

表 1 ベントナイトと母岩の材料パラメータ

Compression index 
Swelling index 
Critical state parameter Μ
Void ratio N (p’=98 kPa on N.C.L.)
Poisson’s ratio νParameter of overconsolidation aParameter of suction b
Parameter of overconsolidationβ
Void ratio Nr (p’=98 kPa on N.C.L.S.)
Thermal expansion coefficient (1/K)
Thermal expansion coefficient of water (1/K)
Thermal conductivity (kJ m-1 K-1Min-1)
Specific heat (kJ Mg-1 K-1)
Specific heat of water (kJ Mg-1 K-1)

Bentonite
Rock0.0050
0.00200.010
0.0001
1.80
1.90
1.04
0.62
0.30
0.30
5.0
5.0
0.00
0.00
1.0
1.0
1.06
0.65
1.0×10-5 3.0×10-6
2.1×10-4
0.06
0.127238744184

表 2 水分特性曲線に含まれるパラメータ

Saturated degrees of saturation S
Residual degrees of saturation Ssrrr

Bentonite

Rock

1.00

1.00

0.40

0.40
21000

Parameter corresponding to drying AEV (kPa) Sd11000

Parameter corresponding to wetting AEV (kPa) Sw
Initial stiffness of scanning curve (kPa) k spe

800

1000

25000

90000

Parameter of shape function c1

0.000001

0.00003

Parameter of shape function c20. 000005

0.00006

Parameter of shape function c3

30.050.0

水の透水性は飽和度に依存することはよく知られている．つまり，透水性は飽和度の増加とともに増加する．ただし，水の透水性は実験では測定されなかった．加熱試験の THMA 解析では，別の方法として，特定の飽和度での透水性の値を使用する内挿法を使用して，飽和による透水性の変化をシミュレートした．加熱実験の結果に合わせるために，水の透水性と解析に使用した飽和度の関係を図 6 に示します．
図 7 にさまざまな位置での飽和度の経時変化を示す．図から，ベントナイトを完全に飽和させるのに必要な時間は約 400 日であることがわかった．尚，図中の D はヒーターからの距離である
図 8 は，周囲の岩石の飽和度の経時変化を示しており，不飽和化と再飽和の過程がはっきりと観察されている．この図から，約 600 日間の給排水後，岩石は再び完全に飽和状態に達することがわかった．
図 9 は，ヒーターから離れたさまざまな位置での温度の変化を示している． THMA 解析では，選択したすべての位置において温度の急激な上昇や急激な低下など，HE-D 実験で測定された温度の変化を十分に説明できることがよくわかる．
図 10 に，異なる位置での計算と実験から得られた温度変化の比較を示す．図から，ヒーターからの距離が近いほど温度が高くなることがわかる．岩石の熱伝導率が低いため，ヒーターから 5m 離れた場所では顕著な温度上昇は見られない．

40 mheaterbentonite

5m2m20 m
図 4.

2D FEM メッシュと境界条件

610-8Permeability of water (m/s)

10Bentonite
5410

Nature rock310

Line: calculation
Dot: experiment-9

Host rock
Bentonite

10-1010-1110-1210-1310-14210

405060708090

100

104050Degree of saturation (%)

6080

図 6. 透水性と飽和度の関係

100

D=0.05m
D=0.10m
D=0.15m

908070070

Degree of saturation (%)

図 5. ベントナイトと母岩の水分特性曲線

Saturation (%)

Suction (MPa)

10200

400

600

800

1000

Time (days)

図 7. 給排水段階でのベントナイトの飽和度の経時変化

90100

120
D=0.05 m
D=0.65 m
D=5.00 m

100oTemperature ( C)

Degree of saturation (%)

100

D=2.00m
D=0.65m
D=0.20m9580
Line: calculation
Dot : experiment

60
40
20

90
0200

400

600

800

0600

1000

Time (days)

図 8. 給排水段階の岩石飽和度の変化

1000

1200

1400

1600

Time (days)
図 9. 加熱段階の各位置での時系列的温度変化

100

Excessive pore water pressure (kPa)

500990 days
1030 days
1200 days

80
oTemperature ( C)800

Line: calculation
Dot: experiment

60402000123

Distance from heater (m)図 10. 温度変化の実験と計算の比較

45D=0.65m

400

D=1.65m
300200100
Line: calculation
Dot : experiment0.00-100800

1000

1200140016001800

Time (days)
図 11. 加熱・冷却時各点の過剰間隙水圧の経時変化

図 11 は，加熱段階と冷却段階のさまざまな位置での過剰間隙水圧（EPWP）の経時変化を示している．この図から，温度の上昇により正の EPWP が大幅に増加することがわかる． EPWP の変化は，主に水の熱膨張係数が岩石の熱膨張係数よりもはるかに高いことに起因することは言うまでもない．
岩石の透水性が低いため，排水が遅く，間隙水の膨張が妨げられ，初期加熱時の間隙水圧が上昇する．その後，間隙水の移動に伴い，EPWP は消散し，結果として減少するようになった．シミュレー
ションは，全体として実験結果をうまく説明できていると言える．
2.2 不飽和シラス地盤の斜面崩壊に関するモデル試験のシミュレーション
本節では Kitamura et al. (2007)が実施したモデル試験をシミュレートして，提案した数値解
析手法が 3 つの異なるパターンで実施したモデル試験で観察された不飽和シラス地盤の斜面崩壊のメカニズムを説明できるかどうかを検証する．
写真 1 は，試験装置とモデル地盤の鳥瞰図を示している（Kitamura et al., 2007）．モデル斜面
地盤は，軽く締固めされた不飽和しらすでできている．モデル斜面は層ごとに慎重に締固めされた（各層の厚さは 5cm である）．地盤の物性を表 3 に示す（Kitamura et al., 2007）．モデル斜面の傾斜角は 45°，高さは 80cm である．排水境界条件を制御するために，透水性板が土壌箱の端に設置されている．
図 12 に，測定ユニット，注水位置，FEM メッシュのレイアウトを示す．ケース 3 では，人工芝シートを斜面の上面に配置し，人工降雨システムを使用して天井から水を浴びせたときに土壌が洗い流されないようにした． No.1〜15 で示される 15 個のテンシオメーターをモデル地盤内に埋めて，サクションを測定した．図 12(b)は，斜面崩壊のシミュレーションで使用された有限要素メッシュを示している． 1681 節点と 1600 の 4 節点アイソパラメトリック要素で構成される FEM メッシュのサイズは，平面ひずみ条件下でのモデル実験のサイズと同じである．境界条件は次のように与えられている：（a）変位条件の場合，x，y 方向の両方で下部に固定され，左側面のみで x 方向の垂直変位が固定され，他の表面は x，y 方向で自由である；（b）排水および排気条件の場合，図 12(a)に示
すように，上面と斜面はすべての場合に透気･透水であり，下面と左側面は水注入の範囲を除いて不透水･不透気である．

Affusion

90cm

180cm

180cm

写真 1. 試験装置とモデル地盤の鳥瞰図（Kitamura et al., 2007）
シラスの物理特性と状態変数e の初期値

表 3Water content in nature** (%)
Void ratio**(p’=10 kPa)
Density of soil particle** (g/cm3)
e (OCR)*

No.2

No.7

No.1

No.6

No.14
No.13

No.12
No.11

1681 nodes and 1600 4-nodeisoparametric elements

10 cm

No.9

10 cm

No.15

10 cm

20cm

No.10

10 cm

20cm

No.5

10 cm

20cm

No.8

Case 3
23.1
1.57
2.45
0.04 (1.5)

25 cm

Range of waterinjection (Case 2)

80 cm

Range of water injection (Case 3)
60 cm

No.3

Case 2
23.3
1.47
2.40
0.14 (4.0)

e = ( −) ln OCR; ** Date from the work by Kitamura et al (2007)*No.4

Case 1
25.6
1.57
2.45
0.04 (1.5)

45o40 cm

Range of waterinjection (Case 1)

140cm(a) Positions of measuring units and water injection(b) FEM mesh
図 12. 測定ユニット, 注水位置, FEM メッシュのレイアウト

シミュレーションで使用したシラスの材料パラメータは，北村らの研究(2007)に基づいて 0.2〜3.0MPa の拘束応力下でのシラスの力学性質に着目して推定した．表 4 に示しているパラメータは，

実験結果に基づいて推定されたものである．水分特性曲線のパラメータを表 5 に示す．
表 4 シラスの材料パラメータ

Compression index λSwelling index κ
Critical state parameter ΜVoid ratio N (p’=10 kPa on N.C.L.)
Poisson’s ratio νParameter of overconsolidation a
Parameter of suction b
Parameter of overconsolidation β
Void ratio Nr (p’=10 kPa on N.C.L.S.)

0.105
0.007
0.47
1.59
0.30
5.00
0.50
1.00
1.61

表 5 シラスの水分特性曲線に含まれるパラメータ

Saturated degrees of saturation Srs

Case 1
0.87

Case 2
0.89

Case 3
0.95

0.20
12.0
0.07
90.0
0.30
0.60
30.0

0.25

0.20

15.0
0.10
20.0

12.0
0.17
50.0

0.30
0.40
50.0

0.33
0.18
5.0

Residual degrees of saturation Srr
Parameter corresponding to drying AEV (kPa) Sd
Parameter corresponding to wetting AEV (kPa) Sw
Initial stiffness of scanning curve (kPa) k spe
Parameter of shape function c1
Parameter of shape function c2
Parameter of shape function c3

No.1
No.4
No.2
No.5
No.3 experiment

PWP (kPa)

0-2-4-6-8-10-12

No.1
No.4
No.2
No.5
No.3 calculation

050

100

150

200

250

No.6
No.9
No.7
No.10
No.8 calculation

0-2-4-6-8-10-12300

050

Elapsed Time (min)PWP (kPa)

PWP (kPa)

数値計算では，まず重力場の初期応力場を計算する．次に，第 2 段階で所定の水頭増分を与える．図 13 は，3 つのケースにおいて個々の測定位置での間隙水圧の実験と計算結果を比較したものを示している．これらの図から，シミュレーションは間隙水圧の発達を全体的によく説明できることが示された．しかし，斜面に沿った各ポイントでは，計算と実験の間にかなり大きな違いがあった．これは，鉛直方向と水平方向で透水係数が同様と仮定し，また初期サクションと飽和度がモデ
ル地盤内で均一と仮定していることが原因と考えられる．

150

200

Elapsed Time (min)
No.11
No.14
No.12
No.15
No.13 calculation

0-2-4-6-8-10-120100

No.6
No.9
No.7
No.10
No.8 experiment

50100

150

No.11
No.14
No.12
No.15
No.13 experiment

200

Elapsed Time (min)2503002503000-2-4-6-8-10-12No.4
No.5
calculation

No.1
No.4
No.2
No.5No.3 experiment

No.1
No.2
No.30
PWP (kPa)

PWP (kPa)(a) Case 1
No.1
No.2
No.34080-2-4-6-850100

1502002500300

20

Elapsed Time (min)
No.6
No.7
No.80-2-4-6-8-10-12

No.9
No.10
calculation

60100120

Elapsed Time (min)

No.6
No.9
No.7
No.10
No.8 experiment

No.6
No.7
No.80
PWP (kPa)

PWP (kPa)

No.1
No.4
No.2
No.5
No.3 experiment

-100

No.9
No.10
calculation

No.6
No.9
No.7
No.10
No.8 experiment4080-2-4-6-8-10050100150200

2500300

20

No.11
No.14
No.12
No.15
No.13 calculation

No.11
No.14
No.12
No.15
No.13 experiment

050100150200

Elapsed Time (min)

250

No.11
No.14
No.12
No.15
No.13 calculation

0PWP (kPa)

0-2-4-6-8-10-1260

100120Elapsed Time (min)Elapsed Time (min)PWP (kPa)

No.4
No.5calculation

300

No.11
No.14
No.12
No.15
No.13 experiment-2-4-6-8-10020

406080100120Elapsed Time (min)(b) Case 2(c) Case 3
図 13. 各測定位置における負の間隙水圧の実験結果と計算結果の比較（試験結果：Kitamura et al., 2007）

Case 2

Case 1Case 3

T=200 min

T=200 min

T=80 min

T=300 min

T=300 min

T=120 min

図 14. 異なる時刻における塑性偏差ひずみ分布の計算結果

図 14 は，3 つのケースについて，塑性せん断ひずみ（偏差塑性ひずみテンソルの第 2 不変量）分布の経時変化を示している．すべての場合において，最初に斜面のつま先で発達し，徐々に上面に向かって後方に伝播していくことが示されており，実験と全く同様である．また，実験の最終段階において，せん断帯はケース 1 とケース 3 ではっきりと発達したが，ケース 2 では顕著に現れなかった．いずれも実験と整合を取れている．
3.

得られた成果本研究では，有限変形アルゴリズムにおける THMA 連成問題の統一場方程式を用いた有限要素法を，最近提案された不飽和/飽和土の熱弾塑性構成モデルに基づいて詳細に導出した．場の方程式では，過剰間隙水圧と間隙空気圧，固相の変位，および温度が未知の変数として用いられる． THMA 連成問題の場の方程式を離散化する場合，平衡方程式とエネルギー保存方程式の空間離散化には FEMが使用され，連続方程式の空間離散化には後方有限差分スキームが用いられている．数値解法の妥当性は，実験室試験と現場試験によって検証された．その結果を以下の結論でまとめられる：
1）提案する数値解析手法では，有限変形アルゴリズムにおける THMA 連成問題は，統一場の方程式を使用している．場の方程式に取り込まれる構成モデルは，状態変数としてビショップ型の骨格応力，飽和度および温度を採用しているため，非等温条件下での飽和/不飽和土の挙動を統一的に記述することができる．
2）提案する数値解析手法の初期検証として，不飽和条件下でのベントナイト･母岩複合体の THMA挙動を調査するために，現場加熱試験（Munoz, 2006）をシミュレートした．シミュレーション結果
より，水の給排水，温度変化，過剰間隙水圧の経時変化などの加熱試験で観察された THMA の挙動を適切に説明できることが証明された．したがって，提案する数値解析手法は，高レベル放射性廃棄物の地層処分場における実際の現場の問題に適用できることはある程度言えることになった．
3）水注入または降雨による斜面崩壊のモデル試験も，等温条件下で同じ数値解析方法でシミュレートした．降雨や注水による水の浸透，飽和度の変化，間隙水圧と間隙空気圧の変化，モデル地盤の変形，斜面崩壊で形成されたせん断帯などの地盤の挙動は計算結果と実験結果を比較すれば，全体として統一された方法で再現できたと言える．飽和/不飽和土の適切な構成モデルが組み込まれているため，シミュレーションでは，モデル実験で観察されたモデル地盤のさまざまな破損パターンを十分に説明できることは特に言及する価値があると言える．また，解析に使用したシラスのすべての材料パラメータがすべてのテストケースで同じであることに注意されたい．これは，数値実験になくてはならない必用条件であり，そのいみでは研究の応用価値があると思われる．
4.

謝

辞

本報告書は「土･水（飽和･不飽和）･熱連成解析手法の研究」というテーマで東電設計株式会社

からの委託研究を実施した成果をまとめたものである．研究支援金に謝意を表する．
参 考 文 献Kitamura R, Sako K, Kato S, Mizushima T, Imanishi, H. (2007): Soil tank test on seepage and failure
behaviors of Shirasu slope during rainfall, Japanese Geotechnical Journal , 2(3): 149-168 (inJapanese).
Munoz J. (2006): Thermo-hydro-mechanical analysis of soft rock, application to a large scale heating
test and large scale ventilation test, PhD thesis, Catalonia, Spain: Polytechnic University ofCatalonia.
Zhang, S. and Zhang, F. (2009): A thermo-elasto-viscoplastic model for soft sedimentary rock. Soilsand Foundations , 49(4), 583-595.
Zhang, F. and Ikariya. T. (2011): A new model for unsaturated soil using skeleton and degree ofsaturation as state variables, Soil and Foundations , Vol. 51(1), 67-81.
Zhang, S., Leng, W.M., Zhang, F. and Xiong, Y. L. (2012): A simple thermo -elastoplastic model forgeomaterials. Plasticity, 34, 93-113.
Xiong Y. L., Bao X. H., Ye B., Zhang F. (2014): Soil-water-air fully coupling finite element analysisof slope failure in unsaturated ground. Soils and Foundations , 54(3), 377-395.