モールの応力円から内部摩擦角 粘着力を求めるためのプログラム 益永八尋 Ⅰ. プログラムの考え方土質試験結果からモールの応力円を描き 内部摩擦角と粘着力を求めるプログラムの開発をおこなった このプログラムを作成するに当って どのような考え方をしているかを以下に技術資料として作成する モールの応力円を作成するプログラム言語は VB とした これは Excel の VBA では描画機能がなく Excel モードのグラフ機能ではモールの応力円 ( 半円 ) を描くことが困難であるためである モールの応力円に包絡線を作成するには少々のテクニックが必要であるため どのようにすれば 包絡線を作図できるかの考え方を以下に記述する 1 第 1 の方法手作業により包絡線を図のように作成する 2 第 2 の方法自動で包絡線を作成する 第 1の方法は プログラムコードの記述がなく 簡単に作成できる 第 2の方法は プログラムコードが複雑となり 確実に描けるのかも不明であるが 自動であるため 誰でもが簡単に内部摩擦角 粘着力を求めることが可能である モールの円が 2 個の場合には 包絡線は 1 本となるので問題ないが モールの円が 3 個以上の場合には 包絡線が考え方の違いにより 2 本以上となる場合がある このため このような場合の包絡線の作成の仕方について明確な基準 ( 考え方 ) を持つ必要がある 包絡線の傾きは 回帰解析を行い求める モールの応力円を 3 以上描ける場合には 回帰分析をおこない 傾きを決定する 傾きを求めるデータは モール円の中心 X 座標値 ( 側圧 ) における Y 座標値 ( モール円の頂点 ) を回帰解析のデータとする この時の解析結果を図 1 に示す 包絡線の傾きと回帰式の傾き ( 一次式の傾き ) は一致するので 包絡線の本数は 1 本にすることができる また 回帰式を求めることから 個人差はなくなり 合理的な内部摩擦角を求めることができる 粘着力は最小となるものを採用する次に 粘着力を求めることになるが これは簡単ではない このプログラムでは 以下のように考えて包絡線を決定する 土の安定や 構造物の設計において粘着力を大きく採れば 一般的には危険側の設計となる このため 土や構造物の設計としては安全性を考え 粘着力が小さくなるように包絡線を決定するのが最善の方法と考える このことから 包絡線の傾き決定したときのモールの応力円と傾きを用いて 粘着力を求める この粘着力が最小となる包絡線を採用する 三軸圧縮試験結果の回帰解析の結果である相関係数 R 2 は一般的には 1 にならない このため 粘着力を決めるためにはモールの応力円の個数分 1
の包絡線を求めることになる 包絡線は一次式で表せるのでこのときの係数 ( 切片 ) を求 めればよいことになる この係数 ( 切片 ) が粘着力となる 包絡線はモールの応力円に外 接する直線であるため 包絡線の式は下記三式を解くことにより求めることができる 包絡線の式 Y=A1 X + B1 ---------------------- 1 包絡線に直交する線の式 Y=A2 X + B2 ---------------------- 2 モールの円の式 ( X - X0) 2 + Y 2 = R 2 --------------- 3 A1: 回帰解析結果より得られる ( 既知 ) A2:-1/A1 ( 既知 ) X0: モールの円の中心 X 座標 ( 既知 ) R: モールの円の半径 ( 既知 ) 包絡線に直交する線はモールの円の中心を通るので 円の中心座標 X=X0 Y=0 を式 2に代入することにより2 式の係数 B2 を決定できる このことにより2 式を3 式に代入すれば X に関する二次式を得ることができる これを示せば B2=-A2 X0 (1+A2 2 ) X 2 +(2 A2 B2-2 X0) X+(X0 2 +B2 2 -R 2 )=0 従って 二次方程式のすべての係数を求めることが可能となり その解は二次方程式の 解を求める公式を用いて得られる 二次方程式の解は下記の条件である このことは図 2 より明らかである X 0 and X X0 ここに a=1+a2 2 b=2 A2 B2-2 X0 c= X0 2 +B2 2 -R 2 上記条件を満足する X を求め この X を2 式に代入すれば Y が得られる 得られた X と Y を1 式に代入すれば B2 が得られる 以上の計算をモールの応力円の個数分繰り返し 得られた B2 が最小となるものが粘着力 c となる 2
Ⅱ. データ入力と結果例題 5-13 p.147 砂 30% シルト 20% 粘土 50% の粒度組成をもつ粘性地盤上に構造物を築造することになった 地盤の支持力を計算するためには土の粘着力と内部摩擦角が必要であるから 三軸圧縮試験を実施した その試験結果は下記の通りである 土の内部摩擦角と粘着力を求めよ データ入力測定番号 含水比 湿潤密度 側圧 圧縮強さ 破壊ヒズミ データの 番号 中心座標 w(%) γ t(g/cm 3 ) σ 3 (kg/cm 2 ) σ 1 (kg/cm 2 ) ε 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 無効 試験結果のデータからモールの応力円を描き 包絡線を求める この包絡線の 傾斜角が内部摩擦角 縦軸との交点が粘着力となる 内部摩擦角 φ = 14.1( ) 粘着力 c = 0.146(kg/cm2) データ入力上の注意点注意点 1 三軸圧縮試験のグラフにおいて近似直線式の相関係数 R2 が小さい場合には 測定データを適宜削除し相関係数 R2 が 1 に近くなるように修正してください 測定データの削除は 試験データを残す必要があるため グラフ作成するデータとして有効であれば データの有効性 の列 (H 列 ) を 有効 にし 無効であれば H 列を 無効 にする 注意点 2 データ入力可能な行数は 10 行までである 一般的にはある地点の土質試験数は 10 以下であれば十分であるとし この設定を行っている これ以上の場合にはプログラム変更で対応可能である 3
圧縮強さ (kg/cm2) 2015 年 3 月 23 日 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 三軸圧縮試験 0 0.2 0.4 0.6 0.8 1 1.2 側圧 (kg/cm2) y = 0.251x + 0.1387 R² = 0.9995 系列 1 線形 ( 系列 1) 図 1 三軸圧縮試験結果 図 2 モールの応力円 4
Ⅲ. プログラムコード Ⅲ-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 '----------------------- 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 5
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 & " 終了しました " 6
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 '-------------------------------------- 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)" Sub Data_Print() Dim Dummy As String '------------------------------------------- 7
'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) '--------------------- Write #1, a Close #1 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 Sub Data_Input_Fai() MyFile = " 内部摩擦角 _ 粘着力.txt" MyPath = ThisWorkbook.Path Open MyPath & " " & MyFile For Input As #1 Input #1, Fai Input #1, c 8
Close #1 Ⅲ-2 VB のコード 1.Form1 Private Sub Command1_Click() End 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 '--------------------------------- 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 '------------------------------------------ ' 目盛線作成 9
'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) '----------------------------------- 'Picture1 に文字 (Text) を出力する Form1.Picture1.CurrentX = 5000 Form1.Picture1.CurrentY = 1000 Form1.Picture1.FontSize = 18 Form1.Picture1.Print " モールの応力円 " Form1.Picture1.FontSize = 10 '----------------------------------- ' 包絡線の式を求める '1 最小 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) 10
'----------------------------------------- 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 '------------------ ' 包絡線の作図 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") 11
'----------------------------- Fai = Atn(a1) * 180 / Pai MsgBox " 粘着力 c = " & B1_Min Call Data_Print(Fai, B1_Min) '---------------------------------- 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) Input #1, a1 Close #1 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 12
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) End Function 13