1
モールの応力円から内部摩擦角、粘着力を求めるためのプログラム
益永八尋 Ⅰ.プログラムの考え方 土質試験結果からモールの応力円を描き、内部摩擦角と粘着力を求めるプログラムの開 発をおこなった。このプログラムを作成するに当って、どのような考え方をしているかを 以下に技術資料として作成する。 モールの応力円を作成するプログラム言語はVB とした。これは、Excel の VBA では描 画機能がなく、Excel モードのグラフ機能ではモールの応力円（半円）を描くことが困難で あるためである。 モールの応力円に包絡線を作成するには少々のテクニックが必要であるため、どのよう にすれば、包絡線を作図できるかの考え方を以下に記述する。 ①第1 の方法 手作業により包絡線を図のように作成する。 ②第2 の方法 自動で包絡線を作成する。 第１の方法は、プログラムコードの記述がなく、簡単に作成できる。第２の方法は、プ ログラムコードが複雑となり、確実に描けるのかも不明であるが、自動であるため、誰で もが簡単に内部摩擦角、粘着力を求めることが可能である。モールの円が2 個の場合には、 包絡線は 1 本となるので問題ないが、モールの円が 3 個以上の場合には、包絡線が考え方 の違いにより 2 本以上となる場合がある。このため、このような場合の包絡線の作成の仕 方について明確な基準（考え方）を持つ必要がある。 包絡線の傾きは、回帰解析を行い求める。 モールの応力円を 3 以上描ける場合には、回帰分析をおこない、傾きを決定する。傾き を求めるデータは、モール円の中心X 座標値（側圧）における Y 座標値（モール円の頂点） を回帰解析のデータとする。この時の解析結果を図 1 に示す。包絡線の傾きと回帰式の傾 き（一次式の傾き）は一致するので、包絡線の本数は 1 本にすることができる。また、回 帰式を求めることから、個人差はなくなり、合理的な内部摩擦角を求めることができる。 粘着力は最小となるものを採用する 次に、粘着力を求めることになるが、これは簡単ではない。このプログラムでは、以下 のように考えて包絡線を決定する。土の安定や、構造物の設計において粘着力を大きく採 れば、一般的には危険側の設計となる。このため、土や構造物の設計としては安全性を考 え、粘着力が小さくなるように包絡線を決定するのが最善の方法と考える。このことから、 包絡線の傾き決定したときのモールの応力円と傾きを用いて、粘着力を求める。この粘着 力が最小となる包絡線を採用する。三軸圧縮試験結果の回帰解析の結果である相関係数Ｒ2 は一般的には 1 にならない。このため、粘着力を決めるためにはモールの応力円の個数分2 の包絡線を求めることになる。包絡線は一次式で表せるのでこのときの係数（切片）を求 めればよいことになる。この係数（切片）が粘着力となる。包絡線はモールの応力円に外 接する直線であるため、包絡線の式は下記三式を解くことにより求めることができる。 包絡線の式 Ｙ＝A1・X + B1 --- ① 包絡線に直交する線の式 Ｙ＝A2・X + B2 --- ② モールの円の式 ( X － X0)2 + Y2 = R2 --- ③ A1：回帰解析結果より得られる（既知） A2：－1/A1 （既知） X0：モールの円の中心X 座標（既知） R：モールの円の半径（既知） 包絡線に直交する線はモールの円の中心を通るので、円の中心座標X=X0、Y=0 を式②に 代入することにより②式の係数B2を決定できる。このことにより②式を③式に代入すれば、 X に関する二次式を得ることができる。これを示せば B2=－A2・X0 （1+A22）・X2+（2・A2・B2－2・X0）・X+（X02+B22－R2）＝0 従って、二次方程式のすべての係数を求めることが可能となり、その解は二次方程式の 解を求める公式を用いて得られる。二次方程式の解は下記の条件である。このことは図 2 より明らかである。 X≧0 and X≦X0 ここに、 a=1+A22 b=2・A 2・B2－2・X0 c= X02+B22－R2 上記条件を満足するX を求め、この X を②式に代入すれば Y が得られる。得られた X と Y を①式に代入すれば、B2が得られる。 以上の計算をモールの応力円の個数分繰り返し、得られたB2が最小となるものが粘着力ｃ となる。
3 Ⅱ.データ入力と結果 例題5-13 ｐ.147 砂30％、シルト 20％、粘土 50％の粒度組成をもつ粘性地盤上に構造物を築 造することになった。地盤の支持力を計算するためには土の粘着力と内部摩 擦角が必要であるから、三軸圧縮試験を実施した。その試験結果は下記の 通りである。土の内部摩擦角と粘着力を求めよ。 試験結果のデータからモールの応力円を描き、包絡線を求める。この包絡線の 傾斜角が内部摩擦角、縦軸との交点が粘着力となる。 内部摩擦角 φ = 14.1(°) 粘着力 c = 0.146(kg/cm2) データ入力上の注意点 注意点1 三軸圧縮試験のグラフにおいて近似直線式の相関係数R2 が小さい場合には、測定データを 適宜削除し相関係数R2 が 1 に近くなるように修正してください。 測定データの削除は、試験データを残す必要があるため、グラフ作成するデータとして有 効であれば”データの有効性”の列（H 列）を”有効”にし、無効であれば H 列を”無効” にする。 注意点2 データ入力可能な行数は 10 行までである。一般的にはある地点の土質試験数は 10 以下 であれば十分であるとし、この設定を行っている。これ以上の場合にはプログラム変更で 対応可能である。 データ入力 測定番号 含水比 湿潤密度 側圧 圧縮強さ 破壊ヒズミ 番号 ｗ(%) γ t(g/cm3 ) σ3(kg/cm2) σ1(kg/cm2) ε 0(%) X座標 Y座標 半径 1 147.5 1.28 0.20 0.502 5.5 有効 1 0.451 0.000 0.251 2 152.6 1.29 0.40 0.643 1.0 有効 2 0.722 0.000 0.322 3 150.2 1.30 0.60 0.770 1.5 有効 3 0.985 0.000 0.385 4 148.7 1.31 0.80 0.635 5.5 無効 中心座標 データの 有効性
4 y = 0.251x + 0.1387 R² = 0.9995 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0 0.2 0.4 0.6 0.8 1 1.2 圧 縮 強さ (kg /c m 2) 側圧(kg/cm2)
三軸圧縮試験
系列1 線形 (系列1) 図2 モールの応力円 図1 三軸圧縮試験結果5 Ⅲ.プログラムコード
Ⅲ-1.Excel（VBA）のコード
Private Sub CommandButton1_Click() Dim Siguma3(10) As Single
Dim Siguma1(10) As Single
ActiveWorkbook.Sheets("Sheet1").Select '--- For I = 1 To 10 AA = Cells(I + 17, 8) Select Case AA Case "有効" J = J + 1 Siguma3(J) = Cells(I + 17, 5) Siguma1(J) = Cells(I + 17, 6) Dn = J Case "無効" Case "" Exit For End Select Next I '--- For J = 1 To Dn Cells(J + 17, 9) = J
Cells(J + 17, 10) = Siguma3(J) + Siguma1(J) / 2 Cells(J + 17, 11) = 0 Cells(J + 17, 12) = Siguma1(J) / 2 Next J '--- Call Module1.Main End End Sub
6
Public Declare Function GetExitCodeProcess Lib "kernel32" _ (ByVal hProcess As Long, lpExitCode As Long) As Long Public Declare Function OpenProcess Lib "kernel32" _
(ByVal dwDesiredAccess As Long, ByVal bInheritHandle As Long, _ ByVal dwProcessID As Long) As Long
Public Const PROCESS_QUERY_INFORMATION = &H400 '---
Public MyEXE As String Public MyEXE_Fie As String Public MyPath As String '--- Public a As Single Public X0(10) As Single Public Y0(10) As Single Public R(10) As Single Public Dn As Integer Public c As Single Public Fai As Single
Sub EXE_RUN()
Dim dwProcessID As Long Dim hProcess As Long Dim lpdwExitCode As Long Dim ret As Long
'--- MsgBox MyPath
dwProcessID = Shell(MyEXE, 1)
hProcess = OpenProcess(PROCESS_QUERY_INFORMATION, True, dwProcessID)
Do
ret = GetExitCodeProcess(hProcess, lpdwExitCode) DoEvents
Loop While lpdwExitCode
MsgBox MyEXE_Fie & "終了しました。" End Sub
7 Sub Main() 'データ取得 Sheets("Sheet1").Select For I = 1 To 10 AA = Cells(I + 17, 9) If AA <> "" Then X0(I) = Cells(I + 17, 10) Y0(I) = Cells(I + 17, 11) R(I) = Cells(I + 17, 12) Else Dn = I - 1 Exit For End If Next I '--- a = Cells(17, 15) '--- Call Module1.Data_Print '--- MyPath = ThisWorkbook.Path
MyEXE = MyPath & "¥" & "Moll_Circle.exe" '--- Call EXE_RUN '--- Call Module1.Data_Input_Fai '--- Sheets("Sheet1").Select
Cells(32, 2) = "内部摩擦角 φ = " & Format(Fai, "#0.0") & "(°)" Cells(33, 2) = "粘着力 c = " & Format(c, "#0.000") & "(kg/cm2)"
End Sub
Sub Data_Print()
Dim Dummy As String
8 'MyPath = "H:¥技術計算 EXCEL¥土質" MyFile = "モールの応力円データ.txt"
MyPath = ThisWorkbook.Path & "¥" & MyFile Open MyPath For Output As #1
Write #1, "'中心座標" Write #1, MyPath
Write #1, "'X 座標 Y 座標 半径" Write #1, Dn
For I = 1 To Dn
Write #1, I, X0(I), Y0(I), R(I) Next I '--- Write #1, a Close #1 End Sub
Public Function Houteishiki_Kai1(a As Single, b As Single, c As Single) As Single Houteishiki_Kai1 = (-b + Sqr(b * b - 4 * a * c)) / (2 * a)
End Function
Public Function Houteishiki_Kai2(a As Single, b As Single, c As Single) As Single Houteishiki_Kai2 = (-b - Sqr(b * b - 4 * a * c)) / (2 * a) End Function Sub Test() X1 = Houteishiki_Kai1(1, -2, 1) X2 = Houteishiki_Kai2(1, -2, 1) End End Sub Sub Data_Input_Fai() MyFile = "内部摩擦角_粘着力.txt" MyPath = ThisWorkbook.Path
Open MyPath & "¥" & MyFile For Input As #1 Input #1, Fai
9 Close #1
End Sub
Ⅲ-2 VB のコード 1.Form1
Private Sub Command1_Click() End
End Sub
Private Sub Command2_Click() Dim X(10) As Single
Dim Y(10) As Single 'モールの応力円作図 '--- Form1.Picture1.Width = 13000 Form1.Picture1.Height = 7000 Form1.Picture1.Line (500, 6000)-(10000, 6000), RGB(255, 0, 0) Form1.Picture1.Line (500, 1000)-(500, 6000), RGB(255, 0, 0) '--- 'X0(I)+R(I)の最大値を求める Dat_Max = -9999 For I = 1 To Dn
If X0(I) + R(I) >= Dat_Max Then Dat_Max = X0(I) + R(I) End If Next I '--- K = Int(Form1.Picture1.Height / Dat_Max) '--- For I = 1 To Dn 'モールの応力円（半円）を描く
Form1.Picture1.Circle (X0(I) * K + 500, 6000), R(I) * K, RGB(255, 0, 0), 0, Pai '---
10 'X 軸目盛
Form1.Picture1.Line (X0(I) * K + 500, 6000)-(X0(I) * K + 500, 6100), RGB(255, 0, 0) Form1.Picture1.CurrentX = X0(I) * K + 500 Form1.Picture1.CurrentY = 6000 Form1.Picture1.Print X0(I) 'Y 軸目盛 YM = 6000 - R(I) * K
Form1.Picture1.Line (400, YM)-(500, YM), RGB(255, 0, 0) '--- Form1.Picture1.CurrentX = 100 Form1.Picture1.CurrentY = YM Form1.Picture1.Print R(I) Next I '--- 'Picture1 に文字（Text)を出力する Form1.Picture1.CurrentX = 5000 Form1.Picture1.CurrentY = 1000 Form1.Picture1.FontSize = 18 Form1.Picture1.Print "モールの応力円" Form1.Picture1.FontSize = 10 '--- '包絡線の式を求める '①最小 2 乗法により回帰曲線の傾きを求める B1_Min = 9999 '--- For I = 1 To Dn '包絡線の切片の値（粘着力）を求める a2 = -1 / a1 b2 = -a2 * X0(I) a = (1 + a2 * a2) b = -2 * X0(I) + 2 * a2 * b2 c = (X0(I) ^ 2 + b2 ^ 2 - R(I) ^ 2)
11 '---
K1 = Houteishiki_Kai1(a, b, c) 'K1>=0 かつ X0(I) K2 = Houteishiki_Kai2(a, b, c) 'K2>=0 かつ X0(I) '---
If K1 >= 0 And K1 <= X0(I) Then X1(I) = K1
Else
If K2 >= 0 And K2 <= X0(I) Then X1(I) = K2
End If End If
'---
Y1(I) = a2 * X1(I) + b2 'Y1(I)>=0 '--- b1 = Y1(I) - a1 * X1(I) '--- If b1 <= B1_Min Then B1_Min = b1 End If Next I '--- '包絡線の作図 b1 = B1_Min X(1) = 0 X(2) = X0(Dn) + R(Dn) Y(1) = b1 Y(2) = a1 * X(2) + b1 '---
Form1.Picture1.Line (500 + X(1), 6000 - Y(1) * K)-(500 + X(2) * K, 6000 - Y(2) * K), RGB(255, 0, 0)
'---
Form1.Picture1.CurrentX = 1000
Form1.Picture1.CurrentY = 6000 - Y(2) * K
Form1.Picture1.Print "包絡線 Y = " & Format(a1, "#0.000") & "・X + " & Format(b1, "##0.000")
12 '---
Fai = Atn(a1) * 180 / Pai
MsgBox "粘着力 c = " & B1_Min Call Data_Print(Fai, B1_Min) '---
End Sub
Private Sub Form_Initialize() Dim Dummy As String
'---
'MyPath = "H:¥技術計算 EXCEL¥土質" MyFile = "モールの応力円データ.txt" MyPath = App.Path
Open MyPath & "¥" & MyFile For Input As #1 Input #1, Dummy
Input #1, MyPath Input #1, Dummy Input #1, Dn For I = 1 To Dn
Input #1, Dummy, X0(I), Y0(I), R(I) Next I
Input #1, a1 Close #1 End Sub
Sub Data_Print(Fai As Single, c As Single) MyFile = "内部摩擦角_粘着力.txt" MyPath = App.Path
Open MyPath & "¥" & MyFile For Output As #1 Write #1, Fai
Write #1, c Close #1 End Sub
13 2.標準モジュール
Public X0(10) As Single Public Y0(10) As Single Public R(10) As Single Public Dn As Single
Public Const Pai As Single = 3.1415926 Public Dat_Max As Single
Public X1(10) As Single Public Y1(10) As Single Public a1 As Single Public a2 As Single Public b1 As Single Public b2 As Single Public a As Single Public b As Single Public c As Single Public B1_Min As Single
Public Fai As Single '内部摩擦角
Public Function Houteishiki_Kai1(a As Single, b As Single, c As Single) As Single Houteishiki_Kai1 = (-b + Sqr(b * b - 4 * a * c)) / (2 * a)
End Function
Public Function Houteishiki_Kai2(a As Single, b As Single, c As Single) As Single Houteishiki_Kai2 = (-b - Sqr(b * b - 4 * a * c)) / (2 * a)
