博士論文

農業用水路トンネルの模型実験に対する
数値解析の適用性に関する研究

2024 年 3 月

田本敏之

岡山大学大学院
環境生命科学研究科

目

第1章

次

序論 ······································································································ 1

1.1 研究の背景 ··································································································· 1
1.2 研究の目的 ································································································· 3
1.3 水路トンネルに関する過去の研究 ····································································· 6
1.4 本研究の基礎となる考え方 ············································································ 8
1.4.1 模型実験に関する研究 ············································································· 8
1.4.2 数値解析に関する研究 ············································································· 11
1.5 本論文の構成 ······························································································ 13

第2章

模型実験 ······························································································· 14

2.1 模型実験の概要 ····························································································· 14
2.1.1 模型実験装置 ························································································· 14
2.1.2 計測方法の詳細 ······················································································ 16
2.1.3 模型の物性値 ························································································· 17
2.1.4 実験方法 ······························································································· 19
2.2 模型実験結果 ································································································ 20
2.2.1 模型実験結果（空洞なし） ······································································· 20
2.2.2 模型実験結果（空洞あり：空洞範囲 45°） ················································· 24
2.2.3 模型実験結果（空洞あり：空洞範囲 90°） ················································· 27
2.2.4 実験結果の比較 ······················································································ 30

第3章

数値解析 ···················································································· 33

3.1 剛体ばねモデル ····························································································· 33
3.1.1 平面要素のモデル化·················································································· 33
3.1.2 本研究における構成モデル ········································································ 37
3.1.3 剛体ばねモデルによる解析結果 ·································································· 42
3.1.4 剛体ばねモデルによる解析結果の比較 ························································· 54
3.2 個別要素法 ··································································································· 56
3.2.1 個別要素法の基礎理論··············································································· 56
3.2.2 衝突検出································································································· 58
3.2.3 ひずみ ···································································································· 61
3.2.4 応力 ······································································································· 64
3.2.5 運動結合································································································· 68
3.2.6 本研究における構成モデル ········································································ 69
3.2.7 個別要素法による解析結果 ········································································ 72
3.2.8 個別要素法による解析結果の比較 ······························································· 75

第4章

模型実験に対する数値解析の適用性······················································ 76

4.1 模型実験に対する剛体ばねモデルの適用性 ························································· 76
4.1.1 模型実験結果と解析結果の比較（空洞なし） ················································ 76
4.1.2 模型実験結果と解析結果の比較（空洞あり：45°） ······································· 81
4.1.3 模型実験結果と解析結果の比較（空洞あり：90°） ······································· 83
4.1.4 剛体ばねモデルの適用性に関する考察 ························································· 85
4.2 模型実験に対する個別要素法の適用性 ······························································· 86
4.2.1 模型実験結果と解析結果の比較 ·································································· 86
4.2.2 個別要素法の適用性に関する考察 ······························································· 89

4.3 数値解析の適用性に関する評価 ········································································ 90
4.3.1 数値解析による模型実験の再現性の特徴 ······················································ 90
4.3.2 実務における数値解析法の適用性 ······························································· 92

第5章

結論 ···································································································· 94

参考文献 ··········································································································· 96

謝辞

第1章

序論

1.1 研究の背景
現在，全国の基幹的な農業水利施設は約 7 千箇所，支線水路を含めた農業用用排水路は約 40 万 km に及び、基幹的農業水利施設の総資産価値（再建設費ベース）で 18.5 兆円のストックが蓄積されている（Fig.1.1）．

農業用用排水路
うち基幹的水路約 40 万 km 以上(地球約 10 周分)
約 5 万 km

ダム、取水堰、用排水機
約 7 千箇所
場等

Fig.1.1 農業水利ストックの現状（出典：H31.2 農林水産省ホームページ）

このうち，水路トンネルは，全国で約 2,000km 以上整備されており，道路トンネル約
4,100km，鉄道トンネル約 2,400km，発電用水路トンネル約 4,700km に匹敵する膨大なストックを有している．しかし，これらのトンネル構造物は昭和以前に矢板工法によって建設されたものが多く，長期供用によって変状が進行し，近年，山陽新幹線福岡トンネルの覆工コンクリート落下事故（平成 11 年），広島県の送水トンネル崩落事故（平成 18 年），
中央自動車道笹子トンネルの天井板落下事故（平成 24 年）などの突発的な事故が発生している．
このような突発的な事故に対するリスク管理の観点から，水路トンネルにおいても，適時適切な診断と対策が求められている．特に，一部の水路トンネルでは，農業用水のみならず，上水，工業用水，発電用水と兼用される重要な社会インフラとなっている場合もあり，計画的な機能保全対策の実施が不可欠である．
そのため，施設の監視・診断等によるリスク管理を行いつつ施設の性能低下状況に応じた補修・補強等を計画的に行うことにより、施設の長寿命化とライフサイクルコストの低

1減を図る戦略的な保全管理の推進が必要とされているが，戦略的な保全管理を推進するためには，ストックマネジメントのプロセスの中で確認された施設の性能低下状況，要求性能，使用環境条件等に応じて，適切な時期に補修や補強等の長寿命化対策を実施し，低下した施設性能の回復又は向上を行い施設の長寿命化を図ることが重要である．農林水産省においても，令和 3 年 6 月に水路トンネルの調査から対策工法の選定，対策工法の設計，対策後のモニタリングに関する内容を網羅した「農業水利施設の補修・補強工事に関するマニュアル【水路トンネル編】」を策定し，取り組みを進めているところである．
上記のような戦略的な保全管理を行っていくためには，トンネルの荷重状態を把握することが不可欠であり，そのためにはひび割れの評価が重要となる．過去の研究において，有限要素法（FEM）を用いた覆工の安全性の評価手法を提案した事例はあるが，有限要素法は，ひび割れの評価が困難である．本研究では，実務レベルをターゲットとして，解析が比較的容易で，ひび割れを良好に再現することが可能と考えられる「剛体ばねモデル（RBSM：Rigid Body Spring Model）」について，不連続体の解析に多く用いられる「個別要素法（DEM : Distinct Element Method）」と比較することで適用性を評価する．

2

1.2 研究の目的
水路トンネルの本来性能である通水性能（水理性能、水利用性能）を維持するためには，構造性能を維持していく必要があり，構造性能の適切な評価のためには，ひび割れの評価が重要である．ひび割れの構造性能に対する影響は，ひび割れの発生位置や幅，走行から，
塑性圧（Fig.1.2），緩み圧（Fig.1.3），偏圧（Fig.1.4）に概ね区分することが可能であるが，覆工背面で発生している事象を扱う特性から，ひび割れが発生した地圧強度やひび割れの進展仮定を把握することは困難である．

Fig.1.2 塑性圧による代表的な変状形態模式図（出典：トンネルの変状メカニズム、土木学会、平成 15 年 9 月、P.55）

3

Fig.1.3 緩み土圧による代表的な変状形態模式図（出典：トンネルの変状メカニズム、土木学会、平成 15 年 9 月、P.93）

Fig.1.4 偏圧による代表的な変状形態模式図（出典：トンネルの変状メカニズム、土木学会、平成 15 年 9 月、P.117）

4

そこで，本研究では，矢板工法で施工された農業用水路トンネルへの適用をイメージし，現場で発生するひび割れ発生過程を数値解析で再現することを目的として，以下の検討を行った．
1）矢板工法で施工された農業用水路トンネルを想定し，アーチ頂部の覆工背面に空洞を有する地山と覆工コンクリートを模擬した 1/10 スケールの模型実験によって塑性
圧によるひび割れ発生過程を再現する．
2）模型実験に対する数値解析には，実務レベルをターゲットとして，解析が比較的容易で，ひび割れを良好に再現することが可能と考えられる「剛体ばねモデル（RBSM：
Rigid Body Spring Model）」について，「個別要素法（DEM : Distinct Element Method）」と比較することで適用性を評価する．

5

1.3 水路トンネルにおける過去の研究
過去に水路トンネルに関する研究として，診断・点検手法、健全度評価手法に関する研究，変状メカニズムに関する研究，補強対策に関する研究等，多数の研究成果が報告されている．
診断・点検手法、健全度評価手法に関する研究に着目すると，佐野ら（2006）は，トン
ネル覆工の残存耐力を予測する手法の確立を目指し，鉛直圧ならびに水平圧に関する模型実験結果により，剛性変化に基づく健全度評価手法を提案した．上原（2009）は，水力
発電所無圧水路トンネルを対象に，評価の個人差を排除することを目的とした健全度評
価手法について検討した．藤原ら（2019）は，通水状態での点検手法の確立を目的として，
フロート式画像撮影装置を開発し，現地実証試験による性能確認を行った．
変状メカニズムに関する研究に着目すると，仲野（1974）は，トンネル支保工に作用す
る地圧のうち，粘土性岩による膨張性押し出し地圧に着目し，地圧の発生状況を計測するとともに，発生メカニズムについて考察した．
補強対策に関する研究として，長束ら（2014）は，発泡ウレタンによる裏込注入工法の有効性について検討した．櫻井ら（2019）は，炭素繊維ストランドシートとセラミック混合型エポキシ樹脂モルタルを用いた内面補強工法を開発し，実物大模型を用いた模型実験によって有効性を確認した．
実験に着目すると，森ら（2015，2021）は，円形および馬蹄形の農業用水路トンネルと地山を模擬した室内実験を行い，覆工及び地山に発生するひび割れの性状について検討した．長束ら（2013），朝倉ら（1994a, 1994b），野城ら（2015）は，エアバッグや載荷用
ボルトを用いて覆工に直接載荷する実験を行っている．横山ら（2001）は，馬蹄形覆工を対象に，二種類の載荷方向による模型実験と数値解析を実施している．また，石井ら（2014）
は，現場載荷試験と有限要素法（FEM）の解析結果に基づき，安全性の評価手法の提案を行っている． Si et al.（2021）は，トンネル覆工を設定せずに 3 次元載荷実験を行い，ほろ形トンネルの剥離による破壊の検討を行っている．
また数値解析に焦点を当てると，Wang et al.（2014）が FEM により背面空洞の位置の影響を検討し，保田ら（2016）および Yasuda et al.（2019）は三次元弾性論に基づき選点法とサブストラクチャー法を組み合わせ，背面空洞の検討を行っている．以上のように地山あるいは覆工のみを考慮している研究も散見されるが，覆工を含めたトンネルの安定を6

検討するには，地山と覆工の相互作用を考慮し，ひび割れの発生を評価するような不連続体を対象とした数値解析手法が必要となる．

7

1.4 本研究の基礎となる考え方
本研究では，塑性圧が発生した水路トンネルを模擬した模型実験を実施し，数値解析による模型実験を再現する．ここでは，本研究の基礎となる考え方として，前項に示した過去の研究のうち，水路トンネルの模型実験に関する研究内容，数値解析（剛体ばねモデル，
個別要素法）に関する研究内容について整理する．

1.4.1 模型実験に関する研究
森ら（2015）は，水路トンネルの覆工
背面の空洞分布に着目し，円形トンネルと地山を模擬し，鉛直荷重のみを作用させた室内実験（Fig.1.5）を行い，周辺地盤と覆工の挙動を，デジタル画像相関法（Digital Image Correlation Method）により計測している．この結果から，円形トFig.1.5 模型実験の概要（円形）ンネルにおける鉛直荷重作用時の背面空洞有り，無しの条件における変形の特徴（Fig.1.6）を明らかにするとともに，デジタル画像相関法により，模擬地盤のひび割れ発生後の挙動や模擬地盤の変形の向きを定量的に計測することが可能であることを明らかにした．本研究では，この実験手法を参考に，模型実験の条件設定を行った．

8

（a）背面空洞なし（b）背面空洞有り（90°）

Fig.1.6 周辺地盤の変形挙動（円形トンネルの実験結果）

9さらに，森ら（2021）は，水路トンネルの覆工背面の空洞分布に着目し，2R 標準馬蹄形
トンネルと地山を模擬し，鉛直荷重のみを作用させた室内実験（Fig.1.7）を行い，周辺地盤と覆工の挙動を計測している．この結果から，
2R 標準馬蹄形トンネルにおける鉛直荷重作

用時の背面空洞有り，無しの条件における変形の特徴（Fig.1.8）を明らかにした．

Fig.1.7 模型実験の概要（2R 標準馬蹄形）

Fig.1.8 周辺地盤の変形挙動（2R 標準馬蹄形における実験結果）

101.4.2 数値解析に関する研究（1）剛体ばねモデル
不連続体の解析には，DEM，マニフォールド法（MM），不連続変形法（DDA）などがあり，動的な変形問題に適用されることも多い．また，FEM を用いてひび割れやき裂などを解析する場合は，メッシュ分割と独立に不連続変形の幾何学性を表現するといった方法がとられる．剛体ばねモデル（RBSM : Rigid Body Spring Model，以後 RBSM と称す）は，要素自身を剛体と仮定し，要素間にばねを設けて解析を行う手法である．解析結果がメッシュ分割に依存するといった特徴があるが，亀裂やひび割れの進展を良好にシミュレートする長所がある．
これまでに，竹内ら（1981）は，トンネル掘削時の坑口付近の斜面安定解析に RBSM を適用し，すべり崩壊に対する RBSM の適用性を示した．城野ら（1994）は，地山と支保工をモデル化した RBSM と FEM の比較によって，RBSM のトンネル問題への適用の可能性と課題を示した．近藤（1996）は，RBSM を用いて間隙水圧を考慮した斜面安定計算を実施し，従来実施されている極限平衡法との比較により，斜面安定問題に対する RBSMの有効性を示した．濱崎ら（2006）が斜面安定問題に三次元 RBSM を適用し，従来の計算法と同様に安全率を算出できることを確認している．
近年では，山本ら（2010，2015）がコンクリートの破壊解析， Yao et al.（2022）により岩石に対する検討が試されている．
これまでの研究において，トンネル施工時の坑口付近のすべり崩壊や，地山と支保工をモデル化した施工時の解析，斜面安定解析，コンクリートや岩石材料に対する解析が行われているが，背面空洞を有する覆工コンクリートと地山を想定した検討が行われた研究は無く，本研究における新しい着眼点である．

11

（2）個別要素法
個別要素法（DEM : Distinct Element Method）は定量的な評価が難しい側面がある一方で，局所的な変形や破壊の再現に優れているという長所があり，例えば，蒋ら（2001），
Funatsu et al.（2008），Boon et al.（2015）は，二次元の DEM を用いた大規模地下空洞およびトンネルの解析を行っている．また，Khanal et al.（2005）や Beckmann et al.（2012）および Riera et al.（2016）は二次元および三次元のコンクリートの破壊解析に DEM を適用している．
本研究で扱うような覆工や地山のひび割れの進展を模擬しようとした場合，三次元のDEM を採用し，ひび割れの進展や内空変位量の把握といった変状の定性的な適用性を明確にすることが非常に重要となる．
DEM のオープンソースとして公開されているソフトウェアとして YADE があるが，本研究においても YADE を用いて解析を行う．
YADE を用いた研究報告には，
Šmilauer（2010），
Kozicki and Donzé （2008），Boon et al. （2015）および Nitka and Tejchman （2015）などがある．
本研究ではほろ形のトンネルを対象とし，地山とトンネル覆工に対する模擬的な載荷実験を実施し，覆工の背面空洞の影響を検討する．横方向（トンネルの水平方向）からの塑性圧のみを模擬するため，荷重をトンネル側部からの載荷に限定し，地山と覆工の相互作用を検討する．DEM による解析は，模型実験と同寸法の三次元 DEM モデルを用いて，
実施する．

12

1.5 本論文の構成
本研究では，矢板工法で施工された農業用水路トンネルへの適用をイメージし，現場で発生するひび割れ発生過程を数値解析で再現することを目的として，以下の検討を行った．
1）矢板工法で施工された農業用水路トンネルを想定し，アーチ頂部の覆工背面に空洞を有する地山と覆工コンクリートを模擬した 1/10 スケールの模型実験によって塑性
圧によるひび割れ発生過程を再現する．
2）模型実験に対する数値解析には，実務レベルをターゲットとして，解析が比較的容易で，ひび割れを良好に再現することが可能と考えられる「剛体ばねモデル（RBSM：
Rigid Body Spring Model）」について，「個別要素法（DEM : Distinct Element Method）」と比較することで適用性を評価する．

本論文は，これらの研究内容を第 1 章の序論から第 5 章の結論の 5 章構成でとりまとめたものである．各章の記載内容は以下のとおりである．
第 2 章は，模型実験の方法，模型寸法，実験に用いた材料の物性値等を整理し，模型実
験の結果について整理した．
第 3 章は，剛体ばねモデル（RBSM）と個別要素法（DEM）の基本的な考え方について
概説し，本研究で用いた構成モデルについて整理したうえで，各解析手法における解析結
果を取りまとめた．
第 4 章では，模型実験結果と各解析手法における解析結果を比較し，各解析手法のひび
割れ発生過程の再現性を検証し，適用性を評価した．
第 5 章では，本研究の内容を総括し，結論を述べた．

13

第2章

模型実験

本章では、本研究で実施した模型実験の概要と実験結果について整理する．

2.1 模型実験の概要
2.1.1 模型実験装置
模型実験は，1/10 スケールの農業用水路トンネルを模擬し，960mm×960mm×400mm の地山を模した供試体の中央に，巻厚 20mm，幅員 240mm，高さ 240mm のほろ形のトンネルを配置する．Fig. 2.1 に模型実験装置の概略図，Fig. 2.2 に模型の寸法と計測位置を示す．
油圧シリンダー

荷重計

地山（低強度モルタル）

覆工（高強度モルタル）

外枠：H 形鋼
背面：鋼板

載荷板

Fig. 2.1 模型実験装置

ポンプ

ここで，地山は背面の鋼板と H 形鋼で作製された実験装置の内側に設置し，地山と背面の鋼板の間には PTFE シートを挿入して，載荷時の摩擦を軽減する．ただし，地山の上部・下部と実験装置との間には 2cm 程度の隙間があり，地山の上下の変形は拘束されていない．

14

（単位：mm）
奥行：400，覆工の巻厚：20
載荷板
地山

240

340

960

280

覆工

240

荷重測定

内空変位量測定
340

動画撮影
載荷板
340

280

340

変位測定

960

Fig. 2.2 模型の寸法と計測位置

Fig. 2.3 に示す背面空洞範囲θは，0°，45°，90°の 3 ケースを設定した．

Fig. 2.3 背面空洞の範囲

15

2.1.2 計測方法の詳細
前述の模型実験装置を用い，地山に作用する塑性圧を表すため右側面から載荷板を介して油圧シリンダーにより荷重を作用させる．油圧シリンダーと載荷板の間には荷重計を据え付け，載荷板変位（接触式変位計による）と覆工の縦／横方向内空変位量（非接触式変位計による）とともに測定する．ここで，荷重計，載荷板変位，内空変位量は模擬的な塑性圧，
模擬地山の変形，覆工の変形を評価するために設置している．また，実験中のひび割れの進展を確認するため，デジタルビデオカメラを用いて，動画撮影を行っている．Fig. 2.4 に計測位置の詳細を示す．

Fig. 2.4 計測位置の詳細

接触式変位計は地山表面（図では上面）から 50mm，図の手前から 10mm の位置に設置しており，非接触式変位計は覆工の片側に固定したレーザー照射板，その反対側に非接触式変位計を固定して計測を行っている．なおこの図では，縦方向と横方向の内空変位量の測定についてそれぞれ示しているが，実際は両者を同時に設置して計測を実施している．また，荷重計は Fig. 2.1 にも示しているように載荷板と油圧シリンダーとの間に配置しており，表面（図では上面）から荷重計の中心までが 200mm，3 個の荷重計の間隔（中心距離）が 250mmになるように配置している．

16

2.1.3 模型の物性値
地山および覆工は低強度モルタル（w/c=63%，単位セメント量 351kg/m3，単位水量 222
kg/m3，単位細骨材量 1,044 kg/m3，細骨材密度 2.61g/cm3，細骨材粗粒率 2.85）と高強度モルタル（w/c=45%，単位セメント量 751kg/m3，単位水量 338 kg/m3，単位細骨材量 1,095 kg/m3，
細骨材密度 2.60 g/cm3，細骨材粗粒率 2.59）を用い，早強ポルトランドセメントにより作製した．養生期間は 7 日としているが，実験装置に地山を設置する業者の都合や天候などの理由により，7 日から数日後に実験を実施することがあった．また，覆工を型枠中央に設置した後にモルタルを打設して地山を作製しているため，両者は分離せず，固着した状態で実験を行う．ここで，覆工は，Fig.2.5 に示す型枠を用いて作製している．覆工上部の円形部は鋼製と EPS 製（養生テープを表面に貼付）の型枠，下部の矩形部はコンクリート型枠用合板の型枠を用いており，いずれも表面は滑らかである．そのため，完成した覆工も滑らかな表面となっている．なお，覆工のモルタルは，作製後 2～3 日が経過してから脱型し，地山供試体を作製している．覆工の打設と地山の打設までの時間には，2～10 日程度の時間の差がある．

Fig. 2.5 覆工の型枠
Table 2.1 に地山と覆工に用いたモルタルの圧縮試験と割裂試験の結果を示す．ここでは，

材齢 7 日の 6 つの供試体（3 回の打設に対し，2 つの供試体を作製）の試験結果の平均値を示しており，後述する解析でもこの値を用いている．なお，割線静弾性係数は，圧縮試験での最大応力の 1/3 の応力を示す点とひずみが 50×10-6 の点とを結ぶ線分の勾配から算出し，ポアソン比は供試体の載荷縦方向および円周方向に張り付けたひずみゲージにより測定した．

17

Table 2.1 圧縮試験および割裂試験の結果

地山

覆工

割線静弾性係数1.21 kN/mm2

5.00 kN/mm2

圧縮強度

8.56N/ mm2

23.2 N/ mm2

引張強度1.37 N/mm2

3.32 N/mm2

1.13×10-3

5.54×10-3

0.20

0.20

圧縮限界ひずみ
ポアソン比

182.1.4 実験方法
実験は，油圧シリンダーは最大載荷荷重 300kN のものを 3 個用いて，載荷変位が 0.1mm/s
以下になるよう油圧シリンダーに荷重を加え，荷重と載荷板変位が一定になるまで載荷を行った．本実験では，変位速度の制御を厳密には行っていないため，載荷荷重が油圧シリン
ダーの最大出力に近づくとその値が一定になり，載荷板変位の値もほぼ変わらなくなる．そのため，この時点で実験を終了としている． 実験状況写真を Fig. 2.6 に示す．

Fig. 2.6 実験状況写真19

2.2 模型実験結果
模型実験は，空洞なし，空洞あり（空洞範囲 45°）
，空洞あり（空洞範囲 90°）の 3 ケースについて実施した．各ケースの模型実験結果について，以下に整理する．

2.2.1 模型実験結果（空洞なし）（1）時間と荷重の関係
空洞無しの実験は，3 回実施した．Fig. 2.7 に 3 回の実験における時間と荷重の関係を示す，ここで，3 回の実験名称を a-1，a-2，a-3 として整理した．後述する数値解析の入力条件
として，模型実験によって得られた時間と荷重の関係を使用する．

荷重 [N][×105]

10

5
a-1

a-2

a-3
0
0

200

400
時間 [s]

600

Fig. 2.7 模型実験における時間と荷重の関係

20

800

（2）荷重と載荷板変位の関係
Fig. 2.8（a）～（c）に 3 回の実験における荷重と載荷板変位の関係を示す．載荷板変位

は，荷重の増加に伴って，ほぼ比例的に増加し，最大の載荷板変位は，2～2.5mm 程度であ載荷板変位 [mm]

った．

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0

Fig. 2.8（a）

載荷板変位 [mm]

←ひび割れ発生

a-1

5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0

5
荷重 [N] [×105]10

模型実験における荷重と載荷板変位の関係（a-1）

a-2←ひび割れ発生

0Fig. 2.8（b）

5
荷重 [N] [×105]10

模型実験における荷重と載荷板変位の関係（a-2）

21載荷板変位 [mm]

5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0

a-3

←ひび割れ発生

0Fig. 2.8（c）

5
荷重 [N] [×105]10

模型実験における荷重と載荷板変位の関係（a-3）

22

（3）載荷板変位と内空変位の関係
Fig. 2.9 に載荷板変位と内空変位の関係を示す．縦軸が縦方向と横方向の内空変位量で，それぞれ正と負の値を示しており，横軸が載荷板変位を表している．側面から荷重を作用させることで，覆工が縦に伸張するように変形し，それにより縦方向の内空変位量は増加している．また，横方向には押しつぶされるように変形するため，内空変位量が減少していることになる． ひび割れ発生時の載荷板変位は，1.4～1.5mm 程度であり，最大の内空変位は0.5mm 程度であった．なお，目視（動画）によるひび割れ判定は，微細なひび割れや模型背面に発生したひび割れを確認することが困難であるため，ここでは，内空変位が急激に増大するタイミング（載荷板変位）をひび割れ発生と定義する．

ひび割れの発生
3

内空変位 [mm]23.E+06

縦方向

背面空洞なし3.E+06a-1

1

2.E+060

2.E+06
1.E+06-1a-2-2

5.E+05

a-3

横方向

-3
0

1

2

3

4

0.E+00載荷板変位 [mm]
a-1

a-2

a-3

Fig. 2.9 模型実験における載荷板変位と内空変位の関係

23

2.2.2 模型実験結果（空洞あり：空洞範囲 45°）
（1）時間と荷重の関係
空洞あり（空洞範囲 45°）の実験は，2 回実施した．Fig. 2.10 に 2 回の実験における時間と荷重の関係を示す，ここで，2 回の実験名称を b-1，b-2 として整理した．実験の結果，荷重は経過時間に応じてほぼ線形的に増加し，7.0～7.5×105 N 付近でほぼ一定値を示した．

荷重 [N][×105]

10

5

b-1

b-2

0
0

200

400
時間 [s]

Fig. 2.10 模型実験における時間と荷重の関係

24

600

（2）荷重と載荷板変位の関係
Fig. 2.11（a），（b）に 2 回の実験における荷重と載荷板変位の関係を示す．載荷板変位は，

載荷板変位 [mm]荷重の増加に伴って，比例的に増加し，最大の載荷板変位は，2.5～3.0mm 程度であった．

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0

Fig. 2.11（a）

載荷板変位 [mm]

←ひび割れ発生

b-1

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

5
荷重 [N] [×105]

模型実験における荷重と載荷板変位の関係（b-1）←ひび割れ発生b-2

0Fig. 2.11（b）

105
荷重 [N] [×105]10

模型実験における荷重と載荷板変位の関係（b-2）

25

（3）載荷板変位と内空変位の関係
Fig. 2.12 に載荷板変位と内空変位の関係を示す．縦軸が縦方向と横方向の内空変位量で，それぞれ正と負の値を示しており，横軸が載荷板変位を表している．側面から荷重を作用させることで，覆工が縦に伸張するように変形し，それにより縦方向の内空変位量は増加している．また，横方向には押しつぶされるように変形するため，内空変位量が減少していることになる．ひび割れ発生時の載荷板変位は，1.5 ～2.0mm 程度であり，最大の内空変位は1.0mm 程度であった．なお，目視（動画）によるひび割れ判定は，微細なひび割れや模型背面に発生したひび割れを確認することが困難であるため，ここでは，内空変位が急激に増大するタイミング（載荷板変位）をひび割れ発生と定義する．

ひび割れの発生

内空変位 [mm]3縦方向3.00E+06

背面空洞45°2

2.50E+06

1

2.00E+06

0

1.50E+06
1.00E+06

-1
b-1

b-2

-2

5.00E+05

横方向
0.00E+00

-3
0

1

2

3

4

載荷板変位 [mm]
b-1

b-2

Fig. 2.12 模型実験における載荷板変位と内空変位の関係

26

2.2.3 模型実験結果（空洞あり：空洞範囲 90°）（1）時間と荷重の関係
空洞あり（空洞範囲 90°）の実験は，2 回実施した．Fig. 2.13 に 2 回の実験における時間と荷重の関係を示す，ここで，2 回の実験名称を c-1，c-2 として整理した．実験の結果，荷重は経過時間に応じてほぼ線形的に増加し，8.0～8.5×105 N 付近でほぼ一定値を示した．

荷重 [N][×105]

10

5

c-1

c-2

0
0

200

400
時間 [s]

Fig. 2.13 模型実験における時間と荷重の関係

27

600

（2）荷重と載荷板変位の関係
Fig. 2.14（a），（b）に 2 回の実験における荷重と載荷板変位の関係を示す．載荷板変位は，荷重の増加に伴って，初期は急激に増加し，徐々に増加は鈍化した．最大の載荷板変位は，
3.5～4.0mm 程度であり，空洞なし及び空洞有り（45°）と比較して，載荷板変位は大きく

載荷板変位 [mm]

なった．

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0

Fig. 2.14（a）

載荷板変位 [mm]

←ひび割れ発生

c-1

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

5
荷重 [N] [×105]

模型実験における荷重と載荷板変位の関係（c-1）←ひび割れ発生c-2

0Fig. 2.14（b）

105
荷重 [N] [×105]10

模型実験における荷重と載荷板変位の関係（c-2）

28

（3）載荷板変位と内空変位の関係
Fig. 2.15 に載荷板変位と内空変位の関係を示す．縦軸が縦方向と横方向の内空変位量で，それぞれ正と負の値を示しており，横軸が載荷板変位を表している．側面から荷重を作用させることで，覆工が縦に伸張するように変形し，それにより縦方向の内空変位量は増加している．また，横方向には押しつぶされるように変形するため，内空変位量が減少していることになる．ひび割れ発生時の載荷板変位は，2.8～3.3mm 程度であった．最大の内空変位は1.0～2.0mm 程度であり，空洞無し，空洞あり（45°）と比較して，内空変位は大きくなった．

ひび割れの発生3.00E+06

3

背面空洞90°内空変位 [mm]2c-12.50E+06

縦方向

12.00E+06

0

1.50E+06-1

1.00E+06

c-2-2

5.00E+05

横方向
0.00E+00

-3
0

1

2

3

4

載荷板変位 [mm]
c-1

c-2

Fig. 2.15 模型実験における載荷板変位と内空変位の関係

29

2.2.4 実験結果の比較
ケース毎の実験結果について，時間と荷重，荷重と載荷板変位，載荷板変位と内空変位の関係に着目して整理した。

（1）時間と荷重の関係
Fig. 2.16 に時間と荷重の関係を示す．比較の結果，a-1 が特異な結果であるものの，空洞範囲が小さい方が，荷重増加の傾きがやや大きくなることが確認された．これは，空洞範囲が小さいほど，模型が変形しにくく，反力が大きくなるためであると考えられる．また，最終的な最大荷重は，実験毎のばらつきはあるものの，ケース毎の特徴は見いだせなかった．

10

a-2

荷重 [N][×105]

a-3

b-1
c-1

b-2

5

c-2
a-1

0
0

200

400
時間 [s]

600

Fig. 2.16 模型実験における時間と荷重の関係

30

800

（2）荷重と載荷板変位の関係
Fig. 2.17 に荷重と載荷板変位の関係を示す．比較の結果，空洞範囲が大きいほど，荷重に対する載荷板変位が大きくなった．これは，側面から荷重を作用させることで，覆工が縦に伸張するように，横方向には押しつぶされるように変形するため，空洞範囲が大きいほど，この変形が顕著になる．この変形に伴い，空洞範囲が大きくなるに伴って，載荷板

載荷板変位 [mm]

変位も大きくなることが確認された．

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

←ひび割れ発生

c-1

c-2

a-3
b-2

b-1

a-1
a-2

0

5
荷重 [N] [×105]

Fig. 2.17 模型実験における荷重と載荷板変位の関係

31

10

（3）載荷板変位と内空変位の関係
Fig. 2.18 に載荷板変位と内空変位の関係を示す．比較の結果，空洞なし，空洞あり（45°），空洞あり（90°）の順に内空変位が急激に増大（ひび割れが発生）する載荷板
変位が大きくなる．また、空洞なし，空洞あり（45°），空洞あり（90°）の順に最大の内空変位が大きくなる．

3

縦方向

c-1

b-1

2

内空変位 [mm]

a-1

a-2

1

a-3

b-2

c-2

0

c-2

b-2

c-1

-1

a-1

a-2

-2

a-3
b-1

横方向
-3
0

1

2

3

4

載荷板変位 [mm]
背面空洞なし

背面空洞45°背面空洞90°Fig. 2.18 模型実験における載荷板変位と内空変位の関係

32

第3章

数値解析

本章では、本研究で用いる「剛体ばねモデル（RBSM：Rigid Bodies-Spring Model）」及
び「個別要素法」の考え方と解析結果について整理する．

3.1 剛体ばねモデル
剛体ばねモデルについては，「鉄筋コンクリート構造の離散化極限解析」（竹内則雄ら,2005）に基づきモデル化した．以下に基本理論を整理する．ここでは，概要を整理したため，詳細は参考文献を参照されたい．

3.1.1 平面要素のモデル化（1）剛体ばねモデルの基礎理論
RBSM では，変位場として剛体変位場を仮定（Fig.3.1）しているため，要素形状は任意多角形を用いることができる．微小変形を仮定した場合，任意多角形内のある点（ xG , yG ）における x, y 方向の平行変位と剛体回転角を（ u , v,）とすれば，点 P（ x, y ）における x 方向， y 方向の移動量（ U , V ）は式 3.1 で表すことができる．

U

u

y

yG

V

v

x xG

(式 3.1)vu（ xG , yG ）Fig.3.1 剛体変位場

33

P（ x , y ）また，変形前後の 2 つの剛三角形要素の関係を示すと Fig.3. 2 のようになるが，このときの任意点 P における 2 要素間の相対変位量 δ は、式 3.2 により求まる．
2

Mi U i（式 3.2）

i 1

t

n

ここに、

, s , M1

1

0

0

1

,M2

1 0
0 1

：2 要素間の相対変位ベクトル

M i ：方向余弦行列

U i ：全体座標系における平行変位（a）変形前（b）変形後

Fig.3. 2 変形前後の要素関係

式 3.2 に対し，式 3.1 を要素ごとの物理量を表すものとして一括表記し，要素境界辺にそった座標系に変換するための座標変換行列 R を用いて変換して代入すれば 1 要素間の相対変位ベクトル は要素重心の剛体変位 ui を用いて式 3.3 のように算出できる。
2

M i R Qi ui
i 1
2

Bi ui

Bu

（式 3.3）i 1

Bi

M i R Qi

ここで，B は，隣接する 2 要素間の重心に関する剛体変位を関連付けるための行列である．

34

（2）ばねのモデル化
ばねの力学では，力と変位の関係は P
に相対変位

n

, s

t

k

で表現されるが，RBSM においても，同様と，単位面積あたりの表面力

n

, s

t

との間に式 3.4 の関係が成立しているものと仮定する．

D
ここで，

n

D

kn

0

0

ks

（式 3.4）, s は，垂直及びせん断方向における単位面積あたりの表面力（ベクトル量）

であり，有限要素法のテンソル量とは異なる性質を持つ．ここでは，D をばね行列とよぶ．

式 3.4 におけるばね定数 kn , ks を決定するために，単位長さあたりの相対変位を式 3.5 のとおり仮定する．

h

ここで，

n

s

h

h1 h2

（式 3.5）であり， n は法線方向の単位長さあたりの相対変位，

s はせん

断方向の単位長さあたりの相対変位を表している．また， h1 , h2 は Fig.3.3 に示す三角形要
素の重心から要素境界辺上に下ろした垂線の長さである．

Fig.3.3 垂線の定義

35一方，平面応力状態を仮定した場合，1 軸状態の応力－ひずみ関係は，弾性係数 E，ポ
アソン比νとして式 3.6 のように与えられる．

E
n

2

1

n

,

E
s

s

1

（式 3.6）

これらの関係から，平面応力状態のばね定数が式 3.7 のように決定される．

kn

E
1

2

h

E

, ks

1

h

（式 3.7）

要素境界辺上に分布しているばねに蓄えられるエネルギーを評価することで，剛性行列を導くことができる．詳細は参考文献を参照されたい．

36

3.1.2 本研究における構成モデル（1）ばねのモデル化
本論文では，地山のばねの構成モデルを Fig.3.4，覆工のばねの構成モデルの Fig.3.5 のとおり仮定した．また，Table 3.1 に実験材料の圧縮試験及び割裂試験の結果，Table 3.2 に解析パラメータを示す．

（a）垂直ばね（b）せん断ばね

Fig.3.4 地山のばねの構成モデル（a）垂直ばね（b）せん断ばね

Fig.3.5 覆工のばねの構成モデル

Fig₋3.4 (a)および Fig₋3.5 (a)では， Table3.1 の圧縮強度 Fc と限界圧縮ひずみεcu を圧縮側のピークとして，線形で応力とひずみの関係を仮定している．竹内ら（2005）は，一軸圧縮試験の結果に基づき，ひずみが圧縮限界ひずみに達した後は，圧縮限界ひずみの 2 倍に達するまで応力を解放するものとし，その後，0.2Fc の応力を保持するものと仮定している．
本論文では，この圧縮試験の結果に基づく設定を参考に構成モデルを決定した．また，圧縮37

側の傾きを引張側に延長し，Table 3.1 に示す引張強度 Ft に達した後にひび割れが発生したものとする．ここで Ft に対応するひずみをひび割れひずみεcr とし，ひび割れが発生する位置を示している．本論文ではεcr を超えた後は，εcr0 まで応力解放を行うとし，線形で仮定した．なおεcr0 は，同実験結果を用いている柴田ら（2020）および Šmilauer（2010）の研究で得た引張部を線形で表すことから算出しており，Fig. 3.4(a)および Fig.3.5(a)でそれぞれεcr0=-3.62×10-4，εcr0=-1.33×10-3 である．
Fig. 3.4 (b)および Fig. 3.5 (b)では竹内ら（2005）の設定と同様に，せん断ひずみγu=0.004

でせん断破壊が発生すると仮定している．ここで，本論文の解析結果に対するせん断ひずみ
γu の感度は低く，0.001～0.005 の範囲で解析を実施しても結果に差は生じなかったため，この値を設定した．また，せん断弾性係数 G は，弾性係数 E とポアソン比νを用いて求めている．

Table 3.1 圧縮試験および割裂試験の結果
地山

覆工

割線静弾性係数1.21 kN/mm2

5.00 kN/mm2

圧縮強度

8.56N/ mm2

23.2 N/ mm2

引張強度1.37 N/mm2

3.32 N/mm2

圧縮限界ひずみ

1.13×10-3

5.54×10-3

ポアソン比

0.20

0.20

38

Table 3.2 解析パラメータ
地山

覆工

決定方法

ポアソン比

0.2

0.2

A

粘着力 c

4.28×106 N/m2

6.63×106 N/m2

A

内部摩擦角

30 °

30 °

B

圧縮強度 Fc

8.56 N/mm2

23.2 N/mm2

A

引張強度 Ft

1.37 N/mm2

3.32 N/mm2

A

1.13×10-3

5.54×10-3

A

5.04×108 N/m2

2.08×109 N/m2

A

-1.81×10-4-6.64×10-4

A

0.004

0.004

C

-3.62×10-4-1.328×10-3

B

圧縮限界ひずみcu

せん断弾性係数 G
ひび割れひずみ

cr

せん断ひずみ u
引張限界ひずみ

cr0

注：A；室内試験およびその値から求めるパラメータ
B；本論文と同じ実験結果を用い，解析的あるいは文献の条件により求めたパラメータ
C；その他の文献から定めたパラメータ39

（2） RBSM における破壊基準
本論文における降伏条件は，Fig.3.6 (a)および(b)のとおり，モール・クーロンの条件による破壊基準により設定した．Fig.3.6 (a)および(b)では，粘着力は同実験結果を用いている柴田ら（2020）の設定より 4.28×106 N/m2（地山）および 6.63×106 N/m2（覆工）とし，内部摩擦角は柴田ら（2020），吉田ら（2004）の報告より，地山・覆工ともに 30°とした．また，竹内ら（2005）の設定を参考に，0.5Fc を圧縮 1 次降伏とし，Fc までの降伏条件ではτが一定の経路とたどると仮定した（山本ら（2010，2015））．ここで，降伏後は塑性流れ則に従って降伏曲面上を移動し，せん断ひずみが γu=0.004 に達したときにせん断破壊すると仮定した．
ここで，せん断ひずみ u は，実験結果を用いず，文献を参考に決定したパラメータである．
また，垂直バネの圧縮領域の応力解放までのひずみ（2 cu）とその後の応力（0.2Fc），降伏条
件の圧縮 1 次降伏の値（0.5 Fc）は実験結果を用いてはいるものの，乗じている数値は文献を参考に決めている．まず，せん断ひずみは本論文の解析モデルに対する感度が低かったた
めに u=0.004 として特に問題はなかったが，せん断変形が卓越するような条件の場合，せん断強度試験等から得たせん断ひずみを用いる必要がある．また，圧縮挙動に関する 2 cu，0.2Fc，0.5Fc の設定についても，本論文の解析においては感度が低く，結果には特に影響を及ぼさない．しかし，解析対象のひずみが圧縮限界ひずみを超え，異なるポストピーク挙動を示す場合は，一軸圧縮試験で得た結果から改めて設定する必要がある．

（a）地山

（b）覆工
Fig.3.6 降伏条件

40

（3）本研究における要素分割の考え方
剛体ばねモデルは，要素分割（要素の大きさ）の影響を大きく受ける解析手法である．そこで，本研究においては，Fig.3.7 に示す 3 ケースの要素分割について解析を実施し，要素分割の影響について考察する．なお，Case#3 は，最も粗く要素を分割したものであり，Case#2
は，Case#3 よりも要素を細かく分割したもの，Case#1 は，Case#2 を基に，覆工の要素を細かく分割したものである．

（a）要素分割 1（Case#1）（b）要素分割 2（Case#2）（c）要素分割 3（Case#3）

Fig.3.7 要素分割

41

3.1.3 剛体ばねモデルによる解析結果
剛体ばねモデルによる解析は，空洞なし，空洞あり（空洞範囲 45°）
，空洞あり（空洞範囲 90°）の 3 ケースについて実施した．各ケースの解析結果について，以下に整理する．

（1）要素分割の影響
1）検討条件
剛体ばねモデルによる解析は，要素分割（要素寸法）の影響が大きい解析手法であるため，
まず要素分割の影響を確認した．要素分割の影響確認には，Fig. 3.8 に示す共通の荷重条件（模型実験における a-1）を用いて検討した．要素分割は，前述の Fig. 3.7（a）～（c）の 3
ケースについて検討した．

荷重 [N][×105]

10

5

0
0

200

400
時間 [s]

600

Fig. 3.8 模型実験における時間と荷重の関係

42

800

2）要素分割の影響に関する検討結果
Fig. 3.9 （a）および（b）に，要素分割数を変化させた際の載荷板変位と内空変位量の解

析結果を示す．ここで，内空変位量は，上下方向を鉛直方向，左右方向を水平方向とし，それぞれ正と負の値で示している．ここで正の値は，内空変位が広がること，負の値は内空変位が狭くなることを表しており，トンネル覆工が縦長の形状に変形していくことを示す．

（a）載荷板変位の比較（b）

内空変位量の比較

Fig. 3.9 要素分割の影響

43

Fig. 3.9 より，載荷板変位及び内空変位量は要素分割に関わらず，ほぼ同じ過程をたどる

ことがわかる．ただし，本解析では降伏曲面の計算の際，塑性変形が大きくなると，収束計
算が発散して，解析が途中で終了する場合がある．Case #2 および Case #3 では約 311s と351s で解析が中断していることから，Case #1 の要素分割を用いることとした．

（2）剛体ばねモデルによる解析結果（空洞なし）
1）荷重条件
荷重条件は，模型実験により得られた Fig. 3.10 の 3 ケースにより検討した．

荷重 [N][×105]

10

5
a-1

a-2

a-3
0
0

200

400
時間 [s]

600

Fig. 3.10 模型実験における時間と荷重の関係

44

800

2）荷重と載荷板変位の関係
Fig. 3.11（a）～（c）に空洞なしの条件における荷重と載荷板変位の関係を示す．載荷板

変位は，荷重の増加に伴って，比例的に増加し，最大の載荷板変位は，2mm 程度であった．
なお，図に示されるように，解析では，ひび割れが検出されても，微少なうちは，載荷板変位−荷重関係が線形弾性体のように振る舞い，多くの場合，塑性ひずみが大きくなった時点

載荷板変位 [mm]で収束計算が発散して解析が終了する事態が見られた．

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0

5
荷重 [N] [×105]
Fig. 3.11（a）

載荷板変位 [mm]

←ひび割れ発生

a-1

5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0

10

荷重と載荷板変位の関係（a-1）

a-2←ひび割れ発生05
荷重 [N] [×105]
Fig. 3.11（b）

荷重と載荷板変位の関係（a-2）

4510載荷板変位 [mm]

5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0

a-3

←ひび割れ発生05
荷重 [N] [×105]
Fig. 3.11（c）

荷重と載荷板変位の関係（a-3）

46103）荷重と内空変位の関係
Fig. 3.12（a）～（c）に載荷板変位と内空変位の関係を示す．縦軸が縦方向と横方向の内空変位量で，それぞれ正と負の値を示しており，横軸が載荷板変位を表している．側面から荷重を作用させることで，覆工が縦に伸張するように変形し，それにより縦方向の内空変位量は増加している．また，横方向には押しつぶされるように変形するため，内空変位量が減
少していることになる．最大の内空変位は鉛直方向 0.2mm，水平方向程度 1.2mm であり，鉛直方向に比べて，水平方向の変位量が大きくなる傾向が得られた．

内空変位量 [mm]

1.5
鉛直方向

a-1

1
0.5
0
-0.5
-1

ひび割れ発生→

-1.5
0

水平方向

5
荷重 [N] [×105]
Fig. 3.12（a）

10

荷重と内空変位量の関係（a-1）

1.5

内空変位量 [mm]

鉛直方向

a-2

1
0.5
0
-0.5
-1

←ひび割れ発生-1.5
05
荷重 [N] [×105]
Fig. 3.12（b）

荷重と内空変位量の関係（a-2）

47

水平方向

10

1.5

内空変位量 [mm]

鉛直方向

a-3

1
0.5
0
-0.5
-1

ひび割れ発生→

水平方向

-1.5
0

5
荷重 [N] [×105]
Fig. 3.12（c）

10

荷重と内空変位量の関係（a-3）（3）剛体ばねモデルによる解析結果（空洞あり：45°）1）荷重条件
荷重条件は，模型実験により得られた Fig. 3.13 の 2 ケースにより検討した．

荷重 [N][×105]

10

5

b-1

b-2

0
0

200

400
時間 [s]

Fig. 3.13 模型実験における時間と荷重の関係

48

600

2）荷重と載荷板変位の関係
Fig. 3.14（a）～（b）に空洞あり（45°）の条件における荷重と載荷板変位の関係を示す．

載荷板変位は，荷重の増加に伴って，比例的に増加し，荷重が 6.0×105 程度に達した段階
で，収束計算が発散して解析が終了した．すなわち破壊状態が検出された．解析が終了した

載荷板変位 [mm]段階における最大の載荷板変位は，1.3～1.4mm 程度であった．

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

b-1

0

5
荷重 [N] [×105]

載荷板変位 [mm]Fig. 3.14（a）

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

←ひび割れ発生

荷重と載荷板変位の関係（b-1）

ｂ-20←ひび割れ発生5
荷重 [N] [×105]

Fig. 3.14（b）

10

荷重と載荷板変位の関係（b-2）

49103）荷重と内空変位の関係
Fig. 3.15（a）～（b）に荷重と内空変位の関係を示す．変形の傾向は，空洞なしのケースと同様であるが，荷重が 6.0×105 程度に達した段階で，収束計算が発散して解析が終了した．すなわち破壊状態が検出された．解析が終了した段階における最大の内空変位は鉛直方
向 0.3mm，水平方向 0.8mm 程度であり，空洞なしに比べて鉛直方向の変位は大きくなった．

1.5
1
内空変位量 [mm]

鉛直方向

b-1

0.5
0
-0.5
-1
ひび割れ発生→

-1.5
0

水平方向

5
荷重 [N] [×105]
Fig. 3.15（a）

10

荷重と内空変位量の関係（b-1）

1.5
b-2

内空変位量 [mm]

1

鉛直方向

0.5
0
-0.5
-1
ひび割れ発生→

-1.5
0

水平方向

5
荷重 [N] [×105]
Fig. 3.15（b）

荷重と内空変位量の関係（b-2）

50

10

（4）剛体ばねモデルによる解析結果（空洞あり：90°）1）荷重条件
荷重条件は，模型実験により得られた Fig. 3.16 の 2 ケースにより検討した．

荷重 [N][×105]

10

5

c-1

c-2

0
0

200

400

600

時間 [s]
Fig. 3.16 模型実験における時間と荷重の関係

2）荷重と載荷板変位の関係
Fig. 3.17（a）～（b）に空洞あり（90°）の条件における荷重と載荷板変位の関係を示す．

載荷板変位は，荷重の増加に伴って，比例的に増加し，荷重が 5.0～6.0×105 程度に達した段階で，収束計算が発散して解析が終了した．すなわち破壊状態が検出された．解析が終了した段階における最大の載荷板変位は，1.2～1.7mm 程度であった．

51

載荷板変位 [mm]

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

c-1

0

5
荷重 [N] [×105]

載荷板変位 [mm]Fig. 3.17（a）

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

←ひび割れ発生c-2

0

荷重と載荷板変位の関係（c-1）←ひび割れ発生5
荷重 [N] [×105]

Fig. 3.17（b）

10

荷重と載荷板変位の関係（c-2）

52103）荷重と内空変位の関係
Fig. 3.18（a）～（b）に荷重と内空変位の関係を示す．変形の傾向は，空洞なしのケースと同様であるが，荷重が 5.0～6.0×105 程度に達した段階で，収束計算が発散して解析が終了した．すなわち破壊状態が検出された．解析が終了した段階における最大の内空変位は鉛
直方向 1.4mm，水平方向 1.1mm 程度であり，水平方向の変位に比べて鉛直方向の変位の方が大きくなった．

1.5
c-1

内空変位量 [mm]

1

鉛直方向

0.5
0
-0.5
-1
ひび割れ発生→

-1.5
0

水平方向

5
荷重 [N] [×105]
Fig. 3.18（a）

10

荷重と内空変位量の関係（c-1）

1.5内空変位量 [mm]

1

c-2

鉛直方向

ひび割れ発生→

水平方向

0.5
0
-0.5
-1
-1.5
0

5
荷重 [N] [×105]
Fig. 3.18（b）

荷重と内空変位量の関係（c-2）

53

10

3.1.4 剛体ばねモデルによる解析結果の比較
剛体ばねモデルによる解析結果について，時間と荷重，荷重と載荷板変位，載荷板変位と内空変位の関係に着目して整理した．

（1）時間と荷重の関係
Fig. 3.19 に入力条件とした模型実験における時間と荷重の関係を示す．

10

a-2

荷重 [N][×105]

a-3

b-1
c-1

b-2

5

c-2
a-1

0
0

200

400
時間 [s]

600

Fig. 3.19 模型実験における時間と荷重の関係

54

800

（2）載荷板変位と内空変位の関係
Fig. 3.20 に載荷板変位と内空変位の関係を示す．比較の結果，横方向変位は，空洞の有無や範囲に関わらず，載荷板変位の変化に対して，概ね一定の傾きで変化すること明らかとなった．また，縦方向変位についても，載荷板変位に対して概ね一定の傾きで変化するが，空洞なし，空洞あり（45°），空洞あり（90°）の順に内空変位が増大する傾向が得られた．これらの傾向は，本解析においては，すべての解析において同一の物性値及びばねモデル，破壊条件による検討を行っており，ばねモデルについては線形で応力と歪みの関係を仮定しているためであると考えられる．
2

鉛直方向

1.5

内空変位 [mm]

1
0.5
0
-0.5
-1
-1.5

水平方向

-2
0

1

2

3

4

載荷板変位 [mm]
空洞なし

空洞あり（45°）

空洞あり（90°）

Fig. 3.20 模型実験における載荷板変位と内空変位の関係

55

3.2 個別要素法
個別要素法（DEM）は，解析対象を粒子の集まりとして表現することで，連続体では表現しきれない亀裂やせん断といった破壊後の挙動をシミュレーションできる点が特徴である．
以下に，DEM の基礎理論について，概要を整理する．詳細は参考文献「Yade」 <https://yadedem.org/doc/>を参照されたい．

3.2.1 個別要素法の基礎理論
個別要素法は，解析対象を粒子として表現することで，不連続体をモデル化してその挙動を解析することが可能であり，粒子間に結合力を適用することでコンクリートのような連続体についてもモデル化しその挙動を解析できる．
個別要素法の解析で使用する粒子には，①剛体であり変形しない，②剛体であるが粒子の重複を許す，③粒子の重複を許すが重複量は微小であることの 3 つの前提条件が存在す
る．また，粒子の力学的な挙動については，フックの法則（粒子間接触力），Newton の運動の第 2 法則（発生した力に応じた粒子の加速度），リープフロッグ法や Verlet 法などの時間積分手法（加速度に基づく粒子の位置情報や方向の情報の更新）によって表現する．本研究では DEM 解析ソフトとして，オープンソースである Yade を用い，解析は以下の手順で行う．以下の（1）～（4）の操作を繰り返し行うことで解析粒子の逐次的な挙動を得られる．

（1）衝突検出
衝突検出は粒子に境界を設定し，その位置座標をリスト化して，衝突の可能性のある粒子を判定する操作である．衝突の可能性のある粒子間においては，更に接触検出を行うことで計算コストの増大を防ぐことができる．衝突検出は解析のステップごとに行う必要はな
いが．定期的に行う必要があり，大まかな衝突検出と詳細な衝突検出（接触検出）の 2 つがある．衝突検出の概要については 3.2.2 衝突検出に示す．

（2）ひずみ
衝突検出において接触判定がなされた粒子間において，重複領域（貫入量）や，粒子の位置情報から，接触粒子の相対的な変位量やひずみを評価する．ただし，粒子間のひずみや変位は法線方向と接線方向とに区分して考える必要がある．ひずみについての概要は 3.2.3
ひずみに示す。

56

（3）応力各粒子に作用する力を求めるために、フックの法則に基き，接触粒子間において発生する接触力を算出する．ただし，粒子間の力は法線方向と接線方向の力とに区分して考える必要がある．応力の概要は 3.2.4 応力に記す．

（4）運動統合
Newton の運動の第 2 法則に基づき，前述した各粒子に作用する力を，各粒子の質量で割り，各粒子の加速度を算出する．ここで求めた加速度を，任意の時間ステップで時間積分し，粒子の変位量を算出して，粒子の位置情報を更新する．Yade においては，時間積分手法としてリープフロッグ法が使用されている．運動統合の概要は 3.2.5 運動統合
に記す。

57

3.2.2 衝突検出
本節では，衝突検出に関して概要を説明する．ここでは，大まかな衝突検出，衝突の再検
出と粒子境界の拡大，詳細な衝突検出（接触検出）の 3 つの項目について説明する．

（1）大まかな衝突検出
個別要素法では対象物を粒子の集まりとして表現するが，粒子 1 つ 1 つと，その他すべての粒子との接触計算を行うには膨大な計算コストと時間がかかる．そこで，スイープとプ
ルーンアルゴリズムを用いる．具体的には，作成した粒子に対して境界ボックスを設定し，それらを座標系に投影し、粒子境界の最小座標と最大座標を各軸に沿ってリストとして保存する。このリストの中身を昇順に並べ替えることで，衝突の可能性のある粒子と無い粒子とに区分することが可能である． 衝突検出の概念を 2 次元で表現するとのとおりとなる．

Fig.3.21 衝突検出の概念

𝑃 ：粒子番号 i の粒子領域

𝑃 ：粒子番号 i の粒子の境界領域

𝑃 , 𝑃 ：粒子番号 i の粒子の x・y 座標の最小値
𝑃 , 𝑃 ：粒子番号 i の粒子の x・y 座標の最大値

58大まかな衝突検出における衝突の判定基準は、以下のとおりである．
①同一粒子境界の始点と終点が隣り合っていたら、その粒子は衝突していない②ある粒子境界の始点の次に別の粒子境界の始点または終点が存在する場合、その粒
子ペアは衝突している可能性がある。
ただし，これらの判定は軸毎に行われており，3 次元の場合は x 軸，y 軸，z 軸の 3 軸について重複があればその粒子ペアは衝突と判定できる点に留意が必要である．

（2）衝突の再検出と粒子境界の拡大「大まかな接触検出」で説明したように個別要素法では，まず，作成したすべての粒子に対して衝突検出が行われる．しかし、時間ステップが経過するごとに粒子が移動することで，作成した粒子境界の外に粒子が出ていくため，境界条件を満たさなくなる．
そのため、粒子が作成した粒子境界の外に出る前に定期的に衝突の再検出を行いリストの更新を行う必要がある．リストの更新には，2 つの方法がある．それは，以下に示すとおり，すべてのステップにおいて衝突検出を行う方法と固定相対長さを基準として定期的に衝突検出を行う方法の 2 つであるが、いずれも計算回数や計算範囲の拡大によって，計算コストの増大の要素があるため，バランスを取りながら設定する必要がある．
① 全てのステップにおいて衝突検出を実行する方法② 固定相対長さ∆Lを基準とした衝突検出の方法

59

（3）詳細な衝突検出（接触検出）「大まかな接触検出」で示した衝突検出は概算であり，大まかな接触検出で衝突判定された粒子の間での詳細な接触判定が行われる．

Fig.3.22 接触検出の概念

𝑑 = |𝐶̅ − 𝐶̅ |𝑑 =𝑟 +

𝑑 −𝑟 −𝑟2𝐶̅ ：粒子 i の中心座標、𝐶̅ ：接触粒子の接触点,𝑑 = 𝑑 − 𝑑𝑑 ：基準距離（接触判定時(粒子作成時)の接触粒子間距離）𝑑 ：𝑑 測定時の接触粒子 i の中心から接触点𝐶̅ までの距離𝑟 ：粒子 i の粒子半径

衝突検出を行った粒子については，あくまで粒子境界の重複の有無を確認しただけであるので，さらに粒子の間の相互作用の有無を確認する必要がある．
一般的に，DEM においては，ある一つの粒子は直接接触している粒子からのみ影響を受け，接触していないより遠方の粒子から影響は受けない．しかし，相互作用半径を拡大する
ことで，直接接触している粒子以外の粒子との間においても相互作用を考慮することが出来るため，コンクリートのような粒子間において結合力が働いているモデルにおいては，相互作用半径を拡大することで直接接していない粒子間においても結合力を考慮し，粒子の接触数を増大できる．これは，圧縮時の反応に影響を与えたりするために利用される．また，接触の非局所性に影響する．相互作用半径の拡大をした場合は粒子境界の拡大もあわせて行う必要があることに留意が必要である．

60

3.2.3 ひずみ
本節では，ひずみ評価に関して概要を説明する．粒子間に接触が発生した場合，粒子の重複，つまり貫入が生じる。この貫入量は、基準距離𝑑 を基準とした変位量から算出される．
この基準長さが基準となり，接触粒子間に引張力や相互作用の力が発生するため，基準長さを基準とした接触粒子間の相対的な変位量やひずみを算出する必要がある．ただし，求める
ひずみや変位量は，接触粒子間の接触面に垂直法線方向と，接触面内での粒子の接線方向とに分けて算出する必要がある．法線方向ひずみ・接線方向ひずみについて，以下に示す．

（1）法線方向ひずみ
粒子間の接触に伴い発生する力のうち反発力により生じる法線方向のひずみは次のように評価される。

𝑛∘ =𝐶∘ − 𝐶∘≡ 𝐶∘ − 𝐶∘|𝐶 ∘ − 𝐶 ∘ |𝐶∘ = 𝐶∘ + 𝑑 −𝑑 − |𝐶 ∘ − 𝐶 ∘ |𝑛2

𝑢 = |𝐶 ∘ − 𝐶 ∘ | − 𝑑|𝐶 ∘ − 𝐶 ∘ |
𝑢
=
−1
𝑑
𝑑

∘
∘
⎧log |𝐶 − 𝐶 |
⎪
𝑑
𝜀 =
∘
∘
⎨ |𝐶 − 𝐶 |
−1
⎪
𝑑
⎩

𝑖𝑓 |𝐶 ∘ − 𝐶 ∘ | < 𝑑𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒𝑛∘ ：現在のステップの法線方向単位ベクトル

真ひずみ
工学ひずみ

𝐶 ∘：現在のステップの接触点位置

𝐶 °：粒子番号 i の粒子の現在のステップにおける中心座標（位置情報）
𝑢 ：法線方向変位 𝜀 ：法線方向ひずみ

𝑑 ：接触検出時の粒子間距離

接触点𝐶 ∘ は，常に粒子の重複領域の中間に位置し，たとえ重複領域が無い場合でも，粒子間の空白空間の中間に存在する．𝑑 は、基準距離で，粒子間の相互作用が発生し始める初期接触時の粒子間距離である．圧縮時のような大変形においては，ひずみがマイナス∞に向かう傾向があるため，変形前を基準としない真ひずみが使用される．

61

（2）接線方向ひずみ
粒子間の接触に伴い発生する摩擦力（せん断力）により生じる接線方向変位の概念図
を Fig.3.9 に示す．せん断力は粒子に回転を生じさせる働きをする．この回転により初期
接触時において接触点であった点にずれが生じてせん断方向の変位となる．つまり，せん断は，粒子間の法線ベクトルに垂直な並進速度の差と，法線ベクトルに垂直な角速度の和により生じるものである．その他にも粒子間の法線ベクトルに平行な角速度の差によりねじれが，粒子間の法線ベクトルに垂直な角速度の差により曲げが生じる．

Fig.3.23 接線方向変位の概念図

62接線方向のひずみは法線方向のひずみに対して垂直になる，これは，せん断変位𝑢 は、常に，𝑢 + (∆𝑢 ) = 𝑢° ⊥ 𝑛°に従って更新されなければならないことを意味している．

つまり，時間ステップの更新に伴い単位法線方向ベクトルが微小変化（𝑛 → 𝑛° ）し，それに伴い接触面も変化する点を考慮して現在のステップにおけるせん断変位を算出しなければならない．また，変位成分が弾性域の変位成分のみになるように，塑性変位成
分を除いてやる必要がある．これらのことを踏まえて，最終的な現在のステップにおけるせん断変位は次のよう計算される．

𝑢° = 𝑢 + (∆𝑢 ) + (∆𝑢 ) + (∆𝑢 )(∆𝑢 )：せん断変位増分

𝑢 ：前のステップでのせん断変位

𝑛 ：前のステップの単位法線方向ベクトル 𝑛°：現在のステップの単位法線方向ベク

トル度∆𝑡：解析における 1 ステップの時間

𝜔⊖：前のステップにおける粒子 i の角速𝑣 ⊖：前のステップにおける粒子 i の速度𝑣 ：粒子 1 と粒子 2 の相対速度𝑑 ：粒子 i の中心から接触点までの距離
𝑢° ：現在のステップでのせん断変位

63𝑣 ：𝑣 の垂直成分（せん断方向成分）

3.2.4 応力
個別要素法では粒子の接触に伴い発生する相互作用を図のように，接触面と垂直な方向（法線方向）と接触面内に沿った方向（接線方向）とに分けて考える．相互作用モデルにはVoigit モデルが使用されており，粒子間相互作用はフックの法則に基づいて形成される．法線方向と接線方向それぞれの相互作用モデルは Fig.3.24～Fig.3.26 のとおりである。バネが粒子間の反発力や摩擦力を表現し，ダッシュポットがエネルギーの散逸としての減衰を表現し，スライダーは摩擦の調整を表現する。なお，この章では、（1）では，剛性 K の概要と接触粒子間の剛性，（2）では、接触粒子間における剛性の換算式，（3）で，応力計算についての説明を記す。

Fig.3.24 法線方向・接線方向図

Fig.3.25 法線方向バネモデル

64

Fig.3.26 法線方向バネモデル

（1）剛性(バネ剛性)の概要と接触粒子の剛性
先に説明したとおり個別要素法では粒子間の相互作用はフックの法則に基づいて形成される．そのため粒子に作用する力を算出するためには，変位量とバネ剛性の数値が必要である．ただしバネ剛性や変位量は，法線方向と接線方向とに分けて考える必要がある．フックの法則は次のように式変形していくことが可能である．

σ = Eε

F：力 K：剛性（バネ定数）
L：元の長さ

→

F=

𝐸𝐴
δ
𝐿

→

F = K𝛿

𝛿 ：変位量 σ：応力 E：ヤング率 ε：ひずみ

A:面積(接触粒子間においては接触断面π × min (𝑟 , 𝑟 ) である．

ここで，ヤング率 E やポアソン比 ν といった特性値は粒子ごとではなく粒子を構成する部材ごと，つまり巨視的な値である．一方で，法線方向，せん断方向の剛性（バネ剛性）K やK といった特性値は，接触粒子レベルの微視的な値であり，巨視的な数値から計算される．Yade では，次のように，接触粒子の法線方向・せん断方向の剛性を計算し
ている．

𝐾 =

𝐸𝐴
𝐿

𝐾 =𝐾 𝜈
𝐴= π × min (𝑟 , 𝑟 )K ：法線方向剛性（バネ定数） K ：接線方向剛性（バネ定数）
𝐴 ：粒子の接触断面積
𝜈：ポアソン比
𝐿：基準長さ𝑑

𝐸 ：ヤング率

𝑟 ：粒子 i の半径

基準距離𝑑 を設定した際の接触点から接触粒子 i の中心までの距離

法線方向の剛性は，フックの法則に基づいて算出し，接線方向の剛性は法線方向の剛性とポアソン比との乗より算出している．

65

（2）接触粒子間の剛性の換算式
Fig.3.27 に示したような接触粒子の剛性を合成することが可能であり，その計算式は次のようになる．ただし，相互作用半径が 1 で基準長さが衝突粒子の半径の和に等しい
状態を表している．

Fig.3.27 接触粒子間の合成バネモデル∆l = ∆𝑙 + ∆𝑙
𝐾 =𝐸𝑙𝑙 = 2𝑟 𝑜𝑟 𝐴 ⁄𝑙

𝐾 ∆l = F = 𝐹 = 𝐹
𝐾 (∆𝑙 + ∆𝑙 ) = 𝐹
𝐾

𝐹
𝐹
+
𝐾
𝐾

𝐾 =

𝐾𝐾
𝐾 +𝐾

𝐾 =𝐾 𝜈

ここで注意する点は𝑙 の長さである．𝑙 は相互作用作成時，つまり接触検出を行い相互作用が形成された瞬間における接触点から接触粒子までの距離であるので粒子半径と等しくならない場合もある．

66

（3）応力計算
接触粒子間の剛性と変位を算出できたのであれば，フックの法則に従い次のように，粒子に作用する力を求める．ただし、𝑢 > 0の場合(接触粒子が引張方向に移動している場合)，接触は何の力も発生させずに削除される．

𝐹 =𝐾 𝑢 𝑛
𝐹 =𝐾 𝑢

最終的な接線方向力は，上記のトライアルせん断力をもとに次のように表現される．

𝐹 =𝐹|𝐹 | tan 𝜑𝐹𝐹

𝑖𝑓 |𝐹 | > |𝐹 | tan 𝜑𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

𝐹 ：法線方向力 𝐾 ：法線方向合成剛性 𝜇 ：法線方向変位（貫入量）

𝑛：単位法線方向ベクトル 𝐹 ：トライアルせん断力 𝐾 ：接線方向剛性
𝜇 ：接線方向変位量 𝐹 ：接線方向力（せん断力）

せん断力が，接触点で作用する最大静止摩擦力と粘着力（設定する場合）の合計を超えた場合には，滑り（塑性流動）が発生する．ただし，滑りが発生すると接触点でせん
断方向に作用する力はそれ以上増加しない．つまり，接触点で作用する最大静止摩擦力と粘着力（設定する場合）の合計がせん断力の上限値となる．

各粒子は各ステップにおいて，トルクと力を蓄積しており，その際にせん断方向力と
法線方向力の和である概要力（合力）は，両方の粒子に適用される．力は接触粒子が離れた後も接触点 C における作用を計算するので，力により発生するトルクを計算する必要があり，次のように計算される．

𝐹 += 𝐹
𝑇 += 𝑑 (−𝑛) × 𝐹

𝐹 += −𝐹
𝑇 += 𝑑 𝑛 × 𝐹

𝐹 ：概要力 𝑇：トルク

𝑑：接触点から流心中心までの距離(回転軸から見た力の加わる点までの距離)67

3.2.5 運動結合
相互作用計算により粒子に作用する力を算出したら，Newton の運動の第 2 法則に従い，求めた力から粒子の加速度や速度を求め，時間ステップにおける粒子の位置情報の更新を行う．時間積分にはリープフロッグ法が使用されている．リープフロッグ法は加速度・変位・力が時間𝑡 において定義されるのに対して，速度は𝑡⁄2ステップずれた時刻において定義される．その計算過程を以下で示すが，∆𝑡 刻みのパラメータについては上付き文字として+，〇，－が使用され，𝑡 + ∆𝑡, 𝑡, 𝑡 − ∆𝑡の時刻におけるパラメータに対応する．∆t⁄2刻みのパラ
メータについては，それぞれ𝑡 + ∆𝑡⁄2 , 𝑡 − ∆𝑡⁄2の自国におけるパラメータに対応する．運動統合の計算手順は以下のよう進める．

まず，Newton の運動の第 2 法則より次のように加速度を求める．

𝑢̈̈ ° = a = 𝐹 ⁄𝑚

加速度は位置情報(変位量)を２回微分したものに等しく，位置情報の２回微分は次のように示される．

𝑢 ̈≅𝑢 − 2𝑢° + 𝑢∆𝑡

解析においては，過去の情報の蓄積が必要であり，そこから現在の平均速度を求め，それを基に次のステップのための位置情報を計算できる．

𝑢̇ ⊕ = 𝑢̇ ⊖ + 𝑢̈ ° ∆𝑡𝑢 = 𝑢° + 𝑢̇ ⊕ ∆𝑡

このようにして，粒子間相互作用の力をもとにして粒子の位置情報の更新が行われる．

68

3.2.6 本研究における構成モデル
本研究で実施する実験では，地山およびトンネル覆工にモルタルを用いるため，DEM の解析に用いる構成モデルは，コンクリートの解析に有効である CPM（Šmilauer，2010）を採用する．
（1）法線方向の構成式
CPM では，要素に作用する法線方向の応力

N （引張）は，下式で表されるように，計算したひずみの値より求める．

H

1

N

N

kN

N

ここで，H はヘビサイド関数，kN は法線方向の要素間の剛性，

N は法線方向のひずみであり，要素の半径と距離から計算される．また， は損傷変数で，以下の式から得られる．

1

f

0

exp
f

max ~
~
ここで， 0 は限界弾性ひずみ，

2

N

2
1

2

T

f は法線方向の応力－ひずみ関係の軟化勾配に関連するパラメータ， T は接線方向のひずみ（せん断ひずみ）， 1 は損傷に対するせん断ひずみの重み係数であり，本論文では

1 =0 とする．

69

（2）接線方向の構成式
要素に作用する接線方向の応力は下式で計算できる．

kT

T

T

Tp

ここで，kT は接線方向の剛性， Tp は接線方向の塑性ひずみである．降伏関数 f は MohrCoulomb の基準を用い，下式で表される．
f

T

r pl

T

cT

cT 0 1cT

N tan

ここで，c は損傷変数により変化する粘着力，cT0 は初期粘着力， は内部摩擦角である．
なお，塑性流動則は下式で表すことができる．
T

Tp

T

ここで， は塑性乗数である．また，接線方向の応力r pl ˆ T

T

T

T は下式で示される．

ˆT
t

ˆ

T
t
T

ここで， T は正規化された接線方向応力， T は試行応力， T は接線方向ひずみ， T は更新された接線方向ひずみを表す．

70

（3）パラメータの決定と解析条件
割線静弾性係数 E とポアソン比

は，法線方向剛性 kN および接線方向剛性 kT と以下の関係がある．

k N kT
4k N kT

1 kT k N
4 kT k N

Ai Li 2 k N

E
kN

3V

4k N

3kT
kT

ここで，Ai は接触点 i の断面積，Li は長さ，V は要素が配置されている空間の体積である．
次に，圧縮試験の結果と文献（Šmilauer（2010））より，地山の相対延性を
tanφ = 0.577 とする．覆工は同様に，

f

0

f

0

= 1.92，= 3.57，tanφ = 0.577 として解析を行う．引張の軟化部は線形とし，粘着力は c=2.74×106 N/m2（地山）および 6.63×106 N/m2（覆工）とした．
要素半径は 5mm（地山）および 3mm（覆工）とし，約 36 万個一様に配置する．地山の右側面を 2500kN/m2 の荷重制御により載荷し，左側面の水平方向の変位を拘束する．なお内空変位は 20mm×20mm×20mm の領域の変位の平均値を用い，上下と左右の差により算出した．

71

3.2.7 個別要素法による解析結果
個別要素法による解析は，載荷板変位を入力条件として与え，空洞なし，空洞あり（空洞範囲 45°），空洞あり（空洞範囲 90°）の 3 ケースについて実施した．各ケースの解析結果について，以下に整理する．

（1）個別要素法による解析結果（空洞なし）
1）載荷板変位と内空変位の関係
Fig. 3.28 に載荷板変位と内空変位の関係を示す．縦軸が縦方向と横方向の内空変位量で，それぞれ正と負の値を示しており，横軸が載荷板変位を表している．側面から荷重を作用させることで，覆工が縦に伸張するように変形し，それにより縦方向の内空変位量は増加している．空洞なしのケースでは，まず，横方向に押しつぶされるように変形し，その後，縦方向の伸長が発生した．

3

内空変位 [mm]23.E+06

縦方向

背面空洞なし3.E+06

解析結果

1

2.E+060

2.E+06-1

1.E+06-2

5.E+05
横方向

-3
0

1

2
3
載荷板変位 [mm]40.E+00Fig. 3.28 載荷板変位と内空変位量の関係（空洞なし）

72

（2）個別要素法による解析結果（空洞あり：45°）1）載荷板変位と内空変位の関係
Fig. 3.29 に載荷板変位と内空変位の関係を示す．側面から荷重を作用させることで，覆工が縦に伸張するように変形し，それにより縦方向の内空変位量は増加している．空洞ありのケースでも，まず，横方向に押しつぶされるように変形し，その後，縦方向の伸長が発生した．変位量は，空洞なしのケースでは横方向い比べて縦方向の変位がかなり小さな値を示したが，空洞があるケースでは変位の縦横比はほぼ等しくなった。

3

縦方向

内空変位 [mm]23.00E+06

背面空洞45°

2.50E+06

解析結果

1

2.00E+06

0

1.50E+06-1

1.00E+06-2

5.00E+05

横方向

-3
0

0.00E+00

1

2

3

4

載荷板変位 [mm]Fig. 3.29 載荷板変位と内空変位量の関係（空洞あり：45°）73

（3）個別要素法による解析結果（空洞あり：90°）1）載荷板変位と内空変位の関係
Fig. 3.30 に載荷板変位と内空変位の関係を示す．側面から荷重を作用させることで，覆工が縦に伸張するように変形し，それにより縦方向の内空変位量は増加している．空洞ありのケースでも，まず，横方向に押しつぶされるように変形し，その後，縦方向の伸長が発生した．変位量は，空洞なしのケースでは横方向い比べて縦方向の変位がかなり小さな値を示したが，空洞があるケースでは変位の縦横比はほぼ等しくなった。

32

内空変位 [mm]3.00E+06

背面空洞90°

解析結果
2.50E+06

縦方向

12.00E+06

0

1.50E+06-1

1.00E+06-2

5.00E+05

横方向-30.00E+00012
載荷板変位 [mm]34Fig. 3.30 載荷板変位と内空変位量の関係（空洞あり：90°）74

3.2.8 個別要素法による解析結果の比較
剛体ばねモデルによる解析結果について，載荷板変位と内空変位の関係に着目して整理した．
比較の結果，Fig. 3.31 に示すとおり，空洞なし，空洞あり（45°），空洞あり（90°）と空洞範囲が大きくなるほどの順に最大の内空変位が大きくなる．

内空変位 [mm]33.00E+06

縦方向

22.50E+06

1

2.00E+06

0

1.50E+06-1

1.00E+06-2

5.00E+05

横方向-30.00E+00012
載荷板変位 [mm]
空洞あり（45°）
空洞あり（90°）

3

4

空洞なし

Fig. 3.31 模型実験における載荷板変位と内空変位の関係

75

第4章

模型実験に対する数値解析の適用性

本章では，模型実験結果と数値解析結果を比較し，模型実験に対する数値解析の適用性について検討する．

4.1 模型実験に対する剛体ばねモデルの適用性
4.1.1 模型実験結果と解析結果の比較（空洞なし）（1）荷重と載荷板変位の関係
Fig. 4.1（a）～（c）に，荷重と載荷板変位の関係について，模型実験結果と解析結果を比較した．比較の結果， 載荷板変位は最大で 0.4mm 程度の差が生じているものの，実験値と載荷板変位 [mm]解析結果は同様の経過をたどっていることが確認できる．

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

実験（a-1）解析（a-1）
ひび割れ発生（解析）→

ひび割れ発生（実験）→

05
荷重 [N] [×105]

Fig. 4.1（a）

荷重と載荷板変位の関係（a-1）

7610載荷板変位 [mm]54.5

実験（a-2）

4

解析（a-2）

3.5←ひび割れ発生（実験）3

ひび割れ発生（解析）→

2.5
2
1.5
1
0.5
0
0

5
荷重 [N] [×105]

載荷板変位 [mm]Fig.4.1（b）

5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0

10

荷重と載荷板変位の関係（a-2）

実験（a-3）解析（a-3）←ひび割れ発生（解析）

ひび割れ発生（実験）→

05
荷重 [N] [×105]

Fig. 4.1（c）

荷重と載荷板変位の関係（a-3）

77

10

（2）荷重と内空変位量の関係
Fig. 4.2（a）～（c）に，荷重と内空変位量の関係について，模型実験結果と解析結果を比較した．比較の結果，鉛直方向の内空変位量について，実験においてひび割れが発生する段階までは，実験値と解析値がほぼ一致している．一方，ひび割れ発生以降は，両者の違いが大きくなった．この両者の差が大きくなることについて，構成式の観点から考察を加える．
垂直バネの圧縮側では，バイリニアや非線形の構成式を用いることもあるが，本論文では，圧縮強度に至るまで線形で仮定した構成モデルを用いている．また，引張側の応力解放に至る過程も線形で仮定している．内空変位量の評価に用いている覆工の上部と下部は圧縮応力が作用し，覆工のスプリングラインにひび割れが生じており，覆工の構成モデルが鉛直方向の内空変位量の実験値と解析値の差を大きくしているものと考えられる．また，本解析は微小変形理論に基づいているため，ひび割れ発生から応力解放までの解析は可能であるが，その後の大変形解析は検討の範囲外となっている．本論文で検討する対象は，載荷板変位が3mm，内空変位量が 1.5mm までの小さな変位の範囲である．

一方，水平方向の内空変位量は，実験値が載荷直後から微増傾向を示しているのに対し，数値解析では大きく増加しており，実験終了時の値は 0.5mm を超える差がある．水平方向の内空変位量は，覆工の左右両端の 2 点の変位の差から求めたものであり，これら 2 点は載荷板からの距離が異なり変形量も異なることから，実験値と解析結果で差が生じたものと考えられる．先に述べたように，垂直バネの構成モデルは，地山に関しても圧縮強度まで線形で仮定している．また，ひび割れ発生以降も線形としていることから，地山の構成モデ
ルが地山の挙動に影響を及ぼし，その結果，水平方向の内空変位量の実験値と解析値の差にも影響を与えているものと考えられる．

78

1.5

実験（a-1）内空変位量 [mm]

1

鉛直方向

解析（a-1）
ひび割れ発生（解析）→

0.5
0
-0.5
-1

ひび割れ発生（実験）→

-1.5
0

水平方向

5
荷重 [N] [×105]
Fig. 4.2（a）

10

荷重と内空変位量の関係（a-1）

1.5
鉛直方向

実験（a-2）内空変位量 [mm]1

解析（a-2）←ひび割れ発生（実験）

0.5
0
-0.5
-1
ひび割れ発生（解析）→

-1.5
0

5
荷重 [N] [×105]

Fig. 4.2（b）

荷重と内空変位量の関係（a-2）

79

水平方向

10

1.5

内空変位量 [mm]1

実験（a-3）鉛直方向

解析（a-3）←ひび割れ発生（解析）

0.5
0
-0.5
-1
ひび割れ発生（実験）→

水平方向

-1.5
0

2

Fig. 4.2（c）

46
荷重 [N] [×105]8

荷重と内空変位量の関係（a-3）

80

10

4.1.2 模型実験結果と解析結果の比較（空洞あり：45°）（1）荷重と載荷板変位の関係
Fig. 4.3（a）～（b）に，荷重と載荷板変位の関係について，模型実験結果と解析結果を比較した．比較の結果， 解析途中で収束計算が発散し，解析が終了したが，解析が終了した段階における載荷板変位量は最大で 1mm 以上の差が発生した．解析においては，破壊を伴載荷板変位 [mm]う急激な変位が発生した段階で解析が終了したものと考えられる．

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

実験（b-1）

ひび割れ発生（実験）→

05
荷重 [N] [×105]

Fig. 4.3（a）

載荷板変位 [mm]

←ひび割れ発生（解析）

解析（b-1）

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

10

荷重と載荷板変位の関係（b-1）

実験（b-2）解析（b-2）←ひび割れ発生（解析）05荷重 [N] [×105]

Fig.4.3（b）

荷重と載荷板変位の関係（b-2）

81

10

（2）荷重と内空変位量の関係
Fig. 4.4（a）～（b）に，荷重と内空変位量の関係について，模型実験結果と解析結果を比較した．比較の結果，空洞がないケースと比較して鉛直方法の解析結果が実験結果よりもかなり大きな数値を示した．これは，空洞があることによって覆工のアーチ部の要素が空洞部に向かって変形しやすく，実験に比べて大きな変位を示したと考えられる．解析途中で収束
計算が発散し，解析が終了したが，破壊を伴う急激な変位が発生した段階で解析が終了したものと考えられる．

内空変位量 [mm]1.5

実験（b-1）1

鉛直方向解析（b-1）

0.5
0
-0.5
-1

ひび割れ発生（実験）→

-1.5
0

5
荷重 [N] [×105]
Fig. 4.4（a）

1.5

10

荷重と内空変位量の関係（b-1）

実験（b-2）1
内空変位量 [mm]

←ひび割れ発生（解析）
水平方向

鉛直方向

解析（b-2）

0.5
0
-0.5
-1
ひび割れ発生（解析）→

-1.5
0

水平方向

5
荷重 [N] [×105]
Fig. 4.4（b）

荷重と内空変位量の関係（b-2）

82

10

4.1.3 模型実験結果と解析結果の比較（空洞あり：90°）（1）荷重と載荷板変位の関係
Fig. 4.5（a）～（b）に，荷重と載荷板変位の関係について，模型実験結果と解析結果を比較した．比較の結果， 解析途中で収束計算が発散し，解析が終了したが，解析が終了した段階における載荷板変位量は最大で 1mm 以上の差が発生した．解析においては，破壊を伴載荷板変位 [mm]う急激な変位が発生した段階で解析が終了したものと考えられる．

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

実験（c-1）解析（c-1）
ひび割れ発生(解析)→

05
荷重 [N] [×105]

Fig. 4.5（a）

載荷板変位 [mm]

←ひび割れ発生（実験）

5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0

10

荷重と載荷板変位の関係（c-1）

実験（c-2）解析（c-2）
ひび割れ発生(解析)→←ひび割れ発生（実験）05
荷重 [N] [×105]

Fig.4.5（b）

荷重と載荷板変位の関係（c-2）

83

10

（2）荷重と内空変位量の関係
Fig. 4.6（a）～（b）に，荷重と内空変位量の関係について，模型実験結果と解析結果を比較した．比較の結果，空洞がないケースと比較して鉛直方法の解析結果が実験結果よりもかなり大きな数値を示した．これは，空洞があることによって覆工のアーチ部の要素が空洞部に向かって変形しやすく，実験に比べて大きな変位を示したと考えられる．解析途中で収束
計算が発散し，解析が終了したが，破壊を伴う急激な変位が発生した段階で解析が終了したものと考えられる．

2

実験（c-1）内空変位量 [mm]1.5

解析（c-1）

1
0.5
0
-0.5
-1

ひび割れ発生(解析)→-1.5←ひび割れ発生（実験） 水平方向-205
荷重 [N] [×105]
Fig. 4.6（a）

2

鉛直方向

解析（c-2）110荷重と内空変位量の関係（c-2）

実験（c-2）1.5
内空変位量 [mm]鉛直方向

←ひび割れ発生（実験）

ひび割れ発生(解析)→

0.5
0
-0.5
-1
水平方向

-1.5
0

5
荷重 [N] [×105]
Fig. 4.6（b）

荷重と内空変位量の関係（c-2）

84

10

4.1.4 剛体ばねモデルの適用性に関する考察
模型実験に対する剛体ばねモデルの適用性を検討した結果，以下のことが言える．

（1）模型実験結果と剛体ばねモデルによる解析結果の比較により，背面空洞なしのケース
においては，載荷板変位を良好に再現することが示された．一方，水平方向の内空変
位量は過大評価することが確認できた．鉛直方向の内空変位量は，実験における損傷
発生前は良好に一致することが示されたが，損傷発生後は実験値と差が認められた．
要因を構成式の観点から考察すると，本研究では，垂直バネの圧縮側を，圧縮強度に至るまで線形で仮定した構成モデル，引張側の応力解放に至る過程も線形で仮定していることから，ひび割れ発生後に変形モードが変わるものの，傾きがわずかに変化す
るのみでほぼ線形的な変形が発生しているものと考えられる．言い換えると，大きな内空変位を伴う破壊は生じていないこととなる．
（2）背面空洞ありのケースにおいては，載荷板変位，内空変位ともに，実験値よりも大きな傾きで推移し，一定の荷重に達した段階で，収束計算が発散し，解析が終了した．
これは，解析過程において，ひび割れが発生し，応力を開放しながら徐々に変形を続けて，破壊を伴う急激な内空変位が発生した段階で解析が終了したものと考えられる．
本研究で適用した剛体ばねモデルによる解析は，載荷板変位が 3mm，内空変位量が1.5mm 程度までの小さな変位の範囲を対象としており，ひび割れ発生から応力解放までの解析は可能であるが，その後の大変形解析は検討の範囲外となるためである．
（3）剛体ばねモデルによる解析では，背面空洞なしのケースにおいて模型実験と同様に，スプリングライン周辺でひび割れの発生が認められた．ひび割れの数やひび割れの長さでは実験とやや異なる結果が得られたが，ひび割れ発生荷重は良好に一致した．

85

4.2 模型実験に対する個別要素法の適用性
4.2.1 模型実験結果と解析結果の比較
Fig. 4.7（a）～（c）に，個別要素法の解析結果と実験結果（載荷板変位と内空変位量の関係）の比較を示す．図中には実験においてひび割れが生じた載荷板変位の発生位置を示した．
模型実験については，ひび割れ発生段階で，一時的に内空変位が急激に増大し，その後の徐々に内空変位が増大していく傾向が得られた．内空変位増大の幅は，空洞範囲が大きいほど大きくなった．また，数値解析との比較では背面空洞なしと背面空洞 45°のケースにおいて，内空変位量の変化を良好に再現できている．しかしながら，背面空洞 90°のケースでは，
内空変位量の急増する変位量と，増加の傾向は著しく異なっている．
また，実験結果では，Fig.4.7（c）に示すように特に c-1 の結果で内空変位量が大きい値を示しているが，c-2 および背面空洞 45°の実験結果（b-1，b-2）とほぼ同じ内空変位量になる．数値解析では，背面空洞 45°と 90°の縦方向の内空変位量でほぼ同じ値を示しており，
載荷板変位との関係は必ずしも一致しないものの，c-1 のような破壊を伴うような大変形を除いては，ほぼ同じ傾向を示した．このことからも，数値解析では緩やかな内空変位量の変化を再現できる一方，急激な内空変位は再現できていないことが確認できる．
これらのことから，解析では内空変位量の変化を良好に再現できる一方で，スプリングラ
イン近傍のひび割れ発生変位と載荷荷重の関係では解析結果と実験結果との間にやや差があり，ひび割れが発生する際の荷重を定量的に判断することは若干，困難であることがわかる．本研究で用いた構成モデルである CPM ではひずみの値から応力を求め，その値から荷重を計算する，というアルゴリズムを繰り返すが，本研究では右側面から荷重を作用させ，発生した荷重（応力）が徐々に左側に伝達され，ひずみの値も変化していくこととなる．荷重（応力）が伝わる過程においては，実験値と伝達の速度差が生じ，結果としてひずみ・ひび割れの発現にも影響が生じている可能性が考えられる．本研究の条件下での CPM を用いた DEM の解析では，急激に進展するひび割れに関しては必ずしも再現できていないが，他のケースでは解析結果との差が小さいことから，実験結果を良好に模擬していることが確
認できる．

86

3

内空変位 [mm]23.E+06

縦方向

背面空洞なし3.E+06a-1

1

2.E+060

2.E+06
1.E+06-1a-2

5.E+05

-2
a-3

-3
0

1

2
3
載荷板変位 [mm]実験値（a-1）
実験値（a-3）
Fig. 4.7（a）

内空変位 [mm]

3

横方向
4

0.E+00

実験値（a-2）
解析結果

載荷板変位と内空変位の関係（空洞なし）

縦方向

3.00E+06

背面空洞45°2

2.50E+06

1

2.00E+06

0

1.50E+06-1

1.00E+06

b-1-2

5.00E+05

b-2

横方向

-3
0

0.00E+00

1

2

3

4

載荷板変位 [mm]
実験値（b-1）

実験値（b-2）

解析結果
Fig.4.7（b）

載荷板変位と内空変位の関係（空洞あり：45°）87

3

背面空洞90°内空変位 [mm]23.00E+06

c-1

2.50E+06

縦方向

12.00E+06

0

1.50E+06-1

1.00E+06

c-2
5.00E+05-2

横方向
0.00E+00

-3
0

1

2
載荷板変位 [mm]実験値（c-1）
解析結果

Fig.4.7（c）

34

実験値（c-2）

載荷板変位と内空変位の関係（空洞あり：90°）88

4.2.2 個別要素法の適用性に関する考察
模型実験に対する個別要素法の適用性を検討した結果，以下のことが言える．

（1）模型実験結果と個別要素法による解析結果の比較により，空洞なし及び空洞 45°のケースでは内空変位量を良好に再現できたものの，空洞 90°のケースではやや乖離が認められた．一方，ひび割れが発生するときの荷重で比較すると，解析結果と実験結果との間で若干の差があった．このことから，数値解析では内空変位量の変化を再現できるが，スプリングライン近傍のひび割れ発生と載荷荷重の関係ではやや差があり，ひび割れが発生する際の荷重を正確に判断することは困難であることが示唆された．また，地山のスプリングライン近傍から生じるひび割れ性状を明瞭に再現することができた．
（2）個別要素法による農業用水路トンネルの解析では，一部のケースでひび割れの発生する時期に差があるものの，定性的に地山と覆工の変状を良好に評価することが可能であることが示された．覆工と地山の変状の検討を定量的に行う場合，例えば，
最大荷重との相関を明示し，条件を限定して評価することなどが有効であると考えられる．
（3）本研究において，構成モデルとして適用した CPM ではひずみの値から応力を求め，その値から荷重を計算するが，本研究では右側面から荷重を作用させ，発生した荷重（応力）が徐々に左側に伝達され，ひずみの値も変化していくこととなる．荷重（応力）が伝わる過程においては，実験値と伝達の速度差が生じ，結果としてひずみ・ひび割れの発現にも影響が生じている可能性が考えられる．

89

4.3 数値解析の適用性に関する評価
4.3.1 数値解析による模型実験の再現性の特徴
4.1 模型実験に対する剛体ばねモデルの適用性及び 4.2 模型実験に対する個別要素法の適用性より，数値解析による再現性の特徴は以下のとおりである．

（1）剛体ばねモデルの再現性の特徴
1）再現可能な範囲
背面空洞なしのケースにおいては，載荷板変位及び実験における損傷発生前の鉛直方向の内空変位量は良好に再現することができた．
2）再現不可能な範囲
背面空洞なしのケースにおいて，水平方向の内空変位量は過大評価し，実験における損傷発生後の鉛直方向の内空変位量についても実験値と差が認められた．
背面空洞ありのケースにおいては，載荷板変位，内空変位ともに，実験値よりも大きな傾きで推移し，一定の荷重に達した段階で，収束計算が発散し，解析が終了
した．
3）再現性の特徴
背面空洞なしのケースにおいては，載荷板変位及び実験における損傷発生前の鉛直方向の内空変位量は良好に再現することができたことから，変形の小さい範囲は良好に再現することが可能である．
ただし，破壊を伴うような急激な内空変位が発生する段階では収束計算が発散することから，破壊を伴うような大変形に対する適用性は低いといえる．
4）今後の課題
背面空洞なしのケースにおいて，実験における損傷発生後の鉛直方向の内空変位量について，変形モードの変化が再現できなかった．要因として，本研究では解析方法の適用性を検討する目的で，垂直バネの圧縮側を，圧縮強度に至るまで線形で仮定した構成モデル，引張側の応力解放に至る過程も線形で仮定するモデルを適用していることが考えられる．より正確にフィッティングしていくためには，構成モ
デルにバイリニアモデルを適用する等の検討が必要となるが，本件については，今後の検討課題である．

90

（2）個別要素法の再現性の特徴
1）再現可能な範囲
空洞なし及び空洞 45°のケースでは内空変位量を良好に再現できた
2）再現不可能な範囲
空洞 90°のケースでは，実験で発生した急激な内空変位量を再現することができな
かった．
3）再現性の特徴
空洞なし及び空洞 45°のケースでは内空変位量を良好に再現できき，空洞 90°の実験で発生した急激な内空変位は再現できなかったことから，緩やかな内空変位の傾向を良好に再現することは可能であるが，破壊を伴うような急激な内空変位の再現は困難であることがいえる．
4）今後の課題
本研究において，構成モデルとして適用した CPM ではひずみの値から応力を求め，その値から荷重を計算するが，本研究では右側面から荷重を作用させ，発生した荷重（応力）が徐々に左側に伝達され，ひずみの値も変化していくこととなる．荷重（応力）が伝わる過程においては，実験値と伝達の速度差が生じ，結果としてひずみ・ひび割れの発現にも影響が生じている可能性が考えられることから，この点が今後の改善点である．

91

4.3.2 実務における数値解析法の適用性
本研究において，剛体ばねモデル及び個別要素法について模型実験の適用性を評価した
結果，実務における適用性は以下のとおりである．

（1）模型実験結果の再現性
剛体ばねモデルは，背面空洞なしのケースにおいては，載荷板変位及び実験における損傷発生前の鉛直方向の内空変位量は良好に再現することができ，ひび割れ発
生荷重も良好に再現できた．背面空洞ありのケースでは，破壊を伴う急激な内空変
位は再現できなかった．
一方，個別要素法は，背面空洞なし及び空洞範囲 45°のケースにおいて，内空変位を良好に再現できたが，空洞範囲 90°のように破壊を伴う大きな変位が発生する
ようなケースでは再現できなかった．また，緩やかな内空変位は良好に再現できたが，ひび割れ発生荷重にはやや乖離があった．
いずれの解析手法も，急激な内空変位や破壊を伴う大変位は再現できなかった．
ひび割れ発生荷重に着目すると，剛体ばねモデルの方が実験に近い値を示した．

（2）解析の所要時間
剛体ばねモデルは，2 次元モデルであることから，1 ケースの解析は数分から十数分程度である．一方，個別要素法は 3 次元の粒子の移動、接触を扱うため，1 回のケースの解析が数時間から十数時間かかる．

（3）実務における適用性
剛体ばねモデルは，破壊を伴うような大変形が発生する場合は，収束計算が発散し，解析が終了する問題があるが，比較対象の個別要素法と比較して，ひび割れ発生荷重の再現性が高く，解析の簡易性及び解析時間ともに有利であることから，補修・補強対策を検討する実務において，ひび割れ発生荷重の同定に適用できる可能性が示唆された．
一方，個別要素法は内空変位と荷重の関係の再現性が低く，特に破壊を伴うような大変形と荷重の関係の再現性が低いが，荷重との関係を考慮しなければ，内空変

92

位量の再現性は高い．解析の簡易性及び解析時間の問題が解消できれば，内空変位量の詳細解析への適用の可能性がある．

93

第5章

結論

第 1 章では，研究の背景として農業水利施設及び水路トンネルのストック量について整理したうえで，水路トンネルの戦略的な保全管理を推進するためには，適切な時期に補修や補強等の長寿命化対策を実施し，低下した施設性能の回復又は向上を行い施設の長寿命化を図ることが重要であることを述べた．
また，水路トンネルの本来性能である通水性能（水理性能、水利用性能）を維持するためには，構造性能を維持していく必要があり，構造性能の適切な評価のためには，ひび割れの評価が重要であることを述べたうえで．本研究の目的として，以下の 2 点を設定した．
1）矢板工法で施工された農業用水路トンネルを想定し，アーチ頂部の覆工背面に空洞を有する地山と覆工コンクリートを模擬した 1/10 スケールの模型実験によって塑性
圧によるひび割れ発生過程を再現する．
2）模型実験に対する数値解析には，実務レベルをターゲットとして，解析が比較的容易で，ひび割れを良好に再現することが可能と考えられる「剛体ばねモデル（RBSM：
Rigid Body Spring Model）」について，「個別要素法（DEM : Distinct Element Method）」と比較することで適用性を評価する．

第 2 章では，過去の研究等を参考に，塑性圧によるひび割れ発生過程を再現するための模型実験方法を整理し，模型実験の結果について整理した．
実験に用いた模型は， 1/10 スケールの農業用水路トンネルを模擬し，断面寸法は960mm×960mm×400mm の地山を模した供試体の中央に，巻厚 20mm，幅 240mm，高さ 240mmのほろ形のトンネルを配置した．矢板工法で施工したトンネルに発生しやすい背面空洞の影響を検証するために，背面空洞範囲を 0°，45°，90°とする 3 ケースを設定した．
実験の結果，荷重と載荷板変位の関係に着目すると，背面空洞範囲が大きいほど，荷重に対する載荷板変位が大きくなった．これは，側面から荷重を作用させることで，覆工が縦に伸張するように，横方向には押しつぶされるように変形するため，空洞範囲が大きいほど，この変形が顕著になる．この変形に伴い，空洞範囲が大きくなるに伴って，載荷板変位も大きくなることが確認された．
一方，載荷板変位と内空変位の関係に着目すると，空洞なし，空洞あり（45°），空洞あり（90°）の順に内空変位が急激に増大する載荷板変位が大きくなり，空洞なし，空洞あり

94

（45°），空洞あり（90°）の順に最大の内空変位も大きくなることが明らかとなった．
第 3 章では，数値解析法としてひび割れ発生過程の再現性を検討する「剛体ばねモデル（RBSM）」と「個別要素法」のモデルについて解説し，解析結果について整理した．
剛体ばねモデルによる解析の結果，横方向変位は，空洞の有無や範囲に関わらず，載荷板
変位の変化に対して，概ね一定の傾きで変化すること明らかとなった．また，縦方向変位についても，載荷板変位に対して概ね一定の傾きで変化するが，空洞なし，空洞あり（45°），
空洞あり（90°）の順に内空変位が増大する傾向が得られた．
個別要素法による解析の結果を載荷板変位と内空変位の関係に着目して整理した結果，空洞なし，空洞あり（45°），空洞あり（90°）と空洞範囲が大きくなるほど最大の内空変位が大きくなることが明らかとなった．

第 4 章では，模型実験に対する剛体ばねモデルと個別要素法の再現性を整理し，実務レベ
ルでの数値解析の適用性について整理した．
剛体ばねモデルは，個別要素法と比較して，ひび割れ発生荷重の再現性が高く，解析の簡易性及び解析時間ともに有利であることから，補修・補強対策を検討する実務において，ひび割れ発生荷重の同定に適用できる可能性が示唆された．
一方，個別要素法は内空変位と荷重の関係の再現性が低く，特に破壊を伴うような大変形と荷重の関係の再現性が低いが，荷重との関係を考慮しなければ，内空変位量の再現性は高いことから，解析の簡易性及び解析時間の問題が解消できれば，内空変位量の詳細解析への適用の可能性があることを示した．

95

参考文献
朝倉俊弘，小島芳之，安東豊弘，佐藤豊，松浦章夫（1994a）：トンネル覆工の力学挙動に関する基礎的研究，土木学会論文集，493/III-27，79-88.
朝倉俊弘，安東豊弘，小俣富士夫，若菜和之，松浦章夫（1994b）：欠陥を有するトンネル覆工の変形挙動と内面補強工の効果，土木学会論文集，493/III-27，89-98.
Beckmann, B., Schicktanz, K., Reischl, D. and Curbach, M. (2012) : DEM simulation of concrete
fracture and crack evolution, Structural Concrete, 13(4), 213-220.
Boon, C.W., Houlsby, G.T. and Utili, S. (2015) : Designing tunnel support in jointed rock masses via
the DEM, Rock Mechanics and Rock Engineering, 48, 603-632.

土木学会エネルギー委員会新技術・エネルギー小委員会水路設備保全技術の実務者育成に関する調査・研究分科会（2021）：水路トンネル維持管理の手引き，土木学会．
Funatsu, T., Hoshino, T., Sawae, H. and Shimizu, N. (2008) : Numerical analysis to better understandthe mechanism of the effects of ground supports and reinforcements on the stability of tunnels using
the distinct element method, Tunnelling and Underground Space Technology, 23, 561-573.

濱崎英作，竹内則雄，大西有三（2006）：三次元斜面安定問題に対する簡易離散化極限解析法の開発，日本地すべり学会誌，42(5)，389-397．
伊佐治

優，山本佳士，中村光，三浦泰人（2016）：有限回転を考慮した 3 次元 RBSM によるコンクリート材料のポストピーク挙動解析，コンクリート工学年次論文集，38(2)，91-96．
石橋孝治，松本嘉司（1983）：ひずみ軟化特性を考慮した極限解析法によるトンネル掘削問題の解析，土木学会論文報告集，331，103-111．
石井将幸，長束勇，緒方英彦，円城寺将貴，上見謙太（2014）：水路トンネルの安全性評価を目的とした現場載荷試験と数値解析によるその評価，農業農村工学会論文集，292，41-48.
石故有生也（2017），トンネル挙動解析における DEM の適用性とロックボルトの補強効果の検証，平成 28 年度岡山大学修士論文
蒋 宇静，肖 俊，山口耕平，棚橋由彦，江崎哲郎（2001）：大規模地下空洞の掘削に伴う局所破壊と補強設計に関する研究，資源と素材，117，639-644.
Jifei, W., Hongwei, H., Xiongyao, X. and Bobet, A. (2014) : Void-induced liner deformation and stressredistribution, Tunnelling and Underground Space Technology, 40, 263-276.
Khanal, M., Schubert, W. and Tomas, J. (2005) : DEM simulation of diametrical compression test on
particle compounds, Granular Matter, 7(2-3), 83-90.

96

近藤観慈（1996）：間隙水圧を考慮した剛体ばねモデルによる地すべり斜面安定解析と全体安全率，地すべり，33(2)，1-9．
Kozicki, J. and Donzé, F.V. (2008) : A new open-source software developed for numerical simulations
using discrete modeling methods, Computer Methods in Applied Mechanics and Engineering, 197,
4429-4443.

森 充広，浅野 勇，渡嘉敷 勝，西原正彦（2013）：トンネルの構造的安全性に関する研究動
向～変状事例に基づく考察～，農工研技報，214，31-76.
森

充広，浅野

勇，渡嘉敷

勝，川上昭彦，川邉翔平（2015）：覆工背面に空洞を有するトンネル覆工および周辺地盤の変形挙動に関する模型実験，農業農村工学会論文集，300，
I_195-I_205.

森

充広，浅野

勇，川邉翔平（2021）：覆工背面に空洞が存在する馬蹄形水路トンネルの破壊試験，農業農村工学会論文集，312，II_17-II_25.
Nagai, K., Sato, Y. and Ueda, T. (2004) : Mesoscopic Simulation of Failure of Mortar and Concrete by
2D RBSM, Journal of Advanced Concrete Technology, 2(3), 359-374.

長束勇，村尾弘道，松本拓，石井将幸（2013）：水路トンネルのひび割れ発生メカニズムとその補強工法，農業農村工学会論文集，288，77-83.
Nitka, M. and Tejchman, J. (2015) : Modelling of concrete behaviour in uniaxial compression andtension with DEM, Granular Matter, 17(1), 145-164.
Riera, J.D., Miguel, L.F.F. and Iturrioz, I. (2016) : Evaluation of the discrete element method (DEM)and of the experimental evidence on concrete behaviour under static 3D compression, Fatigue &Fracture of Engineering Materials & Structures, 39(11), 1366-1378.
柴田俊文，田本敏之，西村伸一，珠玖隆行，福元豊（2020）: 農業用水路トンネルの模型実験に対する個別要素法の適用性，農業農村工学会論文集，310，I_105-I_116.
嶋本敬介，野城一栄，小島芳之，中西祐介，塚田和彦，朝倉俊弘（2013）：地山の塑性化に伴う山岳トンネルの路盤隆起現象とその対策工に関する研究，土木学会論文集 F1（トンネル工学），69 （1），54-72.
城野忠幸，佐藤昌志，佐藤京（1994）：トンネルにおける離散化極限解析に関する一考察，

開発土木研究所月報，496，10-21.
Si, X., Huang, L., Li, X., Ma, C. and Gong, F. (2021) : Experimental Investigation of Spalling Failureof D‑Shaped Tunnel Under Three‑Dimensional High‑Stress Conditions in Hard Rock, Rock97

Mechanics and Rock Engineering, 54, 3017-3038.
Šmilauer, V. (2010) : Cohesive Particle Model using the Discrete Element Method on the Yade
Platform, Ph. D. thesis, Czech Technical University in Prague and Université Joseph Fourier.

竹内則雄，粂田俊男，川井忠彦（1981）：新離散化モデルによる地盤基礎の極限解析（その5）－トンネル問題－，生産研究，33（5），182-185．

竹内則雄，上田眞稔，上林厚志，鬼頭宏明，斉藤成彦，冨田充宏，樋口晴紀（2005）：鉄筋
コンクリート構造の離散化極限解析法，丸善．
Toi, Y. and Yoshida, S. (1991) : Numerical simulation of nonlinear behaviors of two-dimensional block
structures, Computers & Structures, 41(4), 593-603.

「Yade」 <https://yade-dem.org/doc/>
山本佳士，中村 光，黒田一郎，古屋信明（2010）：3 次元 RBSM による横拘束コンクリー
トの 1 軸圧縮破壊解析，土木学会論文集 E，66(4），433-451．
山本佳士，岡﨑宗一郎，別府万寿博，中村光，三浦泰人（2015）：RBSM による高速衝突を受けるコンクリート板の局部破壊挙動評価に関する基礎的研究，構造工学論文集，61A，945-957．
Yao, C., He, C., Jiang, Q., Shao, J. and Zhou, C. (2022) : A modified rigid-body-spring method formodeling damage and failure of brittle rocks subjected to triaxial compression, Computers andGeotechnics, 152, 105046.
野城一栄，平田 亮，岡野法之，小島芳之（2015）: 種々の材料からなる山岳トンネル覆工の変形破壊挙動に関する研究，土木学会論文集 F1（トンネル工学），71(2），78-94.
保田尚俊，塚田和彦，朝倉俊弘（2016）：背面空洞がトンネル覆工に及ぼす力学的影響に関する理論的考察，土木学会論文集 F1（トンネル工学），72(3），I_63-I_73.
Yasuda, N., Tsukada, K. and Asakura, T. (2019) : Three-dimensional seismic response of a cylindricaltunnel with voids behind the lining, Tunnelling and Underground Space Technology, 84, 399-412.

横山幸也，志水義彦，石田毅，水田義明（2001）：応力測定によるトンネルライニングの健全度評価に関する実証的研究，土木学会論文集，687/III-56，65-76.
吉田幸夫，水野英二，畑中重光（2004）：異なる形状比を有する円柱コンクリートの三次元
FEM 一軸圧縮解析，コンクリート工学年次論文集，26(2），19-24.

98謝辞まず，本論文をとりまとめるにあたり，丁寧なご指導を賜りました岡山大学大学院環境生
命科学研究科教授

西村伸一

博士，同教授竹下祐二

博士，同教授小松満

博士に感謝の意を表します．
私は，2014 年 10 月に岡山大学大学院に入学しましたが，仕事との両立が難しく，3 年半の在学及び 3 年の休学を経て，2021 年 3 月に単位取得退学をしました．その間に大病を患ったこともあり，一時は学位取得を諦めかけましたが，先生方の温かいご指導のおかげで何とか論文をまとめることができました．もう一度頑張る気持ちにさせてくれた先生方に心より感謝しています．
また，本論文の審査をお引き受け下さり，貴重なご意見を頂きました岡山大学大学院環境
生命科学研究科准教授

柴田俊文

博士，同准教授珠玖隆行

博士に心より感謝いたします．
大学院生活では博士課程同期入学の金重稔博士，模型実験や解析を一緒に行った石故有生也君，村上諒君，事務連絡等の大学院生活のフォローをしていただいた楢崎良子さんをはじめ，多くの方々にお世話になりました．お世話になった皆様に心から厚くお礼申し上げます．
日化エンジニアリング株式会社の二越文彦社長，若林孝技師長（博士）には学位取得に挑戦するチャンスをいただきました．若林孝技師長とは本研究のターゲットである水路トンネル関連の多くの業務を担当し，議論のなかで多くのヒントを得ることができました．また，
同社中国四国支社の皆さんには，研究期間中の業務遂行にご協力をいただきました．ここに謝意を表します．
最後に，「こんなに田舎で育って博士とは，とぼけたことを言うなあ」と言いながらも，体調を気遣いながら温かく見守ってくれた家族に感謝します．
田本敏之