最新圆弧条分法土坡稳定分析的VB电算解法518.docx
《最新圆弧条分法土坡稳定分析的VB电算解法518.docx》由会员分享,可在线阅读,更多相关《最新圆弧条分法土坡稳定分析的VB电算解法518.docx(31页珍藏版)》请在冰豆网上搜索。
最新圆弧条分法土坡稳定分析的VB电算解法518
圆弧条分法土坡稳定分析的VB电算解法5.18
圆弧条分法土坡稳定分析的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页中。
除边坡稳定分析程序外,本人还编制了赤平投影分析、抗剪强度指标计算、附加应力计算等工程地质与土力学方面的可视化实用小程序,为了方便各位专家和同仁共同研究或探讨,特附上这几个小程序,也可致电邮至:
索取。
下面列出圆弧条分法边坡稳定分析程序的主要程序代码,以供参考:
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'说明延时默认值为零
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
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
……"
lse
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.HtA.Enabled=True
Beep
DangerArc.SetFocus
EndSub
SubSeaCheR()'*****此过程用来即时显示每一试算圆心时的滑弧及分条图形*****
'笛卡尔坐标系:
过图片框底部向右(→)为X轴,过坡面中点向上(↑)为Y轴
j=0
ForB1x=Bx-Bx/4ToBx+Bx/4StepBx/8
IfB1xB1y=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
Picture1.Scale(-3*H,6*H)-(3*H,0)'第一次定义图片框尺寸比例
Callzhq'进入计算安全系数并绘制分条过程
CallFPsPlay'进入在窗体和图片框中显示文字及数据过程
Ua=Atn(Abs((Oy-A1y)/(Ox-A1x)))
IfOxUb=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)))
IfOxUb=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.Circle(Ox,Oy),R,RGB(250,50,140),-(3.1416+Ua),-(6.2832-Ub)'画滑动面圆弧(含径向线,若不要径向线,请将"-"号去掉)
Picture1.CurrentX=A1x:
Picture1.Curr