逆解析手法を用いたトンネル掘削面前方の地質予測
樋川 敦∗ , 川原 睦人∗ , 金子 典由∗∗∗Department of Civil Engineering, Chuo University Kasuga 1-13-27, Bunkyo-ku, Tokyo 112–8551, JAPAN
E-mail : [email protected] E-mail : [email protected]
∗∗SATO KOGYO CO., LTD.
12-20, Nihonbashi-Honcho 4-Chome, Chuo-ku, Tokyo 103-8639, Japan
Abstract
This paper presents a numerical technique to predict a soft stratum ahead of an excavating face of a tunnel. The functional method is a parameter identification of the elastic property which is important coefficient at the actual site. The parameter identification is an useful technique in the reverse analysis to determine an unknown parameter using the observed data. The parameter iden-tification is utilized to minimize the difference between the computed and observed displacements to determine the Young’s modulus of the ground. In order to solve the minimization problem, the conjugate gradient method is applied. As a numerical model, the 3-D linear elastic body with the finite element method is used to excavate the natural ground and to identify the Young’s modulus ahead of a tunnel face. Using the 3-D tunnel, the displacement of the axial direction can be con-sidered.
Key words: Parameter identification, 3-D tunnel, Excavation, Finite element method, Conjugate gradient method
1
はじめに
近年、トンネルの掘削技術は日々、進歩し続けている。しかし、トンネルの掘削現場では、危険と隣合わせで もある。その問題の一つに、地山の強度に関する問題がある。地山は多くの断層で構成されており、その断層は、 強度の異なるものが縦横複雑に入り組んで存在している。もし、トンネル切羽が軟弱断層にぶつかると、トンネ ル掘削面は危険な状態となる。これまでにも、多くの現場でその原因による崩落事故が起こっている。その崩落 事故を防ぐ最も効果的な方法は、事前に切羽前方の地質状況を把握しておくことである。そして、そのような効 果的な方法を容易に行えることが期待されている。 従来では、ボーリング孔やＴＳＰ（弾性波探査）と呼ばれる方法により、切羽前方の地質予測が行われてきた。 しかし、これらの方法は、特殊な計測機器を必要とするため、コストがかかり、施工計画を遅らせるという問題 があり、日常的な施工管理の一環として行えない欠点があった。つまり、特殊な計測機器を使わずに予測が行える 方法が期待されていた。 そこで、新しい切羽前方探査の方法として、数値解析を用いた手法が考えられる。もし、地盤の強度を表す弾 性係数が事前に予測、決定できれば、切羽前方の地質状況をデジタル的に評価できる。その効果的な方法として、 逆解析手法の一つである弾性係数の同定解析がある。同定解析は観測データを用いて、未知係数を決定する方法 である。 この研究報告書では、３次元トンネルモデルを用いた掘削現場における地盤弾性係数の同定解析を提案してい る。３次元トンネルモデルを用いることにより、切羽面の挙動が正確に表現できる。今回の解析で、同定される べきモデルとして、切羽前方の地盤はいくつかの断層に分割されている。また、その断層は切羽面に対して垂直 に分割されている。数値解析例では、この手法の妥当性を検証するために、まず、理想モデルが用いられている。 つまり、順解析によって計算された変位を計測値として用いて、そのモデルにおける切羽前方の断層の弾性係数を 予測している。次に、実際の現場で計測されたトンネル内空変位データを用いて、実際に観測された切羽前方の 軟弱断層を予測している。2
基礎方程式
この報告書では、総和規約によって数式を表記している。 まず、応力のつり合い方程式は次式で表現される。 σij,j− ρbi = 0 (1)ここで、σij, bi, ρ はそれぞれ、全応力、物体力、地盤の単位体積重量を表している。 また、変位・ひずみ方程式は次式で表現される。 εij =− 1 2(ui,j+ uj,i) (2) ここで、εijとui はそれぞれ、ひずみ、変位を表している。また、土の構成則として以下の式が導入される。 σij = Deijklεkl (3) ここで、De ijkl は弾性応力ひずみ行列であり、以下の式で表される。 Deijkl= λδijδkl+ µ(δikδjl+ δilδjk) (4) ここで、δijはクロネッカーデルタ記号であり、λ と µ はラメの定数である。 λ = νE (1− 2ν)(1 + ν) (5) µ = E 2(1 + ν) (6) ここで、E は弾性係数、ν はポアソン比である。 境界S は SU とST に分けられ、その境界上では以下の境界条件が満たされている。 ui= ˆui on SU, (7) ti= σijnj = ˆti on ST, (8) ここで、ˆuiと ˆtiは境界上の既知量を表し、niは境界の法線方向ベクトルである。
3
有限要素方程式
有限要素法を適用することで、線形の四面体要素を用いた離散化式が得られる。 Kαiβkuβk= ˆΓαi (9) Kαiβk= V(Nα,jD e ijklNβ,l)dV (10) ˆ Γαi= V (Nαρ ˆbi)dV − ST (Nαtˆi)dS (11) ここで、Nαiは有限要素法の補間関数である。4
掘削解析
図１はこの本研究での掘削手法の図式を表している。 F F p p (0) (1)σ =
d
σ
(0)u
=
du
(1)σ =
(1)d
σ
(1)+
(0)σ
(2)u
=
du
(1)σ =
(2)d
σ
(2)+
(1)σ
+
(1)u
図１．掘削手法 掘削解析は実際のトンネル現場で計測された変位をシミュレーションするために用いられる。そのアルゴリズム が以下に示されている。ただし、n は掘削の回数を表している。 1. 地山周辺等圧を用いて弾性解析を行い、n = 0 とする。 2. 初期変位 u(0)i = du(0)i と初期応力σ(0)ij = dσij(0)を計算する。 3. 初期応力 σ(0)ij から等価節点力p(n+1)αi を計算する。 4. 掘削境界上に等価節点力 p(n+1)αi を与える。 5. 等価節点力 p(n+1)αi を用いて掘削後の弾性解析を行う。 6. 掘削による増分変位 du(n+1)i と増分応力dσ (n+1) ij を計算する。 7. 全変位 u(n+1)i と全応力σ (n+1) ij を計算する。 8. 全応力 σij(n+1)から等価節点力p (n+2) αi を計算する。 9. n = n + 1 とし、4 へ戻る。 掘削解析では、変位と応力の更新、等価節点力の算出が重要である。 また、等価節点力pαiは以下の式で計算される。 p(n+1)αi = V Nα,jσ (n) ij dV, (12) ここで、σij(n)は前回の掘削までの全応力である。
5
トンネル支保工
実際のトンネル掘削表面は、トンネル支保工により補強されている。そして、トンネル支保工は多くの種類があ る。本研究では、アーチリブ支保工が適用されている。なぜなら、支保工が梁に置き換えて計算できるためであ る。また、３次元の座標変換マトリックスT は梁の部材座標系から全体座標系に変換するものである。図２は、そ の座標変換を表している。ここで、３次元座標変換マトリックスT は以下の式で表現される。 X Y Z A B B’ α β 図２．３次元梁の座標変換図 T = cos α cos β 0 cos α sin β 0 sin α 0 0 cos α cos β 0 cos α sin β 0 sin α . また、図３は３次元梁によって表現されたアーチリブ支保工を表している。X
Z
Y
EXCAVATION 図３．梁によるアーチリブ支保工6
パラメーター同定
弾性係数の推定は、同定解析手法を用いて行われる。つまり、評価関数J の最小化問題である。その評価関数 J は弾性係数に関係する変位の計算値と計測値の残差の２乗和で表現される。 J = V 1 2(u − u ∗)TQ(u − u∗)dV (13) ここで、u と u∗はそれぞれ計算変位と計測変位を示している。つまり、ここでの問題は評価関数を最小にするよ うな弾性係数を求めることである。 最小化問題を解くために、共役勾配法が適用されている。この方法では、ポアソン比が常に一定値であり、弾性係 数が常に更新されて求められる。 {P }T ={E}T ={E 1, E2, E3, · · ·Em}, (14) ここで、m は地盤の断層数、E は各層の弾性係数である。 以下に、共役勾配法のアルゴリズムを示す。 1. 初期パラメーター P(0)を仮定し, 評価関数の許容誤差εjとi = 0 を設定する。 2. 初期状態量 u(P )(0)を計算し、初期評価関数J(0)を計算する。 3. 初期感度行列 ∂u(P ) ∂P (0) を計算する。 4. 評価関数の初期勾配{d}(0)=−∂P∂J (0)=− ∂u(P ) ∂P T (u(P ) − u(P )∗) を計算する。 5. パラメーターを更新するために、α = − {d} T ∂J ∂P {d}T ∂u(P ) ∂P T ∂u(P ) ∂P {d} を計算する。6. P(i+1)= P(i)+ αd(i)によりパラメーターを更新する。 7. 状態量 u(P )(i+1)を計算し、評価関数J(i+1)を計算する。
8. 感度行列 ∂u(P ) ∂P (i+1) を計算する。
9. 評価関数の勾配を更新するために、β = ∂J ∂P (i+1) , ∂J ∂P (i+1) ∂J ∂P (i) , ∂J ∂P (i) を計算する。 10. 評価関数の勾配{d}(i+1)=−∂J ∂P (i+1) + β{d}(i)を更新する。 11. J(i+1)< εj で計算終了。 12. それ以外は i = i + 1 として [5] へ戻る。 ただし、簡便にするために、総和規約は無視している。
7
地質予測システム
ここでは、地山で計測された変位を用いた切羽前方の断層の地質を予測する方法を提案している。 ステップ１ : 順解析 ( 地山の掘削解析 ) 1. 地山全体から計算領域を決定する。 2. その領域の弾性係数を仮定する。 3. 掘削解析によって掘削外力 F を求め、変位と応力を計算する。 ステップ２ : 逆解析 ( 弾性係数の同定解析 ) 1. 各断層の初期弾性係数 Einitを仮定する。 2. トンネル表面上の計測値を入力する。 3. ステップ１で求められた F を用いて同定解析を行う。 4. 各層の弾性係数を決定する。 通常の山岳トンネルでは大きな山で施工されるため、地山と同じスケールのモデルを用いることは困難である。そ のため、トンネルを含む一部分が計算領域に用いられる。図４は、トンネル切羽を含む計算領域を表している。本 研究では、４分の１モデルが有限要素解析に用いられている。地山周辺等圧F = γh を用いた掘削解析がステップ １で行われる。ここで、γ は岩の単位体積重量、h はトンネルの被りを示している。ステップ２では、切羽前方の 地盤がいくつかの断層に分割されている。そして、弾性係数の初期値Einitが各層で仮定されている。計測値とし て、掘削による増分変位が使われている。図５は、掘削による増分変位の算定方法を示している。まず、掘削前の 計測点ＡとＢにおける初期座標値が計測される。そして、一回掘削された後、再び、同じ計測点ＡとＢにおける 座標値が計測される。その座標値の差が計測点ＡとＢにおける掘削増分変位となる。また、各層の弾性係数は共 役勾配法の逆解析アルゴリズムに従って決定される。 Tunnel face Computational domainExcavation Natural ground
A B ( x1,y1,z1 ) ( x1,y1,z1 ) Origin (a) 初期座標値 Origin B A ( x2,y2,z2 ) ( x2,y2,z2 ) Excavation (b) ２回目の座標値 B A Excavation dUx dUy dUz dUx dUx = x2 - x1 dUy = y2 - y1 dUz = z2 - z1 (c) 増分変位 図５．増分変位の算定方法
8
数値解析例１
ここでは、本研究のトンネル切羽前方の断層の地質予測方法の妥当性を検証している。評価関数の既知量はトン ネル内空の増分変位が用いられている。ただし、この解析例では、順解析で掘削計算された増分変位が計測値に 使われている。つまり、順解析で用いた地山モデルの各層の弾性係数を予測されることが結果として求めるべき 答えである。 ここでは、トンネル切羽前方の軟弱断層の予測を検証するために、２つのケースが解析されている。この解析例で は、健全な硬い層と軟弱な層の弾性係数はそれぞれ、3.0 × 105[kN/m2] と 3.0 × 104[kN/m2] に設定した。この値 は、任意に設定したが、実際の地山物性値とあまり差はなく、適当な値である。また、トンネルの径は 4.0[m] に 設定した。図６は、この解析例で用いられた有限要素分割図を表している。図からもわかるように、ここでは手 法の有効性を確かめるために行われているので、モデルは実際の地山とは異なり、小さめのスケールである。そ のため、軟弱断層の影響を反映しやすいモデルとなっている。(a) 初期地山状態の有限要素モデル (b) 掘削された２次地山状態の有限要素モデル 図６．有限要素分割図
8.1 ケース１
このケースでは、切羽前方の軟弱断層が厚く存在している場合を想定している。図７は、ケース１で予測される べきモデルを表している。層厚は６ｍで、トンネル切羽面から軟弱断層までの距離は４ｍである。表１は、各層 の弾性係数の目的値である。各層の弾性係数の初期値は 1.0 × 104[kN/m2] に設定した。 14m 12m 16m 20m 0m 10m E1 E2 E3 E4 図７．予測されるべき地山モデル 弾性係数 地質状況 E(kN/m2) E1 3.0 × 105 健全 ( 0m - 12m ) E2 3.0 × 105 健全 (12m - 14m ) E3 3.0 × 104 軟弱 (14m - 16m ) E4 3.0 × 104 軟弱 (16m - 20m ) 表１．各層の弾性係数の目的値 図８は、繰り返し計算による各層の弾性係数の履歴を表している。トンネル切羽前方の軟弱断層はＥ３とＥ４に 位置している。すべての層の弾性係数は任意に決めた初期値から目的値に収束している。トンネル掘削表面上の 計測値が、トンネル切羽前方の軟弱断層に影響を与え、十分に反応したといえる。0 50000 100000 150000 200000 250000 300000 350000 400000 0 50 100 150 200 250 300 350 400 450 Youngs modulus [kN/m2] Iteration 0m - 12m ( E1 ) 12m - 14m ( E2 ) 14m - 16m ( E3 ) 16m - 20m ( E4 ) 図８．ケース１の各層の弾性係数の履歴
8.2 ケース２
このケースでは、トンネル切羽前方の軟弱断層は薄い層を想定している。図９は、ケース２で予測されるべき地 山モデルを表している。軟弱断層の幅は２ｍに設定した。表２は、各層の弾性係数の目的値を表している。各層 の弾性係数の初期値は 1.0 × 104[kN/m2] に設定した。 14m 12m 16m 20m 0m 10m E1 E2 E3 E4 図９．予測されるべき地山モデルYoung’s modulus Stratum
E(kN/m2) E1 3.0 × 105 Hard ( 0m - 12m ) E2 3.0 × 105 Hard (12m - 14m ) E3 3.0 × 104 Soft (14m - 16m ) E4 3.0 × 105 Hard (16m - 20m ) 表２．各層の弾性係数の目的値 図１０は、各層の弾性係数の繰り返し計算による履歴を表している。ケース２では、薄い軟弱断層であるため、実 際の現場により近い状況を想定している。１４ｍから１６ｍの間の断層の弾性係数が軟弱断層の目的値に収束し、 薄い軟弱断層であっても、トンネル掘削表面上の変位のみで予測できたといえる。この結果、理想モデルでの切 羽前方の軟弱断層の予測が可能であることが確認でき、この手法の有効性を検証できた。
0 50000 100000 150000 200000 250000 300000 350000 400000 0 50 100 150 200 250 Youngs modulus [kN/m2] Iteration 0m - 12m ( E1 ) 12m - 14m ( E2 ) 14m - 16m ( E3 ) 16m - 20m ( E4 ) 図１０．ケース２の各層の弾性係数の履歴
9
数値解析例２
数値解析例２では、実際のトンネル現場で計測された変位データを用いた切羽前方の地質予測を議論している。図 １１は、予測されるべき地山での計算領域の分割を表している。計算領域は２つの部分に分割した。領域１は硬 い断層のみで構成されており、軟弱断層は存在していない。一方、領域２は切羽前方に軟弱断層が存在していた ことが実際の調査で確認されている。表３は、実際の山岳トンネルのデータを表している。この解析例では、土 被りによる地山圧力F = γh が自重の代わりに初期地山状態に載荷されている。２つの領域で、モデルは切羽が 100m から 105m まで掘削されたとき、同定解析が行われている。この解析例で用いられている有限要素分割図は 図１２に示されている。数値解析例１とは異なり、トンネル軸方向のスケールが大きくなり、軟弱断層の影響を考 慮しにくいモデルになっている。 0m 200m 0m 200m Domain1 Domain2 Sand Stone Soft layer 30m Figure 11. 計算領域の分割図(a) 初期地山状態の有限要素モデル (b) 掘削された２次地山状態の有限要素モデル 図１２．有限要素分割図 トンネル径 土被り 岩の単位体積重量 地山圧力 D(m) h(m) γ(t/m3) F (kN/m2) 10.5 24.75 2.2 539 表３．実際のトンネルデータ
9.1 領域１
領域１では、トンネル切羽前方の軟弱断層は存在していない。図１３は、領域１で予測されるべき地山モデルを 示している。表４は、各層の弾性係数の目的値である。健全な硬い断層の弾性係数は 2.45 × 105[kN/m2] に設定 した。各層の弾性係数の初期値は 1.0 × 105[kN/m2] に設定した。この領域では、切羽前方に軟弱断層は存在して いないので、各層のすべての弾性係数は健全な硬い層の目的値に収束するはずである。 0m 105m 110m 130m 150m 200m E1 E2 E3 E4 Measured point : 100m : Hard ground 図１３．予測されるべき地山モデル 弾性係数 地質状況 E(kN/m2) E1 2.45 × 105 Hard ( 0m - 110m ) E2 2.45 × 105 Hard (110m - 130m ) E3 2.45 × 105 Hard (130m - 150m ) E4 2.45 × 105 Hard (150m - 200m ) 図４．各層の弾性係数の目的値 図１４は、領域１での各層の弾性係数の繰り返し計算による履歴を表している。分割したすべての断層の弾性係 数が健全な硬い断層の目的値に近い値に収束していることがわかる。つまり、この領域では観測観測通り、切羽前 方には軟弱断層が存在していなかったと予測できた。0 50000 100000 150000 200000 250000 300000 350000 400000 0 20 40 60 80 100 120 140 Youngs modulus [kN/m2] Iteration 0m - 110m ( E1 ) 110m - 130m ( E2 ) 130m - 150m ( E3 ) 150m - 200m ( E4 ) 図１４．領域１の各層の弾性係数の履歴
9.2 領域２
領域２では、トンネル切羽前方に軟弱断層が存在している。実際の現場では、切羽前方の軟弱断層の影響で崩落 事故が起こっている。図１５は、領域２で予測されるべき地山のモデルを示している。表５は、各層の弾性係数 の目的値を表している。軟弱断層の弾性係数は健全な断層の弾性係数の１０分の１に設定した。各層の弾性係数 の初期値は 1.0 × 105[kN/m2] に設定した。トンネル切羽前方の弾性係数は、計測データに反応して、軟弱断層の 目的値に収束するはずである。 0m 105m 110m 130m 150m 200m E1 E2 E3 E4Measured point : 100m : Soft ground
: Hard ground 図１５．予測されるべき地山モデル 弾性係数 地質状況 E(kN/m2) E1 2.45 × 105 健全 ( 0m - 110m ) E2 2.45 × 105 健全 (110m - 130m ) E3 2.45 × 104 軟弱 (130m - 150m ) E4 2.45 × 104 軟弱 (150m - 200m ) 表５．各層の弾性係数の目的値 図１６は、各層の弾性係数の繰り返し計算による履歴を示している。各層の弾性係数は完全には目的の値に収束 していない。しかし、領域１と比較するとわかるように、切羽前方の断層の弾性係数が初期値から大きく下がり、 軟弱断層の目的値に近づいている。つまり、領域２でのトンネル内空変位データが軟弱断層に影響されており、そ の反応が十分に考慮できたといえる。つまり、計測値に軟弱断層の影響を含んでいれば、必ず切羽前方の軟弱断 層を予測することができるであろう。
0 50000 100000 150000 200000 250000 300000 350000 0 10 20 30 40 50 60 70 80 90 100 Youngs modulus [kN/m2] Iteration 0m - 110m ( E1 ) 110m - 130m ( E2 ) 130m - 150m ( E3 ) 150m - 200m ( E4 ) 図１６．領域２での各層の弾性係数の履歴
10
結論
この報告書では、トンネル切羽前方の断層の弾性係数を予測する新しい手法を提案した。効果的な数値解析手法 として、３次元有限要素モデルを用いた弾性係数の同定解析が検証された。そして、トンネル掘削表面上の計測 データを用いて、切羽前方の地質状況を予測できることが明らかになった。この予測方法で最も重要な点は、正 確な計測データを得る方法である。そのために、解析モデルではトンネル支保工を考慮した掘削解析が行われた。 計測点としては、１断面のみの場所が適用された。実際のトンネル現場では、多くの計測点を設けることが困難 であるため、１断面の計測点のみで解析を行った。今回の解析例では、切羽前方の断層が切羽面に対して、垂直 に仮定されている。今後は、斜めに存在した軟弱断層を考慮し、どのトンネル現場でも適用できる予測システム を開発し、実際のトンネルで適用されることが望まれる。そのために、現在、数箇所のトンネル現場で計測され た変位データを用いて、この手法の改良を行っている。References
[1] Sakurai, H. and Tanigawa, M., Study of some problems on the back analysis of measured displacements in tunnelling, Vol.1, Journal of Geotechnical Engineering, pp85-94(1990).
[2] Sugimoto, M., Sakurazawa, M. and Kageyama, S., Three dimensional reverse analysis on ground properties with time-dependence, 9th International Congress on Rock Mechanics, Vol.2, pp1405-1408(1999).
[3] Y.M. Hsieh and A.J. Whittle, A computational strategy for solving three-dimensional tunnel excavation problems, Second M.I.T. Conference on Computational Fluid and Solid Mechanics (2003).
[4] Nojima, K. and Kawahara, M., Mesh Generation of Three-dimensional Underground Tunnels Based on the Three-Dimensional Delaunay Tetrahedration, Vol.5, Journal of Applied Mechanics JSCE, pp253-262(2002). [5] R.Fletcher and C.M.Reeves., Function Minimization by Conjugate Gradients, Computer J.,
