圆弧条分法土坡稳定分析的VB电算解法
──解题思路说明与讲解
云南省水利水电学校 张华庆
圆弧条分法土坡稳定分析的VB电算解法
云南省水利水电学校 张华庆
【摘要】目前,大多数版本的土力学教材在讲到均质粘性土边坡稳定性分析方法时,均采用圆弧条分法,该法存在大量的作图及试算,若用人工求解,则不但会耗费大量宝贵的时间,而且还难以得到精确解。现就介绍作者自编的粘性土边坡简单圆弧条分法的VB电算解法。
【关键词】土坡稳定 圆弧条分法 VB电算 全自动 可视化
引言:迄今为止,求解均质粘性土边坡稳定最小安全系数Ksmin,通常都采用瑞典圆弧法或圆弧条分法,由于该法存在大量的手工作图以及滑弧试算,因此如果采用过去传统的人工方法求解,则因为效率低下,不但浪费稿纸和笔墨、浪费大量宝贵的时间,而且还难以得出较为精确的解答。为此,笔者利用Visual Basic 6.0编写和开发了可视化的土坡稳定分析程序,该程序及解法具有如下特点:①作图与试算全过程可视化跟踪显示,其优越性是过去的DOS计算程序所无法
比拟的;②速度快,能在一分钟之内(一般仅需10秒左右)找到最危险滑弧以及土坡稳定最小安全系数Ksmin;③充分体现圆弧条分法的精髓,精确度高;④简单易行,并且经过若干教材的例习题演算,证明结果一致;⑤演算过程进行速度可以调节,便于
1
观察和教学;⑥程序大小仅72KB,但却功能完善,完全可以满足教学需要(当然也完全可以用于工程计算);⑦自动化程度高:仅需输入必要的坡高、坡比(或坡度角),土的容重γ以及土的抗剪强度参数φ、c的值,即可自动完成计算,因此为全自动化设计,功能上明显优越于目前有的半自动化边坡稳定分析程序(如有的为DOS程序;有的虽然是可视化程序,但却是半自动化的,需要在AutoCAD里人为画出试算滑弧,然后才能让机器去计算土坡稳定最小安全系数,这自然影响时效及精度,因为你不可能有时间和耐心去画多达100条滑弧!——就别说更多的);⑧分条数可突破人为影响及限制,例如可以达到1000条之多;⑨无论何种条件下,当分条数逐渐增多时,土坡稳定最小安全系数Ksmin均能趋于一个稳定的极限值,证明计算方法以及计算结果是正确的且是可靠的和稳定的;⑩可以考虑渗透力的影响。
下面就本程序的编程思路作一些必要的说明和讲解,以便和有兴趣致力于此方面研究的同仁共同探讨,并希望得到专家的指点,以便使程序变得更加完美和人性化。
一、滑弧圆心的确定变得简单化:试算法首先要确定的是滑弧的圆心,由于过去是采用人工作图计算,工作量或者说劳动强度可想而知,于是必须尽可能地缩小滑动圆心的范围,这样一来才需要定出M点坐标(4.5H,2H)、确定O点所需的参数β
1
=25°,β2=35°,以及利用十字交叉方法确定滑动圆心Oi点的做法就是必需的和可以理解的;但是对于计算机来说,它最不怕的就是计算,因而这种做法就显得多余或者说不必要了,我们可以让它在平面上可能的区域(即所画圆弧至少有一端能交于坡.....面、坡顶或坡底,而另一端也不至于悬空的点)在纵横两个方向每间隔一定的点计算一次去寻找最危险滑弧,问题因此有了更为简单而圆满的解决办法。
二、分条面积计算简单化:由于分条宽度b是固定的,所以分条面积取决于分条高度hi,笔者采用解析几何方法加以处理:
2
①列出坡面及坡顶、坡底的方程,并且列出圆弧的方程,再列出第i个土条的竖线方程;
②用公式法或试算法(迭代原理)求出交点坐标,再用两点间距离公式计算出hi; ③分条面积计算公式:Si=hi*b (“*”号即“×”号);
④由于分条数可以很多,当分条数多到一定程度(一般当取b=R/100或更大)时滑体两端三角形面积部分对总面积的正负或加减影响可以不予考虑。
三、滑弧与坡面的交点分析:程序分别考虑了①滑弧与土坡上游交点的两种情况,即交于坡顶和交于上半坡面;②滑弧与土坡下游交点的两种情况,即交于坡底和交于下半坡面。这几种情况都是圆弧条分法所要求考虑的,均可由解析几何方法得出其交点坐标,然后由上下两交点之间的水平间距d及b=R/n计算出分条数fts。
四、滑弧半径及滑弧长的确定:在滑弧圆心确定后,滑弧半径通过两点间距离公式确定,即R(xOxA1)(yOyA1);至于滑弧长,可通过公式LR*计算滑弧
22总长,即在分条数fts确定后,水平方向正负两边圆心角即可用反三角函数算出,然后
12,这样做的优点比起计算各分条长然后累加的办法来说有二:①避免了误差
累积,提高了精确度;②使计算过程由繁琐变为简单。
五、最小安全系数Ks的得出原理:第二次计算得出的Ks和第一次的进行比较,若小于第一次的值则留下作为Kmin的临时值赋给变量Kt,否则用第一次的值作为临时值,第三次计算得出的再和这个临时值亦即变量Kt作比较,若小于则赋给变量Kt,依此类推,最后即可得出最小安全系数Kmin。
六、渗透力计算原理:由于渗透力的计算比较麻烦(需要确定的边界条件较多),故考虑作简单化处理,即只考虑两个因素:①浸润面积占滑体总面积的百分比;②作用方向:一般按平行于坡面方向考虑;③力的大小:由于渗透力是体积力,按单宽计
3
算就是面积力,即Jj*Sw*i*S。
七、最小安全系数所对应圆心及滑弧的显示:有两种方法可以实现这一目的:①每次计算均用变量来存储滑动圆心、滑弧半径及其两端点的坐标,此法比较麻烦且效率较低,影响计算速度;②总体计算完后,再算一次。由于此次只需计算至Ks=Kmin的那一次即可结束,并且此次计算过程不需要图形显示(只需计算完成后
显示最终结果),故计算速度非常快,一般不超过1至10秒。
八、延时显示功能:考虑到多媒体教学时的需要,编程时用了五种不同的速度来执行计算过程,分别对应键盘上的F1至F5键,并且还可以自定义更慢的速度以方便教
学、讲解或观察,如左图和右图所示。
九、分条宽度b=R/n:在人工试算时,一般取n=10,但是计算机却可以取n=10、20、40、
4
80、100、甚至n=1000来计算(自定义),例如上面左图是n=100的情形。
十、坐标原点的选取:为了能够比较居中地显示图形,选取y轴通过滑面AB中点,x轴通过滑坡底面下距离滑坡底面d=H*(1/3)取整数处。这样,当H=20m时,滑面底部B点纵坐标为7m,而滑面顶、底部A、B两点的横坐标始终为负正相反数,如下图。
十一、交点坐标的得出:通过直线、圆的方程进行代数求解得出,当有两个解时,由边界条件自动判断具体取哪一个值。
十二、安全系数计算公式:
①滑体内无渗流通过时,计算公式为:
Ks(Wcostan)cLbtan(hcos)cL
b(hsin)(Wsin)iiiiiiiii②滑体内有渗流通过时,计算公式为:
Ks
btan(hi1cosi)'btan(hi2cosi)cL
b(hi1sini)satb(hi2sini)以上是编程思路,下面就历界教材上的部分例习题作计算讲解及说明。
【例1】中专教材《土力学》(成都水力发电学校 李道荣)第168页例题:某均质粘性土坡,坡高为5m,坡比为1:2,土的有关指标为:17.66KN/m3,c=9.81kPa,
5
15,试用圆弧条分法求土坡稳定的最小安全系数Kmin。
解:由于计算量太大,为避免出错,教材中采用列表试算,得出其中一次计算的安全系数为Ks=1.55;现采用自编程序计算:得Kmin=1.524091。如图:
【例2】大专教材《土力学》(黄河水利出版社 务新超)第160页例题:一均质粘性土坡,高为20m,坡比为1:3,土的有关指标为:重度17.66KN/m3,粘聚力c=9.81kPa,内摩擦角20,试用太沙基公式(即圆弧条分法)求土坡稳定的最小安
6
全系数Kmin。
解:由于计算量太大,为避免出错,教材中采用列表的方法进行试算,得出其中一次计算的安全系数为Ks=1.83;现采用自编程序计算:得Kmin=1.4878(图略)。
本例题也被收录在高等学校教材《土力学》(中国水利水电出版社 天津大学 杨进良)第195页中。
除边坡稳定分析程序外,本人还编制了赤平投影分析、抗剪强度指标计算、附加应力计算等工程地质与土力学方面的可视化实用小程序,为了方便各位专家和同仁共同研究或探讨,特附上这几个小程序,也可致电邮至:ZhangHqing@sina.com索取。
下面列出圆弧条分法边坡稳定分析程序的主要程序代码,以供参考:
Option Explicit
Private Declaration As String Private Sub bgvuvuti_Click() Form2.Show
Form2.FontBold = False Form2.FontSize = 20
Form2.FontName = \"隶书\"
Form2.ForeColor = RGB(50, 100, 200)
Form2.Print \"本软件适用于圆弧条分法边坡稳定分析\" Form2.Print \" 如需示范演算,请直接单击<开始>按纽\"
Form2.Print \" 如需了解软件功能,请单击最大化按钮,\" Form2.Print \" 软件使用说明书便会呈现于窗体右侧\" Form2.FontSize = 15
Form2.Print \"Copyrihgt : 张华庆 December 9th . 2002\" End Sub
Private Sub ffliyi_Click()
H = 5: a = 26.565: v = 17.66: fia = 15: c = 9.81 End Sub
Private Sub fflioq_Click()
H = 25: a = 24.33: v = 17.66: fia = 17.5: c = 9.81 End Sub
Private Sub fflisf_Click()
H = 120: a = 30: v = 18: fia = 22.5: c = 10 End Sub
Public Sub Form_Load() Show
7
H = 20: a = 30
v = 18.78: fia = 22.58: c = 9.8
fi = fia / 57.2958 '将φ值转化为弧度
t10 = 10 '给分条宽度b=R/t10中的t10赋一初值 jt = 0 '给浸润面积系数jt赋一初值零(不考虑其影响) slptime = 0 '说明延时默认值为零 Form1.DangerArc.Enabled = False Form1.DisPlay.Enabled = False Form1.HtA.Enabled = False Label1.FontSize = 15
Declaration = \" 软件使用说明书 \" & Chr(10) & \"著作人声明:本软件一切版权属 zhq 所有,未经版权所有人允许,任何单位和个人不得以任何形式复制和传播,否则将依法追究其经济赔偿责任。\" & Chr(10) & \"功能键及菜单说明:\" & Chr(10) & \"延时功能:F1~F5(快→慢);\" & \"分条因数:F8、F9、F11、F12(少→多,即分条宽由宽到窄)\" & Chr(10) & \"渗透力:由于影响因素较多,仅考虑浸润面积为滑体面积的1/6、1/4、1/2、3/4等几种情形,作近似计算,且假定A1B1斜率即为渗流水力坡降。\" & Chr(10) & \"滑体面积:自动计算并显示滑动面以上土体的面积。\" & Chr(13) & \"F5:运行(Start)\" & Space(2) & \"F10:暂停/继续(Break/Continue)\" & Chr(10) & \"Ctrl+F6:设置速度(毫秒)\"
Label1.Caption = Declaration End Sub
Private Sub A0_Click() jt = 0 End Sub
Private Sub A2_Click() jt = 0.2 End Sub
Private Sub A4_Click() jt = 0.4 End Sub
Private Sub A6_Click() jt = 0.6 End Sub
Private Sub A8_Click() jt = 0.8 End Sub
Private Sub HtA_Click() Dim S$
S = MsgBox(\"滑体面积为\" & Str$(Sa) & \"(平方米)\渗流影响分析\") End Sub
Private Sub Picture1_MouseDown(Button As Integer, Shift As Integer, x As Single, Y As Single) Unload Form2 '此语句相当于If Button = 1 or 2 Then Unload Form2 If Button = 2 Then Form1.PopupMenu LPRINT '弹出式菜单应用例子 End Sub
Private Sub Fprint_Click()
8
slptime = 0
End Sub '打印窗体暂时的不要,改为延时功能 Private Sub Cprint_Click() slptime = 200
End Sub '打印窗体暂时的不要,改为延时功能 Private Sub ft10_Click(Index As Integer) t10 = 10 End Sub
Private Sub ft20_Click(Index As Integer) t10 = 20 End Sub
Private Sub ft40_Click(Index As Integer) t10 = 40 End Sub
Private Sub ft80_Click(Index As Integer) t10 = 80 End Sub
Private Sub ft100_Click(Index As Integer) t10 = 100 End Sub
Private Sub menu1_Click(Index As Integer) On Error Resume Next
H = InputBox(\"请输入土坡高度H(m):\输入对话框\End Sub
Private Sub Pobi_Click() Dim pbn As Single On Error Resume Next
pbn = InputBox(\"请输入坡比 1:n 中 n 的值:\输入对话框\ a = 57.2957795 * Atn(1 / pbn) End Sub
Private Sub Pojk_Click()
On Error Resume Next
a = InputBox(\"请输入坡度角α(°):\输入对话框\End Sub
Private Sub menu3_Click(Index As Integer) On Error Resume Next
v = InputBox(\"请输入土的容重γ:\输入对话框\End Sub
Private Sub menu4_Click(Index As Integer) On Error Resume Next
fia = InputBox(\"请输入滑面或土的摩擦角φ:\输入对话框\End Sub
Private Sub menu5_Click(Index As Integer) On Error Resume Next
9
c = InputBox(\"请输入滑面或土的粘聚力c:\输入对话框\ End Sub
Private Sub SpeedSet_Click() On Error Resume Next
slptime = InputBox(\"请输入延时毫秒数:\设置速度对话框\End Sub
Sub Start_Click()
Call Clsn_Click
Form1.Label1.Caption = \"程序正在运行中,请稍候„„\" Form1.DangerArc.Enabled = False Form1.DisPlay.Enabled = False
ra = a / 57.2958 '将a值转化为弧度 Ax = -H / Tan(ra) / 2: Bx = H / Tan(ra) / 2
By = CInt(H / 3): Ay = By + H '将坡脚纵坐标四舍五入取整 fi = fia / 57.2958 '将φ值转化为弧度 L = Sqr((Ax - Bx) ^ 2 + (Ay - By) ^ 2) kab = (Ay - By) / (Ax - Bx) Erase Ks: Erase Kn
Kmin = 1000 '为Kmin赋一初值1000,以便比较 Kjin = 1000 '为Kjin赋一初值1000,以便比较
Call SeaCheR '调用过程SeaCheR ,留意参数的传送问题 jn = j
Form1.DangerArc.Enabled = True Form1.HtA.Enabled = True
Form1.Label1.Caption = Declaration Beep
DangerArc.SetFocus End Sub
Sub SeaCheR() '*****此过程用来即时显示每一试算圆心时的滑弧及分条图形***** '笛卡尔坐标系:过图片框底部向右(→)为X轴,过坡面中点向上(↑)为Y轴 j = 0
For B1x = Bx - Bx / 4 To Bx + Bx / 4 Step Bx / 8 If B1x < Bx Then
B1y = g(B1x) ElseIf Bx <= B1x Then B1y = By End If
For Ox = Ax To Bx + H / 2 Step L / 10
For Oy = Ay + H / 4 To Ay + 4.5 * H Step L / 10
j = j + 1 '此句从zhq中提前后便修正了Ksi显示中的错误
Call A1Search '进入用试算法求得A1点坐标过程
10
Form1.Cls
Form1.Picture1.Cls
Picture1.Scale (-3 * H, 6 * H)-(3 * H, 0) '第一次定义图片框尺寸比例 Call zhq '进入计算安全系数并绘制分条过程
Call FPsPlay '进入在窗体和图片框中显示文字及数据过程 Ua = Atn(Abs((Oy - A1y) / (Ox - A1x)))
If Ox < B1x Then '若O点在B1点左侧,则Ub为锐角 Ub = Atn(Abs((Oy - B1y) / (Ox - B1x))) ElseIf Ox = B1x Then
Ub = 3.14159 / 2 '若O点在B1点之上,则Ub =90°(直角) ElseIf Ox > B1x Then
Ub = 3.14159 - Atn(Abs((Oy - B1y) / (Ox - B1x))) '若O点在B1点右侧,则Ub为钝角
End If
Picture1.Line (Ox - H / 12, Oy)-(Ox + H / 12, Oy), RGB(120, 100, 40) '画过圆心的水平线 Picture1.Line (Ox, Oy - H / 12)-(Ox, Oy + H / 12), RGB(120, 100, 40) '画过圆心的铅直线
Picture1.Circle (Ox, Oy), R, RGB(150, 0, 150), -(3.1416 + Ua), -(6.2832 - Ub) '画滑动面圆弧(含径线)
Picture1.DrawStyle = 3
Picture1.Line (Ox, Oy - H / 8)-(Ox, f(Ox)), RGB(0, 200, 40) '画过滑弧圆心的铅垂点划线 Picture1.DrawStyle = 0
DoEvents '当程序运行时允许转让控制权
Sleep (slptime) '增加延时功能^下面三条语句的作用是把已画过的圆心和滑弧及时抹掉
Next Oy Next Ox Next B1x End Sub
Sub DrawPicture()
Picture1.Scale (-3 * H, 6 * H)-(3 * H, 0)
'*************下面语句圈定试算圆心范围**************** Dim Ix As Single Dim Iy As Single
For Ix = Ax To Bx + H / 2 Step L / 10
For Iy = Ay + H / 4 To Ay + 4.5 * H Step L / 10 Picture1.PSet (Ix, Iy) Next Iy Next Ix
'---------------------------------------------------------------------- Picture1.Line (Ax, Ay)-(Bx, By) '画坡面线
Picture1.Line (Ax - L, Ay)-(Ax, Ay) '画坡肩水平线(自左向右) Picture1.Print \"A\"
Picture1.Line (Bx + L, By)-(Bx, By) '画坡底水平线(自右向左)
11
Picture1.Print \"B\"
Picture1.CurrentX = -2.9 * H: Picture1.CurrentY = 5.9 * H
Picture1.Print \"容重γ=\"; v; \"kn/m^3\"; Space(20); \"摩擦角 φ=\"; CInt(fia * 100) / 100; \"°\"; Space(20); \"粘聚力c=\"; c; \"kPa\"
Picture1.CurrentX = -2.9 * H
Picture1.Print \"Ax=\"; CSng(Ax); Space(5); \"Ay=\"; CSng(Ay); Space(13); \"Bx=\"; CSng(Bx); Space(5); \"By=\"; CSng(By); Space(18); \"坡比为 1:\"; CInt(100 * (Bx - Ax) / H) / 100
Picture1.Line (Ox - H / 12, Oy)-(Ox + H / 12, Oy), RGB(0, 0, 0) '画过圆心的水平线 Picture1.Print \"O\"
Picture1.Line (Ox, Oy - H / 12)-(Ox, Oy + H / 12), RGB(0, 0, 0) '画过圆心的铅直线 Ua = Atn(Abs((Oy - A1y) / (Ox - A1x)))
If Ox < B1x Then '若O点在B1点左侧,则Ub为锐角 Ub = Atn(Abs((Oy - B1y) / (Ox - B1x))) ElseIf Ox = B1x Then
Ub = 3.14159 / 2 '若O点在B1点之上,则Ub =90°(直角) ElseIf Ox > B1x Then
Ub = 3.14159 - Atn(Abs((Oy - B1y) / (Ox - B1x))) '若O点在B1点右侧,则Ub为钝角
End If
Picture1.Circle (Ox, Oy), R, RGB(250, 50, 140), -(3.1416 + Ua), -(6.2832 - Ub) '画滑动面圆弧(含径向线,若不要径向线,请将\"-\"号去掉)
Picture1.CurrentX = A1x: Picture1.CurrentY = A1y '注意:A1x、A1y若尚未赋值
Picture1.Print \"A1\" '则A1、B1是坐标错误的
Picture1.CurrentX = B1x: Picture1.CurrentY = B1y Picture1.Print \"B1\"
For y1 = 0 To Oy - H / 10 Step Oy / 20
Picture1.Line (Ox, y1)-(Ox, y1 + R / 25) '画过滑弧圆心的铅垂(虚)线 Next y1 End Sub
Sub A1Search() ' #####(本过程主要以计算A1点坐标为目的)##### OA = Sqr((Ox - Ax) ^ 2 + (Oy - Ay) ^ 2)
R = Sqr((Ox - B1x) ^ 2 + (Oy - B1y) ^ 2) '(取滑弧半径R=OB1)
If OA = R Then
A1x = Ax: A1y = Ay ElseIf OA < R Then
A1y = Ay: A2y = Ay: Avy = Ay '此处曾经漏掉(Av y= Ay)酿成大错
For A1x = Ax - 4 * L To Ax Step 0.4 '此处曾经为Ax - 2* L,结果出现有时找不到交点的毛病,
A1x = A1x: A2x = A1x + 0.4 '致使A1x以Ax输出,Ua计算错误,以致图形错误显示
k: Avx = (A1x + A2x) / 2 '二分法计算A1点坐标开始(坡顶直线段) OA1 = Sqr((Ox - A1x) ^ 2 + (Oy - A1y) ^ 2)
12
OA2 = Sqr((Ox - A2x) ^ 2 + (Oy - A2y) ^ 2) OAv = Sqr((Ox - Avx) ^ 2 + (Oy - Avy) ^ 2) If OA1 = R Then
A1x = A1x: A1y = A1y Exit For
ElseIf OA1 > R And OA2 < R Then If OAv = R Then
A1x = Avx: A1y = A1y Exit For
ElseIf OAv < R Then
A1x = A1x: A2x = Avx: GoTo k ElseIf OAv > R Then
A1x = Avx: A2x = A2x: GoTo k End If End If Next A1x
ElseIf OA > R Then
For A1x = Ax To Bx Step 0.4
A1x = A1x: A2x = A1x + 0.4
q: Avx = (A1x + A2x) / 2 '二分法计算A1点坐标开始(坡面斜线段) A1y = (By - Ay) * (A1x - Ax) / (Bx - Ax) + Ay A2y = (By - Ay) * (A2x - Ax) / (Bx - Ax) + Ay Avy = (By - Ay) * (Avx - Ax) / (Bx - Ax) + Ay OA1 = Sqr((Ox - A1x) ^ 2 + (Oy - A1y) ^ 2) OA2 = Sqr((Ox - A2x) ^ 2 + (Oy - A2y) ^ 2) OAv = Sqr((Ox - Avx) ^ 2 + (Oy - Avy) ^ 2) If OA1 = R Then
A1x = A1x: A1y = A1y Exit For
ElseIf OA1 > R And OA2 < R Then If OAv = R Then
A1x = Avx: A1y = Avy Exit For
ElseIf OAv < R Then
A1x = A1x: A2x = Avx: GoTo q ElseIf OAv > R Then
A1x = Avx: A2x = A2x: GoTo q End If End If Next A1x End If End Sub
Sub FPsPlay() '此过程功能为分别在窗体和图片框中显示文字
13
'以下语句在窗体中显示文字
Form1.ForeColor = RGB(0, 0, 250)
Form1.Print \"坡高为H = \"; Format(H, \"0.00\"); \"米\"; Tab(22); \"坡度角为α=\"; CInt(a * 100) / 100; \"°\"; Tab(46); \"A1点坐标为 (\"; Format(A1x, \"00.0000\"); \",\"; Format(A1y, \"00.0000\"); \")\"
Form1.Print \"坡面长为L=\"; Format(L, \"0.00\"); \"米\"; Tab(46); \"B1点坐标为 (\"; Format(B1x, \"0.0000\"); \",\"; Format(B1y, \"0.0000\"); \")\"
Form1.Print \"试算圆心坐标为Oi (\"; Tab(18); Format(Ox, \"0.0000\"); \",\"; Format(Oy, \"0.0000\"); \")\"; Tab(48); \"计算次数 i=\"; j; \" 次\"
Print \"试算滑弧的半径为 R=\"; CSng(R); Tab(30); \"米\\"试算滑弧长为 ΣL⌒= \"; Format$(SumLi, \"0.0000\"); Tab(70); \"米\"
Form1.Print
Form1.Print \"OA=\"; Format(OA, \"0.0000\") Form1.Print \"OA1=\"; Format(OA1, \"0.0000\") Print
Form1.Print \"Area=\"; Format(Sa, \"0.00\") Form1.Print \"分条数fts=\"; fts Print
If jt < 0.01 Then
Form1.Print \"K(\"; j; \")=\"; Format(Ks(j), \"0.00000\") Form1.ForeColor = RGB(250, 0, 0)
Form1.Print \"Kmin=\"; Format(Kmin, \"0.000000\") ElseIf jt >= 0.01 Then
Form1.Print \"考虑渗流影响时:\"
Form1.Print \"K(\"; j; \")=\"; Format(Kn(j), \"0.00000\") Form1.ForeColor = RGB(250, 0, 180) Print \"Kjin=\"; Format(Kjin, \"0.000000\") End If
'以下语句在图片框中显示文字 Form1.ForeColor = RGB(0, 0, 250) Picture1.Scale (-3 * H, 6 * H)-(3 * H, 0)
'*************下面语句圈定试算圆心范围**************** Dim Ix As Single Dim Iy As Single
For Ix = Ax To Bx + H / 2 Step L / 10
For Iy = Ay + H / 4 To Ay + 4.5 * H Step L / 10 Picture1.PSet (Ix, Iy) Next Iy Next Ix
'---------------------------------------------------------------------- Picture1.Line (Ax, Ay)-(Bx, By) '画坡面线
Picture1.Line (Ax - L, Ay)-(Ax, Ay) '画坡肩水平线 Picture1.Print \"A\"
Picture1.Line (Bx + L, By)-(Bx, By) '画坡底水平线 Picture1.Print \"B\"
14
Picture1.CurrentX = -2.9 * H: Picture1.CurrentY = 5.9 * H
Picture1.Print \"容重γ=\"; v; \"kn/m^3\"; Space(20); \"摩擦角 φ=\"; CInt(fia * 100) / 100; \"°\"; Space(20); \"粘聚力c=\"; c; \"kPa\"
Picture1.CurrentX = -2.9 * H
Picture1.Print \"Ax=\"; CSng(Ax); Space(5); \"Ay=\"; CSng(Ay); Space(13); \"Bx=\"; CSng(Bx); Space(5); \"By=\"; CSng(By); Space(18); \"坡比为 1:\"; CInt(100 * (Bx - Ax) / H) / 100
Picture1.CurrentX = A1x: Picture1.CurrentY = A1y Picture1.Print \"A1\"
Picture1.CurrentX = B1x: Picture1.CurrentY = B1y Picture1.Print \"B1\" End Sub
Sub CountDArc()
Static i As Integer
i = 0 '计算最危险滑弧圆心及半径开始 For B1x = Bx - Bx / 4 To Bx + Bx / 4 Step Bx / 8 If B1x < Bx Then
B1y = g(B1x) ElseIf Bx <= B1x Then B1y = By End If
For Ox = Ax To Bx + H / 2 Step L / 10
For Oy = Ay + H / 4 To Ay + 4.5 * H Step L / 10 i = i + 1 DoEvents
If jt < 0.01 Then '分别考虑并选择有无渗透力影响时的
If i = ju Then Exit Sub '土坡稳定最小安全系数
ElseIf jt >= 0.01 Then If i = jv Then Exit Sub End If Next Oy Next Ox Next B1x End Sub
Private Sub DangerArc_Click() Call CountDArc
Call A1Search '寻找A1点坐标开始--结束 Form1.Cls
Form1.Picture1.Cls Call zhq If jt < 0.01 Then
Print \"总运算次数为jn=\"; jn; \"次\算到其中第i= \"; ju; \"次时获得最危险滑弧\" ElseIf jt >= 0.01 Then
15
Print \"总运算次数为jn=\"; jn; \"次\算到其中第i= \"; jv; \"次时获得最危险滑弧\" End If
Print \"不考虑渗流影响时\";
Form1.ForeColor = RGB(250, 0, 0) Print \"Kmin=\"; \"Ks(\"; ju; \")=\"; Ks(ju), Form1.ForeColor = RGB(0, 0, 250) Print Tab(58); \"坡高为H =\"; H; \"米\"
Print \"最危险滑弧的圆心坐标为 Oi(\"; CSng(Ox); \坡度角α=\"; CInt(a * 100) / 100; \"°\"
Print \"最危险滑弧的半径为 R=\"; CSng(R); Tab(32); \"米\\"最危险滑弧长为 ΣL⌒= \"; CSng(SumLi); Tab(74); \"米\"
Print \"分条宽ftk=R/\"; t10 Print \"分条数fts =\"; fts; \"条\" Print
Form1.Print \"Area=\"; Sa Print
If jt < 0.01 Then
Form1.Print \"K(\"; ju; \")=\"; Ks(ju) Form1.Print \"最小安全系数为:\" Form1.ForeColor = RGB(250, 0, 0) Print \"Kmin=\"; Kmin ElseIf jt >= 0.01 Then
Form1.Print \"考虑渗流影响时:\" Form1.Print \"K(\"; jv; \")=\"; Kn(jv) Form1.ForeColor = RGB(250, 0, 180) Print \"Kjin=\"; Kjin End If
Form1.ForeColor = RGB(0, 0, 250) Call DrawPicture
Form1.DisPlay.Enabled = True Form1.DisPlay.SetFocus End Sub
Sub Clsn_Click() Form1.Cls Picture1.Cls
Picture1.Scale (-3 * H, 6 * H)-(3 * H, 0) End Sub
Sub DisPlay_Click() Form1.Cls
Form1.ForeColor = RGB(0, 0, 250)
Form1.Print \"最危险圆心坐标为Oi (\"; Tab(20); CSng(Ox); \点坐标为 (\"; CSng(A1x); \
Form1.Print \"坡高为H = \"; H; \"米\"; Tab(22); \"坡度角为α=\"; CInt(5729.58 *
16
Atn(Abs(kab))) / 100; \"°\"; Tab(46); \"B1点坐标为 (\"; CSng(B1x); \
If jt < 0.01 Then
Form1.Print \"坡面长为L=\"; L; \"米\"; Tab(44); \"计算次数 i=\"; ju; \" 次\" ElseIf jt >= 0.01 Then
Form1.Print \"坡面长为L=\"; L; \"米\"; Tab(44); \"计算次数 i=\"; jv; \" 次\" End If
Print \"最危险滑弧的半径为 R=\"; CSng(R); Tab(32); \"米\\"最危险滑弧长为 ΣL⌒= \"; CSng(SumLi); Tab(74); \"米\"
Form1.Print
Form1.Print \"OA=\"; OA Form1.Print \"OA1=\"; OA1 Print
Form1.Print \"Area=\"; Sa
Form1.Print \"分条数fts=\"; fts Print
If jt < 0.01 Then
Form1.Print \"K(\"; ju; \")=\"; Ks(ju) Form1.ForeColor = RGB(250, 0, 0) Form1.Print \"Kmin=\"; Kmin ElseIf jt >= 0.01 Then
Form1.Print \"K(\"; jv; \")=\"; Kn(jv) Form1.ForeColor = RGB(250, 0, 180) Print \"Kjin=\"; Kjin End If
Form1.ForeColor = RGB(0, 0, 250) DangerArc.SetFocus End Sub
Private Sub END_Click() End End Sub
Private Sub time_0_Click() slptime = 0 End Sub
Private Sub time_20_Click() slptime = 20 End Sub
Private Sub time_50_Click() slptime = 50 End Sub
Private Sub time_100_Click() slptime = 100 End Sub
Private Sub time_500_Click() slptime = 500
17
End Sub
Private Sub ZDYWIDTH_Click() On Error Resume Next
t10 = InputBox(\"请输入滑体分条宽度\" & Chr$(10) & \" b = R / t 中 t 的值:\分条宽度对话框\
End Sub
Private Sub Suckset_Click() On Error Resume Next
jt = InputBox(\"请输入浸润面积/滑体面积之比例系数\" & Chr$(10) & \"如为25%则输入0.25等等\浸润面积对话框\
End Sub
Option Explicit
Public Ax#, Ay#, Bx#, By#, Ox#, Oy#, kab#, fi#, SumLi#
Public A1x#, A1y#, A2x#, A2y#, Avx#, Avy#, B1x#, B1y#, D# Public R1!, R2!, H!, L!, R!, y1!, OA!, OA1!, OA2!, OAv!, Jrm!, jt! Public j As Integer, jn As Integer, fts As Long, Kmin As Single Public Ks(10000) As Single, slptime As Integer Public a!, ra!, v!, c!, fia As Single, t10 As Integer Public Ua As Double, Ub As Double, Sa As Single Public Kn(10000) As Single, Kjin As Single Public ju As Integer, jv As Integer
Public Declare Sub Sleep Lib \"kernel32\" (ByVal dwMilliseconds As Long) Function f(ByVal x As Double) As Double
f = Oy - Sqr(R ^ 2 - (x - Ox) ^ 2) '滑弧下半圆方程 End Function
Function g(ByVal x As Double) As Double
g = kab * (x - Ax) + Ay '坡面线AB的方程 End Function
Sub zhq() '********本过程计算土坡稳定安全系数Ks********* Dim Xi#, Hi#, ai#, Sumhs#, Sumhc#, A1B1#, b#, q#, u#, XL#, XR# Dim i As Long, SumA As Single Dim Ii As Single, m&, n&
b = R / t10 '分条宽度为滑弧半径的1/t10
m = CLng((Ox - A1x - b / 2) / b): n = CLng((Ox - B1x + b / 2) / b) '-------------------------------------------------------------------------------------- For i = n To m '分条计算各土条高并绘图部分开始
Xi = Ox - b * i '返回第i条土条中心的横坐标Xi XL = Xi - b / 2: XR = Xi + b / 2 'XL,XR分别为第i条土条左、右侧的横坐标 If A1x < Xi And Xi <= Ax Then Hi = Abs(Ay - f(Xi))
18
If A1x < XL Then
Form1.Picture1.Line (XL, f(XL))-(XL, A1y), RGB(20, 20, 240) '前面必须加Form1.(∵是在Module里划直线)
End If
If XR <= Ax Then
Form1.Picture1.Line (XR, f(XR))-(XR, A1y), RGB(20, 20, 240) '前面必须加Form1.(∵是在Module里划直线)
ElseIf Ax < Xi + b / 2 Then
Form1.Picture1.Line (XR, f(XR))-(XR, g(XR)), RGB(20, 20, 240) End If
ElseIf Ax < Xi And Xi < Bx Then Hi = Abs(g(Xi) - f(Xi))
If XR < Bx And XR < B1x Then
Form1.Picture1.Line (XR, f(XR))-(XR, g(XR)), RGB(20, 20, 240) ElseIf Bx <= XR And XR < B1x Then
Form1.Picture1.Line (XR, f(XR))-(XR, By), RGB(20, 20, 240) End If
ElseIf Bx <= Xi And Xi < B1x Then Hi = Abs(By - f(Xi))
If XR < B1x Then
Form1.Picture1.Line (XR, f(XR))-(XR, By), RGB(20, 20, 240) End If End If
ai = Atn(i / Sqr(t10 ^ 2 - i ^ 2)) '∵希腊字母不能作VB变量名∴以ai代表αi SumA = SumA + Hi * b '求拟滑体面积
Sumhc = Sumhc + Hi * Cos(ai): Sumhs = Sumhs + Hi * Sin(ai) '以上计算∑hi*cosαi,∑hi*sinαi两项 Next i
'----------------------------------------------------------------------------------------- fts = m - n + 1 'fts为分条总数(如-1~3条共3-(-1)+1=5条) A1B1 = Sqr((A1x - B1x) ^ 2 + (A1y - B1y) ^ 2)
q = A1B1 / R / 2: u = q / Sqr(1 - q ^ 2) ' tg(α/2)=u,α为滑弧圆心角 ' ∵A1B1 / 2 < R ∴q < 1 ∴1 - q ^ 2 > 0 ∴u > 0 ∴Atn(u)=锐角 SumLi = 2 * R * Atn(u) 'α=2Atn(u) q=sin(α/2) If v * b * Sumhs > 0 Then
Ks(j) = (v * b * Tan(fi) * Sumhc + c * SumLi) / (v * b * Sumhs) '用公式计算Ksi的值
If Ks(j) < Kmin Then '每次计算后将较小的Ks值赋给Kmin Kmin = Ks(j)
ju = j '将最小安全系数所对应的运算次数记录下来 End If
End If '本IF语句的作用是防止除数分母为零
Sa = SumA '把滑体面积传送给全局变量Sa,以便随时可以输出 If v * b * Sumhs > 0 And A1x <> B1x And A1B1 < 2 * R Then
Ii = Abs((A1y - B1y) / (A1x - B1x)) '取水力坡度为A1B1的斜率
19
D = Sqr(R ^ 2 - (A1B1 / 2) ^ 2) '近似取O点到A1B1的距离为力臂 Jrm = SumA * jt '假定浸润面积为滑体面积的jt倍
Kn(j) = (v * b * Tan(fi) * Sumhc + c * SumLi) / (v * b * Sumhs + 9.81 * Ii * Jrm * D / R)
If Kn(j) < Kjin Then Kjin = Kn(j)
jv = j '将最小安全系数所对应的运算次数记录下来 End If End If End Sub
20
因篇幅问题不能全部显示,请点此查看更多更全内容