圆弧条分法土坡稳定分析的VB电算解法5.18

1、圆弧条分法土坡稳定分析的vb电算解法解题思路说明与讲解云南省水利水电学校张华庆圆弧条分法土坡稳定分析的vb电算解法云南省水利水电学校张华庆【摘要】目前，大多数版本的土力学教材在讲到均质粘性土边坡稳定性分析方法时，均采用圆弧条分法，该法存在大量的作图及试算，若用人工求解，则不但会耗费大量宝贵的时间，而且还难以得到精确解。现就介绍作者自编的粘性土边坡简单圆弧条分法的vb电算解法。【关键词】土坡稳定圆弧条分法vb电算全自动可视化引言：迄今为止，求解均质粘性土边坡稳定最小安全系数ksmin，通常都采用瑞典圆弧法或圆弧条分法，由于该法存在大量的手工作图以及滑弧试算，因此如果采用过去传统的人工方法求解，则

2、因为效率低下，不但浪费稿纸和笔墨、浪费大量宝贵的时间，而且还难以得出较为精确的解答。为此，笔者利用visual basic 6.0编写和开发了可视化的土坡稳定分析程序，该程序及解法具有如下特点：作图与试算全过程可视化跟踪显示，其优越性是过去的dos计算程序所无法比拟的；速度快，能在一分钟之内（一般仅需10秒左右）找到最危险滑弧以及土坡稳定最小安全系数ksmin；充分体现圆弧条分法的精髓，精确度高；简单易行，并且经过若干教材的例习题演算，证明结果一致；演算过程进行速度可以调节，便于观察和教学；程序大小仅72kb，但却功能完善，完全可以满足教学需要（当然也完全可以用于工程计算）；自动化程度高：仅需

3、输入必要的坡高、坡比（或坡度角），土的容重以及土的抗剪强度参数、c的值，即可自动完成计算，因此为全自动化设计，功能上明显优越于目前有的半自动化边坡稳定分析程序（如有的为dos程序；有的虽然是可视化程序，但却是半自动化的，需要在autocad里人为画出试算滑弧，然后才能让机器去计算土坡稳定最小安全系数，这自然影响时效及精度，因为你不可能有时间和耐心去画多达100条滑弧！就别说更多的）；分条数可突破人为影响及限制，例如可以达到1000条之多；无论何种条件下，当分条数逐渐增多时，土坡稳定最小安全系数ksmin均能趋于一个稳定的极限值，证明计算方法以及计算结果是正确的且是可靠的和稳定的；可以考虑渗透力

4、的影响。下面就本程序的编程思路作一些必要的说明和讲解，以便和有兴趣致力于此方面研究的同仁共同探讨，并希望得到专家的指点，以便使程序变得更加完美和人性化。一、滑弧圆心的确定变得简单化：试算法首先要确定的是滑弧的圆心，由于过去是采用人工作图计算，工作量或者说劳动强度可想而知，于是必须尽可能地缩小滑动圆心的范围，这样一来才需要定出m点坐标（4.5h，2h）、确定o点所需的参数125，235，以及利用十字交叉方法确定滑动圆心oi点的做法就是必需的和可以理解的；但是对于计算机来说，它最不怕的就是计算，因而这种做法就显得多余或者说不必要了，我们可以让它在平面上可能的区域（即所画圆弧至少有一端能交于坡面、坡

5、顶或坡底，而另一端也不至于悬空的点）在纵横两个方向每间隔一定的点计算一次去寻找最危险滑弧，问题因此有了更为简单而圆满的解决办法。二、分条面积计算简单化：由于分条宽度b是固定的，所以分条面积取决于分条高度hi，笔者采用解析几何方法加以处理：列出坡面及坡顶、坡底的方程，并且列出圆弧的方程，再列出第i个土条的竖线方程；用公式法或试算法（迭代原理）求出交点坐标，再用两点间距离公式计算出hi；分条面积计算公式：sihi*b(“*”号即“”号)；由于分条数可以很多，当分条数多到一定程度（一般当取b=r/100或更大）时滑体两端三角形面积部分对总面积的正负或加减影响可以不予考虑。三、滑弧与坡面的交点分析：程

6、序分别考虑了滑弧与土坡上游交点的两种情况，即交于坡顶和交于上半坡面；滑弧与土坡下游交点的两种情况，即交于坡底和交于下半坡面。这几种情况都是圆弧条分法所要求考虑的，均可由解析几何方法得出其交点坐标，然后由上下两交点之间的水平间距d及b=r/n计算出分条数fts。四、滑弧半径及滑弧长的确定：在滑弧圆心确定后，滑弧半径通过两点间距离公式确定，即；至于滑弧长，可通过公式计算滑弧总长，即在分条数fts确定后，水平方向正负两边圆心角即可用反三角函数算出，然后，这样做的优点比起计算各分条长然后累加的办法来说有二：避免了误差累积，提高了精确度；使计算过程由繁琐变为简单。五、最小安全系数ks的得出原理：第二次计

7、算得出的ks和第一次的进行比较，若小于第一次的值则留下作为kmin的临时值赋给变量kt，否则用第一次的值作为临时值，第三次计算得出的再和这个临时值亦即变量kt作比较，若小于则赋给变量kt，依此类推，最后即可得出最小安全系数kmin。六、渗透力计算原理：由于渗透力的计算比较麻烦（需要确定的边界条件较多），故考虑作简单化处理，即只考虑两个因素：浸润面积占滑体总面积的百分比；作用方向：一般按平行于坡面方向考虑；力的大小：由于渗透力是体积力，按单宽计算就是面积力，即。七、最小安全系数所对应圆心及滑弧的显示：有两种方法可以实现这一目的：每次计算均用变量来存储滑动圆心、滑弧半径及其两端点的坐标，此法比较麻

8、烦且效率较低，影响计算速度；总体计算完后，再算一次。由于此次只需计算至kskmin的那一次即可结束，并且此次计算过程不需要图形显示（只需计算完成后显示最终结果），故计算速度非常快，一般不超过1至10秒。八、延时显示功能：考虑到多媒体教学时的需要，编程时用了五种不同的速度来执行计算过程，分别对应键盘上的f1至f5键，并且还可以自定义更慢的速度以方便教学、讲解或观察，如左图和右图所示。九、分条宽度b=r/n：在人工试算时，一般取n=10，但是计算机却可以取n=10、20、40、80、100、甚至n=1000来计算（自定义），例如上面左图是n=100的情形。十、坐标原点的选取：为了能够比较居中地显示

9、图形，选取y轴通过滑面ab中点，x轴通过滑坡底面下距离滑坡底面d=h*（1/3）取整数处。这样，当h20m时，滑面底部b点纵坐标为7m，而滑面顶、底部a、b两点的横坐标始终为负正相反数，如下图。十一、交点坐标的得出：通过直线、圆的方程进行代数求解得出，当有两个解时，由边界条件自动判断具体取哪一个值。十二、安全系数计算公式：滑体内无渗流通过时，计算公式为：滑体内有渗流通过时，计算公式为：以上是编程思路，下面就历界教材上的部分例习题作计算讲解及说明。【例1】中专教材土力学（成都水力发电学校李道荣）第168页例题：某均质粘性土坡，坡高为5m，坡比为1：2，土的有关指标为：，c=9.81kpa，试用圆

10、弧条分法求土坡稳定的最小安全系数kmin。解：由于计算量太大，为避免出错，教材中采用列表试算，得出其中一次计算的安全系数为ks1.55；现采用自编程序计算：得kmin1.524091。如图：【例2】大专教材土力学（黄河水利出版社务新超）第160页例题：一均质粘性土坡，高为20m，坡比为1：3，土的有关指标为：重度，粘聚力c=9.81kpa，内摩擦角，试用太沙基公式（即圆弧条分法）求土坡稳定的最小安全系数kmin。解：由于计算量太大，为避免出错，教材中采用列表的方法进行试算，得出其中一次计算的安全系数为ks1.83；现采用自编程序计算：得kmin1.4878（图略）。本例题也被收录在高等学校教材

11、土力学（中国水利水电出版社天津大学杨进良）第195页中。除边坡稳定分析程序外，本人还编制了赤平投影分析、抗剪强度指标计算、附加应力计算等工程地质与土力学方面的可视化实用小程序，为了方便各位专家和同仁共同研究或探讨，特附上这几个小程序，也可致电邮至：zhangh索取。下面列出圆弧条分法边坡稳定分析程序的主要程序代码，以供参考：option explicitprivate declaration as stringprivate sub bgvuvuti_click() form2.show form2.fontbold = false form2.fontsize = 2

12、0 form2.fontname = 隶书 form2.forecolor = rgb(50, 100, 200) form2.print 本软件适用于圆弧条分法边坡稳定分析 form2.print 如需示范演算,请直接单击按纽 form2.print 如需了解软件功能,请单击最大化按钮， form2.print 软件使用说明书便会呈现于窗体右侧 form2.fontsize = 15 form2.print copyrihgt : 张华庆 december 9th . 2002end subprivate sub ffliyi_click() h = 5: a = 26.565: v = 1

13、7.66: fia = 15: c = 9.81end subprivate sub fflioq_click() h = 25: a = 24.33: v = 17.66: fia = 17.5: c = 9.81end subprivate sub fflisf_click() h = 120: a = 30: v = 18: fia = 22.5: c = 10end subpublic sub form_load() show h = 20: a = 30 v = 18.78: fia = 22.58: c = 9.8 fi = fia / 57.2958 将值转化为弧度 t10 =

14、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) &

15、功能键及菜单说明： & chr(10) & 延时功能：f1f5（快慢）； & 分条因数：f8、f9、f11、f12（少多,即分条宽由宽到窄） & chr(10) & 渗透力：由于影响因素较多，仅考虑浸润面积为滑体面积的16、14、12、34等几种情形，作近似计算，且假定a1b1斜率即为渗流水力坡降。 & chr(10) & 滑体面积：自动计算并显示滑动面以上土体的面积。 & chr(13) & f5：运行（start） & space(2) & f10：暂停继续（breakcontinue） & chr(10) & ctrl+f6：设置速度（毫秒） label1.caption = decla

16、rationend subprivate sub a0_click() jt = 0end subprivate sub a2_click() jt = 0.2end subprivate sub a4_click() jt = 0.4end subprivate sub a6_click() jt = 0.6end subprivate sub a8_click() jt = 0.8end subprivate sub hta_click() dim s$ s = msgbox(滑体面积为 & str$(sa) & (平方米), 64, 渗流影响分析)end subprivate sub p

17、icture1_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 subprivate sub fprint_click() slptime = 0end sub 打印窗体暂时的不要，改为延时功能private sub cprint_click() slptime = 200end

18、 sub 打印窗体暂时的不要，改为延时功能private sub ft10_click(index as integer) t10 = 10end subprivate sub ft20_click(index as integer) t10 = 20end subprivate sub ft40_click(index as integer) t10 = 40end subprivate sub ft80_click(index as integer) t10 = 80end subprivate sub ft100_click(index as integer) t10 = 100end

19、subprivate sub menu1_click(index as integer) on error resume next h = inputbox(请输入土坡高度h（m）:, 输入对话框, 20, 8000, 7500)end subprivate sub pobi_click() dim pbn as single on error resume next pbn = inputbox(请输入坡比 1：n 中 n 的值:, 输入对话框, 2, 8000, 7500) a = 57.2957795 * atn(1 / pbn)end subprivate sub pojk_click

20、() on error resume next a = inputbox(请输入坡度角（）:, 输入对话框, 30, 8000, 7500)end subprivate sub menu3_click(index as integer) on error resume next v = inputbox(请输入土的容重:, 输入对话框, 18.78, 8000, 7500)end subprivate sub menu4_click(index as integer) on error resume next fia = inputbox(请输入滑面或土的摩擦角:, 输入对话框, 22.58,

21、 8000, 7500)end subprivate sub menu5_click(index as integer) on error resume next c = inputbox(请输入滑面或土的粘聚力c:, 输入对话框, 9.8, 8000, 7500) end subprivate sub speedset_click() on error resume next slptime = inputbox(请输入延时毫秒数:, 设置速度对话框, 1000, 8000, 7500)end subsub start_click() call clsn_click form1.label1

22、.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: era

23、se 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.setfocusend subsub seacher() *此过程用来即时显示每一试算圆心时的滑弧及分条图形* 笛卡尔坐标系：过图片框底部向右（）为x轴，过坡面中点向

24、上（）为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点坐标过程 form1.cls form1.picture

25、1.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 ub = 3.14159 - atn(abs(oy - b1y) / (ox - b1x) 若o点在b1点右侧,则ub为钝角end ifpicture1.line (ox - h / 12, oy)-(ox + h / 12, oy), rgb(120, 100

26、, 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 = 3picture1.line (ox, oy - h / 8)-(ox, f(ox), rgb(0, 200, 40) 画过滑弧圆心的铅垂点划线picture1.drawstyle = 0 doev

27、ents 当程序运行时允许转让控制权 sleep (slptime) 增加延时功能下面三条语句的作用是把已画过的圆心和滑弧及时抹掉 next oy next ox next b1xend subsub 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

28、.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 picture1.currentx = -2.9 * h: picture1.currenty = 5.9 * h picture1.print 容重=; v; kn/m3; space(20);

29、 摩擦角 ; 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) / 100picture1.line (ox - h / 12, oy)-(ox + h / 12, oy), rgb(0,

30、0, 0) 画过圆心的水平线 picture1.print opicture1.line (ox, oy - h / 12)-(ox, oy + h / 12), rgb(0, 0, 0) 画过圆心的铅直线 ua = atn(abs(oy - a1y) / (ox - a1x) if ox b1x then ub = 3.14159 - atn(abs(oy - b1y) / (ox - b1x) 若o点在b1点右侧,则ub为钝角 end ifpicture1.circle (ox, oy), r, rgb(250, 50, 140), -(3.1416 + ua), -(6.2832 - u

31、b) 画滑动面圆弧（含径向线,若不要径向线,请将-号去掉） picture1.currentx = a1x: picture1.currenty = a1y 注意：a1x、a1y若尚未赋值 picture1.print a1 则a1、b1是坐标错误的 picture1.currentx = b1x: picture1.currenty = b1y picture1.print b1for y1 = 0 to oy - h / 10 step oy / 20picture1.line (ox, y1)-(ox, y1 + r / 25) 画过滑弧圆心的铅垂(虚)线next y1end subsu

32、b 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 and oa2 r then if oav = r then a1x = avx: a1y = a1y exit for elseif oav r then a1x = avx: a2x = a2x: goto k end if end if next a1x elseif

33、oa r then for a1x = ax to bx step 0.4 a1x = a1x: a2x = a1x + 0.4q: 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 -

34、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 = avx: a2x = a2x: goto q end if end if next a1x end ifend subsub fpsplay() 此过程功能为分别在窗体和图片框中显示文字以下

35、语句在窗体中显示文字 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); ) for

36、m1.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=; fo

37、rmat(sa, 0.00) form1.print 分条数fts=; fts print if 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) *下

38、面语句圈定试算圆心范围* 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)-(b

39、x, by) 画坡底水平线 picture1.print b picture1.currentx = -2.9 * h: picture1.currenty = 5.9 * h picture1.print 容重=; v; kn/m3; 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); spac

40、e(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 b1end subsub countdarc() static i as integer i = 0 计算最危险滑弧圆心及半径开始 for b1x = bx - bx / 4 to bx + b

41、x / 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 = jv then exit sub end if next oy next ox next b1xend subprivate sub dangerarc_click() call countdarc call a1search 寻找a1点坐标开始-结束 form1.cls form1.picture1.cls call zhq if jt = 0.01 then print 总运算次数为jn; jn; 次, 算到其中第i= ; jv; 次时获得最危险滑弧 end if print 不考虑渗流影响时