論文

地盤構造の安定性評価に用いるヤング率の同定解析に関する基礎的研究
倉橋 貴彦1・吉原 健太2・井山 徹郎3
1 長岡技術科学大学大学院

機械創造工学専攻（Department of Mechanical Engineering, Graduate School of Nagaoka University of Technology）
2 長岡技術科学大学大学院

機械創造工学専攻（現：2 株式会社メイテック）（Department of Mechanical Engineering, Graduate School of Nagaoka University of Technology）
3 機械工学科（Department of Mechanical Engineering, National Institute of Technology, Nagaoka College）

Fundamental Study on Identification Analysis of The Young's Modulus
Used in Stability Evaluation of Ground Structure
Takahiko KURAHASHI1, Kenta YOSHIHARA2 and Tetsuro IYAMA3Abstract
In this paper, we present results of the fundamental study for identification analysis of the Young'smodulus used in stability evaluation of ground structure. In general, though the stability evaluation ofground structure is carried out by the shear strength reduction method, the value of the Young's modulus
used in this evaluation is estimated by the result of the standard penetration test. Therefore, it is necessary
to know the suitability of the Young's modulus, because the result of the shear strength reduction method
depends on the prescribed Young's modulus. As for the introduction to obtain the appropriate Young'smodulus, we investigate about the identification method of the Young's modulus based on the adjoint
variable method.
Key Words : identification, finite element method, adjoint variable method, ground structure

１．はじめに般にヤング率は，標準貫入試験によって得られたN値を用いる等して算定されるものであり，実際に斜地盤の安定性評価として斜面の安全率等を計算

面の安全率を算定する際に適切な値であるとは限らする方法の一つとしてせん断強度低減法 1)が使用さない．そのため，斜面全体が崩壊した際の変位情報れている．この方法は，斜面の強度定数を少しずつを用いて地盤の安定性評価に用いるヤング率を推定低減することにより，斜面全体が崩壊した時点で全する方法について検討を行う．本検討では，随伴変体安全率を定義する方法である．この解析を行う場

数法 2)を適用し，反復計算法によるヤング率の推定合，ヤング率も事前に解析条件として入力する．一解析を行う．地盤は弾完全塑性体として変形解析を6行う．地盤の支配方程式の離散化は有限要素法 3)に𝑒𝑒𝑒𝑒基づき行い，塑性変形の解析は修正Newton-Raphson法により行う4)．

２．ヤング率同定問題に対する定式化
２．１{Δ𝜎𝜎(𝑒𝑒) } = �𝐷𝐷(𝑒𝑒) � �Δ𝜖𝜖(𝑒𝑒) �(3)′�𝜎𝜎(𝑒𝑒)� = �σ(𝑒𝑒) � + {Δσ(𝑒𝑒) }(4)∗�𝜎𝜎(𝑒𝑒)� = �1 + 𝛽𝛽(𝑒𝑒) � ��σ′(𝑒𝑒) � − �𝜎𝜎(𝑒𝑒) �� + �𝜎𝜎(𝑒𝑒) �𝑚𝑚

弾塑性変形の解析

𝑚𝑚(5)弾塑性変形の解析を行うにあたり，まず線形弾性

また，上記の右辺の残差力ベクトル{Δf}が零ベクト体の支配方程式を離散化した式を式(1)のように書ルになり，{Δu}が零になるまで上記の計算を続けき表す．[K]は剛性行列，{u}は変位ベクトル，{f}はることになるが，斜面の安全率を計算するためには，修正Newton-Raphson法による計算の終了後，地盤の外力および反力によるベクトルを示す．

せん断強度を低減させ，上記の計算を繰り返すこと(1)[𝐾𝐾]{𝑢𝑢} = {𝑓𝑓}になる．安全率Fは式(6)により与えられ，土のせん断強度τr（滑りに抵抗する力）と滑ろうとする力τs塑性変形の場合は，外力と変位の関係が線形の関係の比により表される．土のせん断強度τrは拘束圧σで表せないことから，反復計算による変位の算定をとの関係から式(7)のように与えられる．式(7)にお行う．本検討では修正 Newton-Raphson 法を導入し，いてcは粘着力， φ はせん断抵抗角（内部摩擦角）剛性行列を固定した上で反復計算を行う．各反復過を示す．式(7)は式(8)のように変形することができ，程では，残差力ベクトル{Δf}を外力・反力のベクト

土のせん断強度τr を少しずつ下げていく（安全率Fル{f}から差し引き，各状態での変位増分ベクトルの値を少しずつ上げていく）ことで，最終的に安全{Δu}を足し合わせることにより，変位場を求める

率Fを求める．具体的な解析の中での操作としては，ことになる．

土のせん断強度τrを下げる場合は，せん断抵抗角，ダイレタンシー角，粘着力を低減し，弾塑性解析を(2)[𝐾𝐾]{Δ𝑢𝑢} = {𝑓𝑓} − {Δ𝑓𝑓}行う．計算プログラム詳細は文献 5)を参照して頂くことにする．また，せん断強度低減法では，ある安

まず，解析の最初のステップでは，残差力ベクトル{Δf}は零ベクトルとし，変位増分ベクトル{Δu}が算全率Fの際に修正Newton-Raphson法よる計算が収束せず，ある程度の反復回数を続けても計算が終了し定される．その変位ベクトルに対して個々の要素においてひずみ増分ベクトル{Δε(e)}を算定し，応力増ない場合に斜面崩壊に至ったとして，その時のFの分ベクトル{Δσ(e)}を算定する．ここに[Dep(e)]は弾塑

では500回を修正Newton-Raphson法における繰り返値を斜面安全率と定義している．本論文に示す解析

4)

し回数の上限値としている．

性状態における弾性係数行列であり，文献 を参照して頂くこととする．式(4)に示すように，現状での応力場に応力増分ベクトルを加えることで，塑性変形時における応力場{σ’(e)}が求まる．しかし，こτ𝐹𝐹 = r(6)𝜏𝜏𝑟𝑟 = 𝑐𝑐 + 𝜎𝜎 tan 𝜙𝜙(7)𝜏𝜏𝑠𝑠 = r(8)𝜏𝜏𝑠𝑠の応力場が計算ステップごとに降伏条件式から外れていないか修正Newton-Raphson法における各ステップで補正をかけながら計算を行う必要がある．降伏条件式から外れている場合，降伏曲面上に戻すための 処理 として ，式 (5)に 示す 応力 の補正 を行 う．

τ𝐹𝐹{σ(e)m}は平均主応力であり，β(e)は補正係数を示す．
ここで求まった補正後の応力{σ*(e)}を用いて，個々以下に弾塑性解析の計算手順を示す．の要素における体積力を計算することで，残差ベク
1.

トル{Δf}を求めて，式(2)の右辺のにように外力ベクトルの修正を行う．その後でまた変位増分ベクトル{Δu}を求めて，繰り返し計算が続くことになる．

安全率F0および，弾塑性解析における計算条件を入力する．

72.

弾性解析（式(1)）を行う．

3.

弾性解析により得られた変位場に対して塑性ひずみを計算し，塑性時の応力場を算定する（式(3)，式(4)）．
4.

解析では最急降下法によるヤング率の更新を行い，その際に使用するステップ長さは1としている．

式(3)によって求まった応力を式(5)により補正[𝐾𝐾]𝑇𝑇 {𝜆𝜆} = −[𝑄𝑄]𝑇𝑇 {𝑢𝑢 − 𝑢𝑢𝑜𝑜𝑜𝑜𝑜𝑜. }し，式(6)により残差力ベクトルを算定する．
5.

式(2)の右辺の荷重ベクトルを残差力ベクトルに6.

より修正し，変位増分ベクトル{Δu}を算定する．
変位場を{u}n+1={u}n +{Δu}nにより更新する．

7.

修正Newton-Raphson法の繰り返し計算の回数が𝜕𝜕𝐽𝐽∗𝜕𝜕𝜕𝜕= {𝜆𝜆}𝑇𝑇𝜕𝜕[𝐾𝐾]
𝜕𝜕𝜕𝜕{𝑢𝑢}

(11)

(12)

ヤング率の同定計算のアルゴリズムを以下に示す．

上限値を上回った場合には，弾塑性解析を終了
し，その際算定された安全率を斜面の安全率とする．上記でない場合は，次のステップへ進む．
8.

1.

修正Newton-Raphson法の上限回数を上回っている．

ない場合は次の収束判定を行う．収束判定
Σ|{Δu}n|/|{u}n+1|≦εを満たした場合は次のステッ

2.

プへ進み，満たさない場合はステップ2へ戻る．

3.

弾塑性解析により得られた変位場の値を用いて4.

収束判定を行う．|J(l+1)- J(l)|<εなら計算を終了．

与えられたヤング率を使用して弾塑性解析（2.1
節参照）を行う．
評価関数（式(9)）の計算を行う．

ここにΣは全節点での足し合わせを示す．
9.

ヤング率の初期値E(0) ，収束判定定数を設定す

安全率Fk+1=Fk+ΔFにより更新する．せん断抵抗角，ダイレタンシー角，粘着力を安全率Fk+1にそうでなければ，次のステップへ進む．
5.

より除し，ステップ2へ戻る．

弾塑性解析で得られた変位場を用いて随伴方程式（式(11)）の解析を行う．（本解析では，随２．２

伴方程式の計算も，2.1節に示す弾塑性解析の計

ヤング率の同定に対する定式化

本論文では，地盤内の変位量が観測されている状算の流れと同様に計算を行うことにしている．

態を想定し，その状態におけるヤング率を同定する

剛性行列が対称であることから，外力ベクトル問題を考える．まず，式(9)に示す評価関数を定義を式(11)のものに変更し，変位ベクトル{u}を随する．uは弾塑性解析によって算定された変位，uobs.

伴変数{λ}とすることで解析を実施する．）

は観測変位，[Q]は重み対角行列を示す．ここに目6.

的は式(9)の評価関数Jを最小とするようなヤング率Eを同定することである．

ラグランジュ関数のヤング率に対する勾配（式(12)）を計算する．

7.

ヤング率を更新し，ステップ2へ戻る．ヤング率の更新はE(l+1)=E(l)－∂J*(l)/∂Eにより行う．

1

𝐽𝐽 = {𝑢𝑢 − 𝑢𝑢𝑜𝑜𝑜𝑜𝑜𝑜. }𝑇𝑇 [𝑄𝑄]{𝑢𝑢 − 𝑢𝑢𝑜𝑜𝑜𝑜𝑜𝑜. }2(9)３．数値実験

式(10)の変位は，式(1)の結果より計算されることか

解析モデルおよび境界条件を図-1 に，また計がある．そのため，式(1)を考慮した式(9)の最小化
算条件の表を図-1 に示す．六面体二次要素によ問題を解くために，随伴変数ベクトル{λ}を導入し， る要素分割を行い，弾塑性解析における数値式(9)を式(10)のように拡張する．式(10)はラグラン積分に対して，ガウス・ルジャンドル法を適ジュ関数または拡張評価関数と呼ばれる．
用する．外力は自重のみ考慮する．本論文では，ヤング率が 10,000kPa の場合の変位分布を∗𝑇𝑇{𝜆𝜆}{𝑓𝑓})(10)([𝐾𝐾]{𝑢𝑢} −𝐽𝐽 = 𝐽𝐽 +弾塑性解析により求め，鉛直方向の変位を解
析モデル上の全節点において測定した値と仮式(10)の停留条件を誘導するために，式(10)の第一定する．ヤング率を 1,000kPa をヤング率の同定変分を計算すると，各変数の勾配から以下の式が誘解析の初期値として，正解とする 10,000kPa が導される．式(10)の変位ベクトル{u}に対する勾配算定されるか検証を行う．
から式(11)に示す随伴変数に対する方程式（随伴方上記の条件のもと同定解析を実施した結果程式）が得られる．また，同定するヤング率Eに対を以下に示す．図-2 に 0 イタレーション目におする勾配として式(12)が得られる．本論文におけるける鉛直方向の変位分布を示す．ヤング率をら，式(9)の制約条件式として式(1)を考慮する必要

81,000kPa と与えた場合にも適切に変位分布が得られていることがわかる．ヤング率の正解値を 10,000kPa としているため，変位分布に違いがあることから，鉛直方向の変位分布の差を外力とした随伴方程式を解くことで随伴変位を 得 る こ と と な る ．本解 析 で は ， 状 態 変数（x,y,z 方向の変位）の求解に用いた弾塑性解析と同様の計算の流れで随伴方程式の解析を行った．図-3 に 0 イタレーション目における鉛直方向の随伴変位分布を示す．状態変数の結果と同様に，発散することなく随伴変数が適切に得られていることがわかる．この結果を用いてラグランジュ関数のヤング率に対する勾配を計算し，ヤング率の同定を行った．0 イタ
レーション目における評価関数の値で正規化した評価関数の履歴を図-4 に示す．前後の評価関数の差の絶対値が 0.01を下回った場合に計算を終了するようにした結果，反復回数 19 回にて計算を終了した．ヤング率の履歴を図-5 に示す ． 反 復 回 数 19 回 目 に お け る ヤ ン グ 率 は1,446.27kPa であり，解析初期の 1,000kPa に比べ，正解とする 10,000kPa に近づく結果となったが，評価関数の収束判定を厳しくするあるいは，最急降下法におけるステップ長さを適切に設定する等，同定解析における工夫が必要であると考えられる．

ABEF面は完全にフリーE

上面ADHE面は完全にフリー

表-1 解析条件

節点数
要素数
数値積分におけるガウス点の点数
せん断抵抗角ダイレタンシー角
ヤング率
ポアソン比
粘着力
単位体積重量収束判定定数 ε406
（6 面体 2 次要素）
60
810.0 deg.
10.0 deg.
初期値 1,000 kPa（正解値 10,000 kPa）
0.30
20.0 kN/m2
20.0 kN/m30.010

奥行面EFGHは対称境界（x方向のみ変位0)

図-2 鉛直方向の変位分布（0 イタレーション時）

5mAHD4mF

BG

手前面ABCDは（x,y方向 変位0)C

2.5m
（対称モデルのため半分）

底面BCGF面x,y,z方向固定（x,y,z方向 変位0）

図-1 解析モデルおよび境界条件

図-3 鉛直方向の随伴変位分布（0 イタレーション時）

9初期の値で正規化した評価関数

1.20

謝辞：本論文の解析結果は，九州大学情報基盤研究

1.00

開発センターの高性能アプリケーションサーバを使0.80わせていただいた上で計算を行ったものである．セ

0.60

ンターのスタッフの皆様に対して，ここに謝意を表0.40

す．

0.20
0.00

05101520

参考文献

反復回数1) 鵜飼恵三, “弾塑性 FEM による斜面の全体安全率の計算法”，土質工学会論文報告集，Vol.29, No.2, pp.190-195,

図-4 正規化した評価関数の履歴

ヤング率,kPa

1989.
2) 吉原健太，倉橋貴彦，村上祐貴，外山茂浩，池田富士

1500
1450
1400
1350
1300
1250
12001150
11001050
1000

雄，井山 徹郎，井原 郁夫， “随伴変数法に基づくコンクリート構造内部における欠陥深度の推定解析 （打音検査による観測変位量を想定した数値実験）”, 日本機械学会論文集, 第85巻，第869号，pp.1-13, 2019.
3) 竹内則雄，寺田賢二郎，樫山和男，計算力学 : 有限要素法の基礎，森北出版，2003．
0510154) 地盤工学会地盤技術者のための FEM シリーズ編集委員

20反復回数会編, “地盤技術者のための FEM シリーズ 2 弾塑性有限要素法がわかる”, 地盤工学会, 2003.

図-5 ヤング率の履歴

5) 地盤工学会地盤技術者のための FEM シリーズ編集委員会編, “地盤技術者のための FEM シリーズ 3 弾塑性有限要素法をつかう”, 地盤工学会, 2003.

４．おわりに(2020. 8. 1 受付)

本論文では，斜面の安全性評価に使用される弾塑性解析におけるヤング率を変位場により同定する方法について数値実験による検討結果を整理した．ヤング率の同定の定式化に際しては，随伴変数法を導入し，弾塑性解析においては修正Newton-Raphson法による計算を行った．また随伴方程式の計算は，弾塑性解析と同様の計算手順により実施した．本論文において示した計算アルゴリズムによりヤング率の同定解析を行った結果，正解とするヤング率には近づく傾向が得られた．しかしながら，収束判定の基準によっては，正解値から離れた箇所で解が得られる場合や，ステップ長さの設定の仕方によっては同定計算における反復回数が膨大になる可能性もあるため，この点は今後の課題として考えている．本解析法が確立した場合，実際に斜面崩壊をした現場において変位の値を測定することで，適切なヤング率を推定することができる．実際の斜面崩壊時におけるヤング率を用いて斜面の安定性解析を行うことで，将来的には，より適切な安全率を把握することができると考えている．
10