圆弧条分法土坡稳定分析的VB电算解法518.docx

资源描述

圆弧条分法土坡稳定分析的VB电算解法518.docx

《圆弧条分法土坡稳定分析的VB电算解法518.docx》由会员分享，可在线阅读，更多相关《圆弧条分法土坡稳定分析的VB电算解法518.docx（25页珍藏版）》请在冰豆网上搜索。

圆弧条分法土坡稳定分析的VB电算解法518.docx

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

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

──解题思路说明与讲解

云南省水利水电学校　张华庆

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

云南省水利水电学校　张华庆

【摘要】目前，大多数版本的土力学教材在讲到均质粘性土边坡稳定性分析方法时，均采用圆弧条分法，该法存在大量的作图及试算，若用人工求解，则不但会耗费大量宝贵的时间，而且还难以得到精确解。

现就介绍作者自编的粘性土边坡简单圆弧条分法的VB电算解法。

【关键词】土坡稳定　圆弧条分法　VB电算　全自动　可视化

引言：

迄今为止，求解均质粘性土边坡稳定最小安全系数Ksmin，通常都采用瑞典圆弧法或圆弧条分法，由于该法存在大量的手工作图以及滑弧试算，因此如果采用过去传统的人工方法求解，则因为效率低下，不但浪费稿纸和笔墨、浪费大量宝贵的时间，而且还难以得出较为精确的解答。

为此，笔者利用VisualBasic6.0编写和开发了可视化的土坡稳定分析程序，该程序及解法具有如下特点：

①作图与试算全过程可视化跟踪显示，其优越性是过去的DOS计算程序所无法比拟的；②速度快，能在一分钟之内（一般仅需10秒左右）找到最危险滑弧以及土坡稳定最小安全系数Ksmin；③充分体现圆弧条分法的精髓，精确度高；④简单易行，并且经过若干教材的例习题演算，证明结果一致；⑤演算过程进行速度可以调节，便于观察和教学；⑥程序大小仅72KB，但却功能完善，完全可以满足教学需要（当然也完全可以用于工程计算）；⑦自动化程度高：

仅需输入必要的坡高、坡比（或坡度角），土的容重γ以及土的抗剪强度参数φ、c的值，即可自动完成计算，因此为全自动化设计，功能上明显优越于目前有的半自动化边坡稳定分析程序（如有的为DOS程序；有的虽然是可视化程序，但却是半自动化的，需要在AutoCAD里人为画出试算滑弧，然后才能让机器去计算土坡稳定最小安全系数，这自然影响时效及精度，因为你不可能有时间和耐心去画多达100条滑弧！

——就别说更多的）；⑧分条数可突破人为影响及限制，例如可以达到1000条之多；⑨无论何种条件下，当分条数逐渐增多时，土坡稳定最小安全系数Ksmin均能趋于一个稳定的极限值，证明计算方法以及计算结果是正确的且是可靠的和稳定的；⑩可以考虑渗透力的影响。

下面就本程序的编程思路作一些必要的说明和讲解，以便和有兴趣致力于此方面研究的同仁共同探讨，并希望得到专家的指点，以便使程序变得更加完美和人性化。

一、滑弧圆心的确定变得简单化：

试算法首先要确定的是滑弧的圆心，由于过去是采用人工作图计算，工作量或者说劳动强度可想而知，于是必须尽可能地缩小滑动圆心的范围，这样一来才需要定出M点坐标（4.5H，2H）、确定O点所需的参数β1＝25°，β2＝35°，以及利用十字交叉方法确定滑动圆心Oi点的做法就是必需的和可以理解的；但是对于计算机来说，它最不怕的就是计算，因而这种做法就显得多余或者说不必要了，我们可以让它在平面上可能的区域（即所画圆弧至少有一端能交于坡面、坡顶或坡底，而另一端也不至于悬空的点）在纵横两个方向每间隔一定的点计算一次去寻找最危险滑弧，问题因此有了更为简单而圆满的解决办法。

二、分条面积计算简单化：

由于分条宽度b是固定的，所以分条面积取决于分条高度hi，笔者采用解析几何方法加以处理：

①列出坡面及坡顶、坡底的方程，并且列出圆弧的方程，再列出第i个土条的竖线方程；

②用公式法或试算法（迭代原理）求出交点坐标，再用两点间距离公式计算出hi；

③分条面积计算公式：

Si＝hi*b　（“*”号即“×”号）；

④由于分条数可以很多，当分条数多到一定程度（一般当取b=R/100或更大）时滑体两端三角形面积部分对总面积的正负或加减影响可以不予考虑。

三、滑弧与坡面的交点分析：

程序分别考虑了①滑弧与土坡上游交点的两种情况，即交于坡顶和交于上半坡面；②滑弧与土坡下游交点的两种情况，即交于坡底和交于下半坡面。

这几种情况都是圆弧条分法所要求考虑的，均可由解析几何方法得出其交点坐标，然后由上下两交点之间的水平间距d及b=R/n计算出分条数fts。

四、滑弧半径及滑弧长的确定：

在滑弧圆心确定后，滑弧半径通过两点间距离公式确定，即

；至于滑弧长，可通过公式

计算滑弧总长，即在分条数fts确定后，水平方向正负两边圆心角即可用反三角函数算出，然后

，这样做的优点比起计算各分条长然后累加的办法来说有二：

①避免了误差累积，提高了精确度；②使计算过程由繁琐变为简单。

五、最小安全系数Ks的得出原理：

第二次计算得出的Ks和第一次的进行比较，若小于第一次的值则留下作为Kmin的临时值赋给变量Kt，否则用第一次的值作为临时值，第三次计算得出的再和这个临时值亦即变量Kt作比较，若小于则赋给变量Kt，依此类推，最后即可得出最小安全系数Kmin。

六、渗透力计算原理：

由于渗透力的计算比较麻烦（需要确定的边界条件较多），故考虑作简单化处理，即只考虑两个因素：

①浸润面积占滑体总面积的百分比；②作用方向：

一般按平行于坡面方向考虑；③力的大小：

由于渗透力是体积力，按单宽计算就是面积力，即

。

七、最小安全系数所对应圆心及滑弧的显示：

有两种方法可以实现这一目的：

①每次计算均用变量来存储滑动圆心、滑弧半径及其两端点的坐标，此法比较麻烦且效率较低，影响计算速度；②总体计算完后，再算一次。

由于此次只需计算至Ks＝Kmin的那一次即可结束，并且此次计算过程

不需要图形显示（只需计算完成后显示最终结果），故计算速度非常快，一般不超过1至10秒。

八、延时显示功能：

考虑到多媒体教学时的需要，编程时用了五种不同的速度来执行计算过程，分别对应键盘上的F1至F5键，并且还可以自定义更慢的速度以方便教学、讲解或观察，如左图和右图所示。

九、分条宽度b=R/n：

在人工试算时，一般取n=10，但是计算机却可以取n=10、20、40、80、100、甚至n=1000来计算（自定义），例如上面左图是n=100的情形。

十、坐标原点的选取：

为了能够比较居中地显示图形，选取y轴通过滑面AB中点，x轴通过滑坡底面下距离滑坡底面d=H*（1/3）取整数处。

这样，当H＝20m时，滑面底部B点纵坐标为7m，而滑面顶、底部A、B两点的横坐标始终为负正相反数，如下图。

十一、交点坐标的得出：

通过直线、圆的方程进行代数求解得出，当有两个解时，由边界条件自动判断具体取哪一个值。

十二、安全系数计算公式：

①滑体内无渗流通过时，计算公式为：

②滑体内有渗流通过时，计算公式为：

以上是编程思路，下面就历界教材上的部分例习题作计算讲解及说明。

【例1】中专教材《土力学》（成都水力发电学校　李道荣）第168页例题：

某均质粘性土坡，坡高为5m，坡比为1：

2，土的有关指标为：

，c=9.81kPa，

，试用圆弧条分法求土坡稳定的最小安全系数Kmin。

解：

由于计算量太大，为避免出错，教材中采用列表试算，得出其中一次计算的安全系数为Ks＝1.55；现采用自编程序计算：

得Kmin＝1.524091。

如图：

【例2】大专教材《土力学》（黄河水利出版社　务新超）第160页例题：

一均质粘性土坡，高为20m，坡比为1：

3，土的有关指标为：

重度

，粘聚力c=9.81kPa，内摩擦角

，试用太沙基公式（即圆弧条分法）求土坡稳定的最小安全系数Kmin。

解：

由于计算量太大，为避免出错，教材中采用列表的方法进行试算，得出其中一次计算的安全系数为Ks＝1.83；现采用自编程序计算：

得Kmin＝1.4878（图略）。

本例题也被收录在高等学校教材《土力学》（中国水利水电出版社　天津大学　杨进良）第195页中。

除边坡稳定分析程序外，本人还编制了赤平投影分析、抗剪强度指标计算、附加应力计算等工程地质与土力学方面的可视化实用小程序，为了方便各位专家和同仁共同研究或探讨，特附上这几个小程序，也可致电邮至：

ZhangHqing@索取。

下面列出圆弧条分法边坡稳定分析程序的主要程序代码，以供参考：

OptionExplicit

PrivateDeclarationAsString

PrivateSubbgvuvuti_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:

张华庆December9th.2002"

EndSub

PrivateSubffliyi_Click（）

H=5:

a=26.565:

v=17.66:

fia=15:

c=9.81

EndSub

PrivateSubfflioq_Click（）

H=25:

a=24.33:

v=17.66:

fia=17.5:

c=9.81

EndSub

PrivateSubfflisf_Click（）

H=120:

a=30:

v=18:

fia=22.5:

c=10

EndSub

PublicSubForm_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赋一初值

jt=0'给浸润面积系数jt赋一初值零（不考虑其影响）

slptime=0'说明延时默认值为零

Form1.DangerArc.Enabled=False

Form1.DisPlay.Enabled=False

Form1.HtA.Enabled=False

Label1.FontSize=15

Declaration="软件使用说明书"&Chr（10）&"著作人声明：

"&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

EndSub

PrivateSubA0_Click（）

jt=0

EndSub

PrivateSubA2_Click（）

jt=0.2

EndSub

PrivateSubA4_Click（）

jt=0.4

EndSub

PrivateSubA6_Click（）

jt=0.6

EndSub

PrivateSubA8_Click（）

jt=0.8

EndSub

PrivateSubHtA_Click（）

DimS$

S=MsgBox（"滑体面积为"&Str$（Sa）&"（平方米）",64,"渗流影响分析"）

EndSub

PrivateSubPicture1_MouseDown（ButtonAsInteger,ShiftAsInteger,xAsSingle,YAsSingle）

UnloadForm2'此语句相当于IfButton=1or2ThenUnloadForm2

IfButton=2ThenForm1.PopupMenuLPRINT'弹出式菜单应用例子

EndSub

PrivateSubFprint_Click（）

slptime=0

EndSub'打印窗体暂时的不要，改为延时功能

PrivateSubCprint_Click（）

slptime=200

EndSub'打印窗体暂时的不要，改为延时功能

PrivateSubft10_Click（IndexAsInteger）

t10=10

EndSub

PrivateSubft20_Click（IndexAsInteger）

t10=20

EndSub

PrivateSubft40_Click（IndexAsInteger）

t10=40

EndSub

PrivateSubft80_Click（IndexAsInteger）

t10=80

EndSub

PrivateSubft100_Click（IndexAsInteger）

t10=100

EndSub

PrivateSubmenu1_Click（IndexAsInteger）

OnErrorResumeNext

H=InputBox（"请输入土坡高度H（m）:

","输入对话框",20,8000,7500）

EndSub

PrivateSubPobi_Click（）

DimpbnAsSingle

OnErrorResumeNext

pbn=InputBox（"请输入坡比1：

n中n的值:

","输入对话框",2,8000,7500）

a=57.2957795*Atn（1/pbn）

EndSub

PrivateSubPojk_Click（）

OnErrorResumeNext

a=InputBox（"请输入坡度角α（°）:

","输入对话框",30,8000,7500）

EndSub

PrivateSubmenu3_Click（IndexAsInteger）

OnErrorResumeNext

v=InputBox（"请输入土的容重γ:

","输入对话框",18.78,8000,7500）

EndSub

PrivateSubmenu4_Click（IndexAsInteger）

OnErrorResumeNext

fia=InputBox（"请输入滑面或土的摩擦角φ:

","输入对话框",22.58,8000,7500）

EndSub

PrivateSubmenu5_Click（IndexAsInteger）

OnErrorResumeNext

c=InputBox（"请输入滑面或土的粘聚力c:

","输入对话框",9.8,8000,7500）

EndSub

PrivateSubSpeedSet_Click（）

OnErrorResumeNext

slptime=InputBox（"请输入延时毫秒数:

","设置速度对话框",1000,8000,7500）

EndSub

SubStart_Click（）

CallClsn_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）

EraseKs:

EraseKn

Kmin=1000'为Kmin赋一初值1000，以便比较

Kjin=1000'为Kjin赋一初值1000，以便比较

CallSeaCheR'调用过程SeaCheR，留意参数的传送问题

jn=j

Form1.DangerArc.Enabled=True

Form1.HtA.Enabled=True

Form1.Label1.Caption=Declaration

Beep

DangerArc.SetFocus

EndSub

SubSeaCheR（）'*****此过程用来即时显示每一试算圆心时的滑弧及分条图形*****

'笛卡尔坐标系：

过图片框底部向右（→）为X轴，过坡面中点向上（↑）为Y轴

j=0

ForB1x=Bx-Bx/4ToBx+Bx/4StepBx/8

IfB1x

B1y=g（B1x）

ElseIfBx<=B1xThen

B1y=By

EndIf

ForOx=AxToBx+H/2StepL/10

ForOy=Ay+H/4ToAy+4.5*HStepL/10

j=j+1'此句从zhq中提前后便修正了Ksi显示中的错误

CallA1Search'进入用试算法求得A1点坐标过程

Form1.Cls

Form1.Picture1.Cls

Picture1.Scale（-3*H,6*H）-（3*H,0）'第一次定义图片框尺寸比例

Callzhq'进入计算安全系数并绘制分条过程

CallFPsPlay'进入在窗体和图片框中显示文字及数据过程

Ua=Atn（Abs（（Oy-A1y）/（Ox-A1x）））

IfOx

Ub=Atn（Abs（（Oy-B1y）/（Ox-B1x）））

ElseIfOx=B1xThen

Ub=3.14159/2'若O点在B1点之上,则Ub=90°（直角）

ElseIfOx>B1xThen

Ub=3.14159-Atn（Abs（（Oy-B1y）/（Ox-B1x）））'若O点在B1点右侧,则Ub为钝角

EndIf

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）'增加延时功能^下面三条语句的作用是把已画过的圆心和滑弧及时抹掉

NextOy

NextOx

NextB1x

EndSub

SubDrawPicture（）

Picture1.Scale（-3*H,6*H）-（3*H,0）

'*************下面语句圈定试算圆心范围****************

DimIxAsSingle

DimIyAsSingle

ForIx=AxToBx+H/2StepL/10

ForIy=Ay+H/4ToAy+4.5*HStepL/10

Picture1.PSet（Ix,Iy）

NextIy

NextIx

'----------------------------------------------------------------------

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/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）））

IfOx

Ub=Atn（Abs（（Oy-B1y）/（Ox-B1x）））

ElseIfOx=B1xThen

Ub=3.14159/2'若O点在B1点之上,则U

展开阅读全文