圆弧条分法土坡稳定分析的vb电算解法课件

圆弧条分法土坡稳定分析的 VB 电算解法 解题思路说明与讲解 云南省水利水电学校 张华庆 1 圆弧条分法土坡稳定分析的圆弧条分法土坡稳定分析的 VB 电算解法电算解法 云南省水利水电学校 张华庆 【摘要】目前，大多数版本的土力学教材在讲到均质粘性土边坡稳定性分析方法时，均采用 圆弧条分法，该法存在大量的作图及试算，若用人工求解，则不但会耗费大量宝贵的时间，而且 还难以得到精确解。现就介绍作者自编的粘性土边坡简单圆弧条分法的 VB 电算解法。 【关键词】土坡稳定 圆弧条分法 VB 电算 全自动 可视化 引言：迄今为止，求解均质粘性土边坡稳定最小安全系数 Ksmin，通常都采用瑞典 圆弧法或圆弧条分法，由于该法存在大量的手工作图以及滑弧试算，因此如果采用过 去传统的人工方法求解，则因为效率低下，不但浪费稿纸和笔墨、浪费大量宝贵的时 间，而且还难以得出较 为精确的解答。为此， 笔者利用 Visual Basic 6.0 编写和开发了可视 化的土坡稳定分析程序， 该程序及解法具有如下 特点：作图与试算全 过程可视化跟踪显示， 其优越性是过去的 DOS 计算程序所无法 比拟的；速度快，能在一分钟之内（一般仅需 10 秒左右）找到最危险滑弧以及土 坡稳定最小安全系数 Ksmin；充分体现圆弧条分法的精髓，精确度高；简单易行， 并且经过若干教材的例习题演算，证明结果一致；演算过程进行速度可以调节，便 2 于观察和教学；程序大小仅 72KB，但却功能完善，完全可以满足教学需要（当然 也完全可以用于工程计算） ；自动化程度高：仅需输入必要的坡高、坡比（或坡度 角） ，土的容重 以及土的抗剪强度参数 、c 的值，即可自动完成计算，因此为全自 动化设计，功能上明显优越于目前有的半自动化边坡稳定分析程序（如有的为 DOS 程序；有的虽然是可视化程序，但却是半自动化的，需要在 AutoCAD 里人为画出试 算滑弧，然后才能让机器去计算土坡稳定最小安全系数，这自然影响时效及精度，因 为你不可能有时间和耐心去画多达 100 条滑弧！就别说更多的） ；分条数可突 破人为影响及限制，例如可以达到 1000 条之多；无论何种条件下，当分条数逐渐 增多时，土坡稳定最小安全系数 Ksmin均能趋于一个稳定的极限值，证明计算方法以 及计算结果是正确的且是可靠的和稳定的；可以考虑渗透力的影响。 下面就本程序的编程思路作一些必要的说明和讲解，以便和有兴趣致力于此方面 研究的同仁共同探讨，并希望得到专家的指点，以便使程序变得更加完美和人性化。 一、滑弧圆心的确定变得简单化：试算法首先要确定的是滑弧的圆心，由于过去 是采用人工作图计算，工作量或者说劳动强度可想而知，于是必须尽可能地缩小滑动 圆心的范围，这样一来才需要定出 M 点坐标（4.5H，2H） 、确定 O 点所需的参数 125，235，以及利用十字交叉方法确定滑动圆心 Oi点的做法就是必需的 和可以理解的；但是对于计算机来说，它最不怕的就是计算，因而这种做法就显得多 余或者说不必要了，我们可以让它在平面上可能的区域（即所画圆弧至少有一端能交 于坡面、坡顶或坡底，而另一端也不至于悬空的点）在纵横两个方向每间隔一定的点 计算一次去寻找最危险滑弧，问题因此有了更为简单而圆满的解决办法。 二、分条面积计算简单化：由于分条宽度 b 是固定的，所以分条面积取决于分条 高度 hi，笔者采用解析几何方法加以处理： 3 列出坡面及坡顶、坡底的方程，并且列出圆弧的方程，再列出第 i 个土条的竖 线方程； 用公式法或试算法（迭代原理）求出交点坐标，再用两点间距离公式计算出 hi； 分条面积计算公式：Sihi*b (“*”号即“”号)； 由于分条数可以很多，当分条数多到一定程度（一般当取 b=R/100 或更大）时 滑体两端三角形面积部分对总面积的正负或加减影响可以不予考虑。 三、滑弧与坡面的交点分析：程序分别考虑了滑弧与土坡上游交点的两种情况， 即交于坡顶和交于上半坡面；滑弧与土坡下游交点的两种情况，即交于坡底和交于 下半坡面。这几种情况都是圆弧条分法所要求考虑的，均可由解析几何方法得出其交 点坐标，然后由上下两交点之间的水平间距 d 及 b=R/n 计算出分条数 fts。 四、滑弧半径及滑弧长的确定：在滑弧圆心确定后，滑弧半径通过两点间距离公 式确定，即；至于滑弧长，可通过公式计算滑 22 )()( 11 AOAO yyxxR*R L 弧总长，即在分条数 fts 确定后，水平方向正负两边圆心角即可用反三角函数算出， 然后，这样做的优点比起计算各分条长然后累加的办法来说有二：避免 21 了误差累积，提高了精确度；使计算过程由繁琐变为简单。 五、最小安全系数 Ks的得出原理：第二次计算得出的 Ks和第一次的进行比较， 若小于第一次的值则留下作为 Kmin的临时值赋给变量 Kt，否则用第一次的值作为临时 值，第三次计算得出的再和这个临时值亦即变量 Kt作比较，若小于则赋给变量 Kt，依 此类推，最后即可得出最小安全系数 Kmin。 六、渗透力计算原理：由于渗透力的计算比较麻烦（需要确定的边界条件较多） ， 故考虑作简单化处理，即只考虑两个因素：浸润面积占滑体总面积的百分比；作 4 用方向：一般按平行于坡面方向考虑；力的大小：由于渗透力是体积力，按单宽计 算就是面积力，即。SiSjJ w * 七、最小安全系数所对应圆心及滑弧的显示：有两种方法可以实现这一目的： 每次计算均用变量来存储 滑动圆心、滑弧半径及其 两端点的坐标，此法比较 麻烦且效率较低，影响计 算速度；总体计算完后， 再算一次。由于此次只需 计算至 KsKmin的那一次 即可结束，并且此次计算 过程不需要图形显示（只 需计算完成后显示最终结果） ，故计算速度非常快，一般不超过 1 至 10 秒。 八、延时显示功能：考虑到多媒体教学时的需要，编程时用了五种不同的速度来 执行计算过程，分别对应键盘上的 F1至 F5键，并且还可以自定义更慢的速度以方便 教学、讲解或观察，如左图和右图所示。 九、分条宽度 b=R/n：在人工试算时，一般取 n=10，但是计算机却可以取 n=10、20、40、80、100、甚至 n=1000 来计算 （自定义） ，例如上面左图是 n=100 的情形。 十、坐标原点的选取：为了能够比较居中地显 示图形，选取 y 轴通过滑面 AB 中点，x 轴通过滑 5 坡底面下距离滑坡底面 d=H*（1/3）取整数处。这样，当 H20m 时，滑面底部 B 点 纵坐标为 7m，而滑面顶、底部 A、B 两点的横坐标始终为负正相反数，如下图。 十一、交点坐标的得出：通过直线、圆的方程进行代数求解得出，当有两个解时， 由边界条件自动判断具体取哪一个值。 十二、安全系数计算公式： 滑体内无渗流通过时，计算公式为： )sin( )cos(tan )sin( )tancos( ii ii ii iii s hb Lchb W LcW K 滑体内有渗流通过时，计算公式为： )sin()sin( )cos(tan)cos(tan 21 21 iisatii iiii s hbhb Lchbhb K 以上是编程思路，下面就历界教材上的部分例习题作计算讲解及说明。 【例 1】中专教材土力学 （成都水力发电学校 李道荣）第 168 页例题：某均 质粘性土坡，坡高为 5m，坡比为 1：2，土的有关指标为： ，c=9.81kPa，试用圆弧条分法求土坡稳定的最小安全系数 3 /66.17mKN15 Kmin。 解：由于计算量太大，为避免出错，教材中采用列表试算，得出其中一次计算的 6 安全系数为 Ks1.55；现采用自编程序计算：得 Kmin1.524091。如图： 【例 2】大专教材土力学 （黄河水利出版社 务新超）第 160 页例题：一均质 粘性土坡，高为 20m，坡比为 1：3，土的有关指标为：重度，粘聚 3 /66.17mKN 力 c=9.81kPa，内摩擦角，试用太沙基公式（即圆弧条分法）求土坡稳定的最 20 小安全系数 Kmin。 解：由于计算量太大，为避免出错，教材中采用列表的方法进行试算，得出其中 7 一次计算的安全系数为 Ks1.83；现采用自编程序计算：得 Kmin1.4878（图略） 。 本例题也被收录在高等学校教材土力学 （中国水利水电出版社 天津大学 杨进良）第 195 页中。 除边坡稳定分析程序外，本人还编制了赤平投影分析、抗剪强度指标计算、附加 应力计算等工程地质与土力学方面的可视化实用小程序，为了方便各位专家和同仁共 同研究或探讨，特附上这几个小程序，也可致电邮至：ZhangH 索取。 下面列出圆弧条分法边坡稳定分析程序的主要程序代码，以供参考： 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 H = 20: a = 30 v = 18.78: fia = 22.58: c = 9.8 fi = fia / 57.2958 将 值转化为弧度 t10 = 10 给分条宽度 b=R/t10 中的 t10 赋一初值 8 jt = 0 给浸润面积系数 jt 赋一初值零（不考虑其影响） slptime = 0 说明延时默认值为零 Form1.DangerArc.Enabled = False Form1.DisPlay.Enabled = False Form1.HtA.Enabled = False Label1.FontSize = 15 Declaration = “ 软件使用说明书 “ v; “kn/m3“; Space(20); “摩擦角 “; CInt(fia * 100) / 100; “; Space(20); “粘聚力 c=“; c; “kPa“ 12 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 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 And OA2 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 A1x = Avx: A2x = A2x: GoTo q End If End If Next A1x End If End Sub Sub FPsPlay() 此过程功能为分别在窗体和图片框中显示文字 以下语句在窗体中显示文字 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“); 14 “)“ 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 “考虑渗流影响时:“ 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“ 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“ 15 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 = 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= “; jv; “次时获得最危险滑弧“ End If Print “不考虑渗流影响时“; 16 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); “,“; CSng(Oy); “)“; Tab(58); “坡度角 “; CInt(a * 100) / 100; “ Print “最危险滑弧的半径为 R=“; CSng(R); Tab(32); “米“, “最危险滑弧长为 L= “; CSng(SumLi); Tab(74); “米“ Print Print “分条宽 ftk=R/“; t10 Print “分条数 fts =“; fts; “条“ Print Form1.Print “Area=“; Sa Print If 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(Oy); “)“; Tab(46); “A1 点坐标为 (“; CSng(A1x); “,“; CSng(A1y); “)“ Form1.Print “坡高为 H = “; H; “米“; Tab(22); “坡度角为 “; CInt(5729.58 * Atn(Abs(kab) / 100; “; Tab(46); “B1 点坐标为 (“; CSng(B1x); “,“; CSng(B1y); “)“ If 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(“; 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 End Sub Private Sub ZDYWIDTH_Click() 18 On Error Resume Next t10 = InputBox(“请输入滑体分条宽度“ & Chr$(10) & “ b = R / t 中 t 的值:“, “分条宽度对 话框“, 200, 8000, 7500) End Sub Private Sub Suckset_Click() On Error Resume Next jt = InputBox(“请输入浸润面积滑体面积之比例系数“ & Chr$(10) & “如为 25则输入 0.25 等等“, “浸润面积对话框“, 0.25, 8000, 7500) 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