第11章 程 序 设 计
11,1 前言
从本书介绍的内容可以发现随着近代数值分析技术的进展边坡稳定分析的极限平衡分析方法已经发展到十分成熟的程度但是随之而出现的问题是这一分析步骤必须通过计算机程序才能得以实现而现状是尽管计算机的硬件已经获得了飞速的发展尽管土力学的理论和实现这些理论的工具和手段均已十分先进但岩土工程许多应用软件仍然处于起步状态远不能满足工程实际的要求且大部分岩土工程师对边坡稳定分析计算的实际应用水平还不很高使用近代科学技术的成果来提高岩土工程设计的效率和水平是工程师的强烈愿望也是保证建筑物的安全和减少投资非常重要手段当前岩土工程软件的落后状况是有多方面因素造成的
首先开发岩土工程软件需要投入大量的人力和物力研究一个新的数值分析方法并编制相应的科研程序是一回事将其变为一个能够在工程设计中广泛应用的软件是另一回事其源代码的数量起码增加十倍在Windows的界面上编制一个直观方便的具有处理图形对话和纠错功能的软件不仅是源代码数量上的增加而且意味着对编程人员素质的新要求而这一要求是由所支付的高额报酬体现在程序价格中的
其次岩土工程软件的市场极为有限对于这种专业性的软件单靠市场回报来实现其自身的发展是不现实的而这一问题并未得到共识众所周知设计人员在图纸上多画一根预应力锚索多画一根抗滑桩则是十万甚至百万元的投资但对于开发这些有关软件的研究经费(尽管无法与施工费用相比)却难以得到批准
再次我们的软件管理还远未达到规范化市场化的水平盗版现象普遍存在因此营造一个保护软件的法律和道德约束环境是软件作为一个产业兴旺发达的最基本的条件
工程应用软件开发的另一个重要问题是确认其可靠性和准确性工程软件应用的对象是关系到人民生命财产安全的建筑物因此绝非儿戏对程序进行鉴定和考核自然是一个必须屡行的手续但笔者认为对一个包含成千上万行指命的又需要面临复杂多变的实际情况的工程应用软件单靠鉴定和考核是不可能真正发现其全部问题的这样的程序只有在长期的大量的实际应用过程中才能变得成熟成为可以信赖的设计工具
作者希望通过第11章和第12章的内容推动边坡稳定计算机程序的普及和推广本章将介绍编制一个边坡稳定分析程序的核心步骤和源程序将这几个子程序串联起来就是一个计算任意几何形状边坡包含第2章和第3章介绍的各种分析方法的程序本章还将讨论与程序推广使用相关的问题在第12章作者提供STAB程序的一个简化版使用这一版本读者可以进行本书介绍的部分稳定分析工作需要说明的是本书所列的所有源程序
336 土质边坡稳定分析?原理
方法
程序
只供科研和教学使用作者不对其承担法律责任他人也不得将其应用于商业目的
11,2 几何图形的识别和分析
11,2,1 边坡剖面图形的数字信息
通常坝体并非均一而是由若干个物理力学指标不同的土质区域组成图11.1所示坝体由6个土质不同的区域组成控制点17个边界线22条图中滑弧上阴影所示的土条中心线和3条边界线相交将该土条分成了三段其高度分别为H
1
H
2
和H
3
并分别位于
I区V区和VI区所有这些在图解法中依靠人的直观操作是很容易判断和量测得的编制程序的关键是如何让计算机代替人判断土条的中心线究竟与哪几个线条相交各段高度分别为多少我们的办法是让计算机对所有的边界线循环一遍逐个判断每一条线是否与土条中心线相交边界线用所压土层号的数组来表示输入这些数据如交上则记下该交点的X Y坐标值及该线段所压的土层号这样根据相邻两交点的Y坐标值的差值即可算得H
1
H
2
和H
3
根据土层编号即可查得相应的物理力学指标随即可算得土条重量地震力孔隙水压力等这一对几何图形进行处理的方法充分利用了计算机可进行大量重复的判断忌讳过多的逻辑判断的特点为编制一个效率高鲁棒性 (Robust) 强的边坡稳定分析程序创造了条件
图 11,1 由6个土质不同的区域组成坝体
11,2,2 构造圆弧滑裂面的步骤
设某一坝坡坡面为具有n个折点的折线如图11.2所示(其中起点1和尾点n在极远处)
建立坐标系统X轴水平指向滑坡方向为正Y坐标轴垂直向下为正某一半径为R
0
的滑弧圆心O′ 至坝坡各折点距离分别为R
1
,R
2
,…,R
n
该圆弧自左至右在连接i i+1两点的线段进入坝体上交点在连接j j+1线段逸出坝体下交点
不难证明某一半径为R
0
的圆弧与线段A
1
A
2
有一个且只有一个交点的必要充分条件是
R
2
R
0
R
1
其中R
1
R
2
为圆心至线段端点的距离如果R
0
同时大于R
1
R
2
则圆弧与A
1
A
2
没有交点如果R
0
同时小于R
1
R
2
则圆弧与A
1
A
2
不是没有交点就是有两个交点,见图第11章 程序设计 337
11.3
则根据上述原理可知i j应满足下面条件
图 11,3 具有n个折点的边坡线与圆弧滑裂面相交的关系
图 11,2 圆弧滑裂面与边坡相交的关系
(11.1)
ii
RRR <≤
+ 01
(11.2)
10 +
<≤
jj
RRR
找到了上下交点所在的线段号i j之后就可以根据该线段两端点的坐标值和圆弧的圆心坐标值x
0
y
0
和半径值R算得其交点的坐标值对于点m (m = i,j)计算公式为
(11.3)
22
0
2
0
)()( Ryyxx =?+?
(11.4)
sxxyy
mm
= )/()(
(11.5)
)/()(
11 mmmm
xxyys=
++
联立以上方程组求解得
])()1([
)1(
1
2
00
2
02
TsxRssTx
s
x+±+
+
=
(11.6)
(11.7)
0
yysxT
mm
+?=
此联立方程组有两个根计算下交点时取+号计算上交点时取?号
11,2,3 构造任意形状滑裂面的步骤
如果要将若干点用光滑的曲线相连可以采用样条函数予以实现在STAB程序中采用了第二种边界条件的三次样条函数
第二种边界条件三次样条函数曲线两端点的二阶导数值为已知值即)(
112
xyy ′′=
)(
2 nn
xyy ′′=在构造滑裂面时取0)(
112
=′′= xyy 0)(
2
=′′=
nn
xyy由此确定的样条称为
338 土质边坡稳定分析?原理
方法
程序
自然样条根据此边界条件计算n个节点上的一阶导数值m nj
j
,,2,1,L=相应的计算公式如下
(11.8) 5.0
1
=a
)(2
)()(
3
12
12
1
xx
xyxy
b
×= (11.9)
(11.10)
1,,2,1
1
=?=
+
njxxh
jjj
L
)(
11 jjjj
hhh +=
α
(11.11)
])())(1[(3
111 jjjjjjjjj
hyyhyy?+=
+
ααβ (11.12)
])1(2[
1?
+?=
jjjj
aa αα
(11.13)
1,,3,2])1(2[])1([
11
=?+=
njabb
j
j
jjjj
Lααβ (11.14)
)2/(])(3[
1111
+=
nnnnnn
abhyym (11.15)
(11.16)
1,,2,1
1
L=+=
+
nnjbmam
jjjj
于是对于,利用数值导数 (j )
1+
<≤
jj
xxx
j
m n,,2,1L=计算插值节点处函数值的公式如下
1
3
3
2
2
3
1
3
2
1
2
1
3
3
2
2
3
1
3
2
1
2
)(
2
)(
1
)(
1
)(
1
)(
2
)(
3
)(
2
)(
3
)(
+++
+++
+
+
=
jj
j
j
j
fjj
j
j
j
f
jj
j
j
j
jj
j
j
j
mxx
h
xx
h
hmxx
h
xx
h
h
yxx
h
xx
h
yxx
h
xx
h
xy
(11.17)
如果插值区间被等距划分为l段),(
1+jj
xx则位于第s (0<s<l) 段的被插值节点
jj
hlsxx )(+=处的函数值的计算公式如下
(11.18)
)]([)(
1234
κκκκ
sjsjsj
hhhxy +++=
其中
lshxxh
jjsj
=?= )(
(11.19)
(11.20)
jjjjj
hmmxyxy )()]()([2
111 ++
++?=κ
(11.21)
jjjjj
hmmxyxy )2()]()([3
112 ++
+=κ
(11.22)
jj
hm?=
3
κ
(11.23))(
4 j
xy=κ
第11章 程序设计 339
11,3 边坡稳定分析的源程序
11,3,1 概述
本章介绍边坡稳定分析程序STAB的核心步骤和程序程序 CIRCLE.FOR具体实施
11.2.2节的计算步骤对给定的边坡剖面信息构造一个圆弧滑裂面程序SPLINE.FOR具体实施11.2.3节的计算步骤对给定的滑面控制点信息构造一个任意形状滑裂面程序
PROFILE.FOR具体实施11.2.1节的计算步骤 计算土条重量和物理力学参数程序S.FOR
计算安全系数包含了第2章和第3章介绍的各种分析方法的源程序即极限平衡法的通用条分法瑞典法毕肖普法陆军工程师团法和简化法将上述几个子程序串联起来就是一个计算复杂几何形状边坡稳定的安全系数的程序
在介绍这些程序时将配以图11.4所示的一个例子这一例子是将在11.4节中介绍的澳大利亚ACADS程序考核的EX1(c)该例已在第2章和第4章中多次使用读者在移植这些程序时可通过这些考核题进行考核
图 11,4 滑裂面和条分
11,3,2 构造圆弧滑裂面的程序 CIRCLE.FOR
1,说明
本子程序对给定的边坡剖面信息构造一个圆弧滑裂面(图11.4)通过子程序SEAR找到圆弧滑裂面与边坡的上下交点所在线段子程序DIVI计算上下交点的坐标并划分土条
2,源程序
REAL YTENSION,DS,CX,CY,XN(99),YN(99),X(80),Y(80),ALF(80)
INTEGER IN1,LSL,N,LNUM(80),IC(80,3),IWR5,L,NN,NLOW,NUPP
OPEN(5,FILE='',STATUS='UNKNOWN')
OPEN(6,FILE='',STATUS='UNKNOWN')
NLOW=0
NUPP=0
READ(5,*)IWR5,YTENSION,DS,CX,CY
C IWR5=0,1分别相应顶部有无拉力缝YTENSION=拉力缝底Y坐标
C DS,CX,CY=滑弧深度圆心坐标
READ(5,*)N!土条总数
READ(5,*)IN1!线段总数
DO 302 I=1,IN1
READ(5,*)(IC(I,J),J=1,3)!线段信息
302 CONTINUE
340 土质边坡稳定分析?原理
方法
程序
READ(5,*)LSL!外边坡线总数
READ(5,*)(LNUM(I),I=1,LSL)!外边坡线编号
READ(5,*)NN!点总数
READ(5,*)(L,XN(I),YN(I),I=1,NN)!点坐标
WRITE(6,703)
703 FORMAT(T25,'******************'//)
WRITE(6,704)
704 FORMAT(10X,'THEABSCISSAVALUESOFTHEUPPERANDLOWER')
CALL SEAR(LNUM,IC,DS,CX,CY,XN,YN,LSL,NLOW,NUPP)!寻找上下交点所在线段
CALL DIVI(N,IC,IWR5,DS,CX,CY,XN,YN确定上下交点坐标条分
$,YTENSION,NLOW,NUPP,X,Y,ALF)
WRITE(6,711)
711 FORMAT(T5,'NO.',T17,'X',T32,'Y',T47,'ALF')
DO 309 I=1,N
WRITE(6,710)I,X(I),Y(I),ALF(I)
710 FORMAT(1X,I5,3F15.6)
309 CONTINUE
12 CALLCLOSEFL
END
SUBROUTINE CLOSEFL
CLOSE(5)
CLOSE(6)
RETURN
END
SUBROUTINE DIVI(N,IC,IWR5,DS,CX,CY,XN,YN
$,YTENSION,NLOW,NUPP,X,Y,ALF)
INTEGERN S,N,K,IC(80,3),K1,KTEMP,IWR5,J,J1,JTEMP,NLOW,NUPP
REAL RO,DS,CY,XN(80),YN(80),SK,T,XL,CX,XU1,YTENSION
REAL YU1,SJ,XU,ROO,YU,X2(80),Y2(80),X(80),Y(80),ALF(80)
S(E1,F1,E2,F2)=(F2-F1)/(E2-E1)
NS=N-1
RO=DS-CY
K=IC(NLOW,1)
K1=IC(NLOW,2)
IF(XN(K).GT.XN(K1))GOTO190
KTEMP=K
K=K1
K1=KTEMP
190 IF(ABS(XN(K)-XN(K1)).LT.0.00005)GOTO23
SK=S(XN(K),YN(K),XN(K1),YN(K1))
T=SK*XN(K)-YN(K)+CY
XL=(CX+SK*T+SQRT((1+SK**2)*RO**2-(SK*CX-T)**2))/(1+SK**2)
!按式(11.6)计算下交点的横坐标
GOTO28
23 XL=XN(K)
28 CONTINUE
IF(IWR5.NE.0)THEN
XU1=CX-SQRT(RO*RO-(YTENSION-CY)**2)
YU1=CY+SQRT(RO**2-ABS(XU1-CX)**2)
ENDIF
J=IC(NUPP,1)
J1=IC(NUPP,2)
IF(XN(J).GT.XN(J1))GOTO191
JTEMP=J
J=J1
J1=JTEMP
第11章 程序设计 341
191 F(ABS(XN(J)-XN(J1)).LT.0.00005)GOTO82
SJ=S(XN(J),YN(J),XN(J1),YN(J1))
T=SJ*XN(J)-YN(J)+CY
XU=(CX+SJ*T-SQRT((1+SJ**2)*RO**2-(SJ*CX-T)**2))/(1+SJ**2)
!XU上交点的横坐标
ROO=RO**2-ABS(XU-CX)**2
IF(ROO.LT.0)RETURN!无交点,需修改RO,重新计算
YU=CY+SQRT(RO**2-ABS(XU-CX)**2)
IF(CY.GT.YN(J).AND.CY.GT.YN(J1))RETURN!凸弧,否定,重算
IF(IWR5.NE.0.AND.YU.LT.YU1)XU=XU1
GOTO20
82 XU=XN(J)
20 X2(1)=XU
DO541I=1,NS!划分土条,求各点坐标
X2(I)=(I-1)*(XL-XU)/(NS-1)+XU!土条宽度B=(XL-XU)/(NS-1)
Y2(I)=CY+SQRT(RO**2-ABS(X2(I)-CX)**2)
541 CONTINUE
IF(ABS(Y2(NS)-Y2(1)).LT.0.01)THEN
WRITE(6,*)'THE SLIP SURFACE INTERCEPTS A HORIZONTAL SURFACE'
RETURN
ENDIF
N=NS+1
X2(N)=X2(NS)
Y2(N)=Y2(NS)
C 计算土条中点坐标
DO544I=2,NS
X(I)=(X2(I)+X2(I-1))/2,
Y(I)=(Y2(I)+Y2(I-1))/2,
ALF(I)=ATAN((Y2(I)-Y2(I-1))/(X2(I)-X2(I-1)))
544 CONTINUE
ALF(1)=ALF(2)
ALF(N)=ALF(NS)
X(N)=X2(N)
Y(N)=Y2(N)
X(1)=X2(1)
Y(1)=Y2(1)
RETURN
END
SUBROUTINE SEAR(LNUM,IC,DS,CX,CY,XN,YN,LSL,NLOW,NUPP)
REAL R,DS,CY,R2,CX,XN(80),YN(80),R1,D,XXL,SLO
INTEGER I,I1,K,LNUM(80),J1,J2,IC(80,3),IN,NLOW,NUPP,LSL
I=1 !NLOW=下交点所在直线段的线段号
I1=0 !NUPP=上交点所在直线段的线段号
12K=LNUM(I)
J1=IC(K,1)
J2=IC(K,2)
R=DS-CY!R=圆弧半径
R2=SQRT(ABS(CX-XN(J2))**2+ABS(CY-YN(J2))**2)
R1=SQRT(ABS(CX-XN(J1))**2+ABS(CY-YN(J1))**2)
!R1,R2为圆心至该线段两端点的距离
IF(R.EQ.R2.AND.I1.NE.0)GOTO11
IF(R.GT.R1.AND.R.GT.R2)GOTO11 !无交点
IF(R.LT.R1.AND.R.LT.R2)GOTO19 !有两个或零个交点
GOTO20
19 CALL DD(XN(J1),YN(J1),XN(J2),YN(J2),CX,CY,D,IN)
XXL=XN(J2)-XN(J1)
342 土质边坡稳定分析?原理
方法
程序
IF(ABS(XXL).LT.0.0001)THEN
SLO=1.0E-6
ELSE
SLO=(YN(J2)-YN(J1))/XXL !斜率
ENDIF
IF(IN.EQ.1.OR.SLO.LE.0)GOTO11
C 第一条外边坡线为水平,同时圆弧与其有两个交点,不接受,继续搜索
IF(R.LT.D.OR.R.EQ.D)GOTO11 !无交点
NLOW=LNUM(I)
NUPP=NLOW
RETURN
20 IF(I1.NE.0)GOTO14
NLOW=LNUM(I)
I1=I1+1
11 CONTINUE
IF(I.EQ.1.AND.R.LT.R1.AND.R.GT.R2)RETURN
I=I+1
IF(I.GT.LSL)GOTO14
GOTO12
14 NUPP=LNUM(I)
15 CONTINUE
RETURN
END
SUBROUTINE DD(X1,Y1,X2,Y2,XC,YC,D,IN)
REAL X1,Y1,X2,Y2,XC,YC,D
REAL XXL,AK,AKF,X,Y,XT,YT
INTEGER IN
D=9999,
IN=0
XXL=X2-X1
IF(ABS(XXL).LT.0.001)XXL=1.0E-6
IF(ABS((Y2-Y1)/XXL).LT.0.001)RETURN
AK=(Y2-Y1)/XXL
AKF=1./AK
X=(AK*X1-Y1+AKF*XC+YC)/(AK+AKF)
Y=AK*X-AK*X1+Y1 !从圆心作垂直与线段的线,获与线段的交点
XT=(X-X1)*(X-X2)
YT=(Y-Y1)*(Y-Y2)
IF(XT.LT.0.OR.YT.LT.0)GOTO12
IN=1 !交点不在线段内
RETURN
12 D=SQRT((X-XC)**2+(Y-YC)**2) !线段与圆心的垂直距离
RETURN
END
3,例题
数据文件CIRCLE.DAT
0 0,-24.235 -35.164 -42.046
21
9
1 2 3
2 3 1
3 4 1
2 5 3
5 6 2
6 7 2
第11章 程序设计 343
7 8 2
5 9 3
9 10 3
3
1 2 3
10
1,000 -25.000
2 -30.000 -25.000
3 -50.000 -35.000
4 -70.000 -35.000
5 -40.000 -27.000
6 -50.000 -29.000
7 -54.000 -31.000
8 -70.000 -31.000
9 -52.000 -24.000
10 -70.000 -24.000
计算成果CIRCLE.C
THE ABSCISSA VALUES OF THE UPPER AND LOWER
NO,X Y ALF
1 -51.522050 -35.000000 1.094647
2 -50.955680 -33.901810 1.094647
3 -49.822940 -31.978910,969015
4 -48.690200 -30.491210,863838
5 -47.557460 -29.278430,770531
6 -46.424730 -28.265840,685123
7 -45.291990 -27.411070,605365
8 -44.159260 -26.687410,529825
9 -43.026520 -26.076840,457513
10 -41.893780 -25.566700,387709
11 -40.761050 -25.147810,319849
12 -39.628310 -24.813490,253490
13 -38.495570 -24.558860,188260
14 -37.362830 -24.380470,123833
15 -36.230100 -24.275990,059923
16 -35.097360 -24.244130 -.003743
17 -33.964630 -24.284500 -.067425
18 -32.831890 -24.397580 -.131382
19 -31.699150 -24.584810 -.195887
20 -30.566420 -24.848610 -.261233
21 -30.000050 -25.000020 -.261233
11,3,3 构造任意形状滑裂面的程序SPLINE.FOR
1,说明
本程序刘德贵1980用来计算已知第二种边界条件(曲线两端点处的二阶导数值)
三次样条函数的一阶导数值并利用样条函数对一元函数进行插值进而构筑包含直线段的复合滑裂面对于滑裂面上控制点区间为直线段时则直接线性插值而对于非直线段区间则运用上述介绍的方法进行三次样条插值滑裂面被NS1个点分为NS1?1段从上交点向下交点编号为1,2,,NS1?1程序用LS代表此NS1?1段中为直线段的线段总数在数组
LSA(I)(I=1,2,,LS)中存入这些直线段的编号
LS,滑裂面上NS1?1条线段中为直线段的总数
LSA(LS)滑裂面上线段为直线段的编号
344 土质边坡稳定分析?原理
方法
程序
AI(I,4)滑裂面上控制点处四个与样条插值有关的数
1
,κ
2
,κ
3
,κ
4
κ
A(I)滑裂面上控制点处的一阶导数值
j
m
B(I)一维数组工作单元
2,源程序
C SUBROUTINE SP1(LS,LSA,X2,Y2)
DIMENSION KQ2(50),X1(50),Y1(50),AI(50,4),LSA(60),A(50),B(50)
DIMENSION X2(80),Y2(80)
COMMON /A23/NS1,KQ2,X1,Y1/A13/IDD,IWR1,IWR2,IWR3
COMMON/WALL/IWALL,GWALL,HMW,EWALL,ETA
INTEGER*4 OPTION(6),OP1(6)
COMMON/OPP/OPTION,OP1
OPEN(5,FILE=' ',STATUS='UNKNOWN')
OPEN(6,FILE=' ',STATUS='UNKNOWN')
READ(5,*)NS1
DO I=1,NS1
READ(5,*)KQ2(I),X1(I),Y1(I)
ENDDO
C 如果滑裂面上控制点总数NS1=2,则按直线段处理并计算滑裂面与土条侧面边界线交点的x,y值;
READ(5,*)LS
DO I=1,NS1-1
LSA(I)=0
ENDDO
IF(LS.NE.0)READ(5,*)(LSA(I),I=1,LS)
LSA(LS+1)=NS1
C 滑裂面上第NS1个控制点对应的土条侧面边界线的编号
C X1(I) Y1(I)滑裂面上第I个控制点的x y值
N1=0
DO 100 IL=1,LS+1
L1=N1+1
N1=LSA(IL)
N2=N1-1
IF(LS.EQ.0)GOTO 83 ! 滑裂面上没有直线段
IF(N1.GE.NS1) GOTO 40
AI(N1,1)=0,! 计算直线段左端点处κ
1
,κ
2
,κ
3
,κ
4
AI(N1,2)=0,
AI(N1,3)=Y1(N1+1)-Y1(N1)
AI(N1,4)=Y1(N1)
40 IF(N1-L1.EQ.0)GOTO 100 ! 滑裂面上没有直线段
83 B(L1)=1.5*(Y1(L1+1)-Y1(L1))/(X1(L1+1)-X1(L1))
A(L1)=-0.5 ! 插值区间左端点处的一阶导数值
CN=3.0*(Y1(N1)-Y1(N2))/(X1(N1)-X1(N2))
IF(N1-L1.EQ.1)GOTO 55 ! 插值区间只具有左右两个端点
L2=L1+1
DO 50 J=L2,N2
H1=X1(J)-X1(J-1)
H=X1(J+1)-X1(J)
AF=H1/(H+H1) ! 计算α值
BT=3.*((1.-AF)*(Y1(J)-Y1(J-1))/H1+AF*(Y1(J+1)-Y1(J))/H)
!计算β值
A(J)=-AF/(2.+(1.-AF)*A(J-1)) ! 计算插值区间内插值节点处的一阶导数值m
j
50 B(J)=(BT-(1.-AF)*B(J-1))/(2.+(1.-AF)*A(J-1))
55 A(N1)=(CN-B(N2))/(2.+A(N2))
60 A(N2)=A(N2)*A(N2+1)+B(N2) ! 计算插值区间右端点处的一阶导数值m
n
第11章 程序设计 345
N2=N2-1
IF(N2.GE.L1)GOTO 60
N2=N1-1
DO 80 I=L1,N2
H=X1(I+1)-X1(I)
AI(I,4)=Y1(I) ! 计算插值区间内插值节点处的κ
1
,κ
2
,κ
3
,κ
4
AI(I,3)=A(I)*H
AI(I,2)=3.*(Y1(I+1)-Y1(I))-(2.*A(I)+A(I+1))*H
80 AI(I,1)=2.*(Y1(I)-Y1(I+1))+(A(I)+A(I+1))*H
100 CONTINUE
DO 200 I=1,NS1-1
L1=KQ2(I)
N1=KQ2(I+1)
STEP=1./(N1-L1)
SX=X1(I+1)-X1(I)
S1=0,
DO 200 J=L1,N1
X2(J)=X1(I)+S1*SX ! 计算滑裂面与土条侧面边界线交点的X Y值
Y2(J)=AI(I,4)+S1*(AI(I,3)+S1*(AI(I,2)+S1*AI(I,1)))
200 S1=S1+STEP
NW=KQ2(NS1)
WRITE(6,*)'N= ',NW
WRITE(6,105)(M,X2(M),Y2(M),M=1,NW)
105 FORMAT(1X,I5,2F10.3)
STOP
END
3,光滑曲线滑裂面例题
(1) 数据文件Spline-1.dat
5
1 -51.522 -35.000
6 -46.314 -28.235
10 -41.685 -25.511
15 -35.898 -24.224
20 -30.000 -25.000
0
计算成果Spline-1.c
N= 20
1 -51.522 -35.000
2 -50.480 -33.466
3 -49.439 -31.978
4 -48.397 -30.580
5 -47.356 -29.317
6 -46.314 -28.235
7 -45.157 -27.280
8 -44.000 -26.545
9 -42.842 -25.975
10 -41.685 -25.511
11 -40.528 -25.107
12 -39.370 -24.765
13 -38.213 -24.495
14 -37.055 -24.311
15 -35.898 -24.224
16 -34.718 -24.243
17 -33.539 -24.352
18 -32.359 -24.530
19 -31.180 -24.754
20 -30.000 -25.000
346 土质边坡稳定分析?原理
方法
程序
4,曲线-直线组合滑裂面例题
(1) 数据文件Spline-2.dat
5
1 -51.522 -35.000
6 -46.314 -28.235
10 -41.685 -25.511
15 -35.898 -24.224
20 -30.000 -25.000
1
3
(2) 计算成果Spline-2.c
N= 20
1 -51.522 -35.000
2 -50.480 -33.459
3 -49.439 -31.965
4 -48.397 -30.565
5 -47.356 -29.306
6 -46.314 -28.235
7 -45.157 -27.300
8 -44.000 -26.583
9 -42.842 -26.011
10 -41.685 -25.511
11 -40.528 -25.511
12 -39.370 -25.511
13 -38.213 -25.511
14 -37.055 -25.511
15 -35.898 -25.511
16 -34.718 -25.409
17 -33.539 -25.307
18 -32.359 -25.204
19 -31.180 -25.102
20 -30.000 -25.000
11,3,4 计算土条重量和物理力学参数的程序 PROFILE.FOR
1,说明
本子程序根据11.2.1介绍的方法计算每个土条的重量条底的强度指标和孔隙水压力
2,源程序
DIMENSION SH(20),M4(20),SHNE(20),M4NE(20)
DIMENSION XU(80),HSL(80),X3(80),Y3(80)
COMMON/CRU/RUS(10),LINES(80),YSF(80)/B3/IKD
COMMON/A1/IN,IN1,PF(20),PC(20),PDW(20),PDS(20),IPH,ICPH(40,3)
COMMON/A2/NN,NJ(20),NHY,DHRP,WL(80),BDA(80),CA(80),FAA(80)
#,LRU(20)
COMMON/A3/RW,RU1,UWL/SWEAK/IWEAK,KQ3(50),LWK,IC3(40,3)
COMMON/A6/X(80),Y(80),XN(99),YN(99),X2(80),Y2(80),IC(80,3),N
COMMON/A7/W(80),RU(80),C(80),F(80)/A8/MYK,KXSX,MTE,KL
COMMON/A15/UWL1,NCOUN,MLA(80),PF1(20),PC1(20),MID(20),PNOR(20)
#,NHT,MRD
INTEGER*4 OPTION(6),OP1(6)
COMMON/OPP/OPTION,OP1
OPEN(5,FILE=' ',STATUS='UNKNOWN')
OPEN(6,FILE=' ',STATUS='UNKNOWN')
C IWR5=0,1 分别相应顶部有无拉力缝YTENSION=拉力缝底Y坐标
C DS,CX,CY=滑弧深度圆心坐标
第11章 程序设计 347
READ(5,*) N ! 土条总数
READ(5,*) RW,UWL ! 水容重坡外水位
READ(5,*) IN1 ! 线段总数
DO 302 I=1,IN1
READ(5,*) (IC(I,J),J=1,3) !线段信息
302 CONTINUE
READ(5,*)NN ! 点总数
READ(5,*) (L,XN(I),YN(I),I=1,NN) ! 点坐标
DO I=1,N
READ(5,*) I,X(I),Y(I)
ENDDO
READ(5,*)IPH
IF(IPH.EQ.0)GOTO 18
DO I=1,IPH
READ(5,*) (ICPH(I,J),J=1,3) !线段信息
ENDDO
18 READ(5,*)IN,(L,PF(I),PC(I),PF1(I),PC1(I),PDW(I),PDS(I),
%LRU(I),I=1,IN)
CONTINUE
C 对土条信息重新处理保证IC(I,1)的x坐标比IC(I,2)的小
DO 37 I=1,IN1
M1=IC(I,1)
M2=IC(I,2)
IF(XN(M1).LT.XN(M2))GO TO 37
IC(I,1)=IC(I,2)
IC(I,2)=M1
37 CONTINUE
IF(IPH.EQ.0)GO TO 672
C REARRANGE ICPH(I,J) IN THE SAME WAY AS IC(I,J)
DO 38 I=1,IPH
M1=ICPH(I,1)
M2=ICPH(I,2)
IF(XN(M1).LT.XN(M2))GO TO 38
ICPH(I,1)=ICPH(I,2)
ICPH(I,2)=M1
38 CONTINUE
672 CONTINUE
GTENSION=0
C FIND THE LINES WHICH INTERCEPT THE CONSIDERED SLICE
MX=IC(1,1)
YYY=YN(MX)
DO 813 I=1,IN1
DO 813 II=1,2
IY=IC(I,II)
IF(YN(IY).GT.YYY)YYY=YN(IY)
813 CONTINUE
YYY=YYY+100,
X3(1)=X2(1)
Y3(1)=Y2(1)
DO 211 K=1,N
MF=0
QN=0,
W(K)=0,
DHR=0
WL(K)=YYY
348 土质边坡稳定分析?原理
方法
程序
LPHS=IC(IN1,3)
IF(IPH.EQ.0)GO TO 739
DO 735 I1=1,IPH
M1=ICPH(I1,1)
M2=ICPH(I1,2)
IF(X(K).LT.XN(M1))GO TO 735
IF(X(K).GE.XN(M2))GO TO 735
IF(ABS(XN(M2)-XN(M1)).LT.0.001)GOTO 735
WL(K)=YN(M1)+(X(K)-XN(M1))*(YN(M2)-YN(M1))/(XN(M2)-XN(M1))
LPHS=ICPH(I1,3)
735 CONTINUE
739 CONTINUE
PP=0,
CT=0,
FT=0,
K1=0
DO 207 I1=1,IN1
M1=IC(I1,1)
M2=IC(I1,2)
M3=I1
IF(X(K).LT.XN(M1))GO TO 207
IF(X(K).GE.XN(M2))GO TO 207
SH2=YN(M1)+(X(K)-XN(M1))*(YN(M2)-YN(M1))/(XN(M2)-XN(M1))
152 FORMAT(1X,'**',I5,2F15.6)
DSH2=ABS((Y(K)-SH2))
DSH1=SH2-Y(K)
C ASSURE THAT THE STRENGTH PARAMETERS OF THE UNDERLYING SOIL
IF(DSH2.GT.0.00005.AND.DSH1.GT.0.001)GO TO 207
K1=K1+1
SH(K1)=SH2
M4(K1)=M3
207 CONTINUE
IF(K1.NE.1)GO TO 707
LX=M4(1)
GO TO 708
707 LX=M4(1)
HN=SH(1)
DO 45 L1=2,K1
IF(SH(L1).GT.HN)GO TO 45
LX=M4(L1)
HN=SH(L1)
45 CONTINUE
708 M1=IC(LX,1)
M2=IC(LX,2)
DELX=XN(M2)-XN(M1)
IF(ABS(DELX).GT.0.001)BBK=(YN(M2)-YN(M1))/(XN(M2)-XN(M1))
X3(K)=X2(K)
Y3(K)=YN(M1)+(X3(K)-XN(M1))*BBK
139 FORMAT(1X,'K1,LX,M2,M1,XM2,XM1',4I5,2F10.3)
LINES(K)=LX !96.12.29
47 DO 46 L1=1,K1
LSB=M4(L1)
M4(L1)=IC(LSB,3)
46 CONTINUE
IF(WL(K).GT.Y(K))GO TO 734
DO 530 ISP=1,K1
SDIF=ABS(WL(K)-SH(ISP))
IF(SDIF.LT.0.00001)GO TO 734
第11章 程序设计 349
530 CONTINUE
K1=K1+1
SH(K1)=WL(K)
M4(K1)=LPHS
734 J=1
27 MSP=1
DO 209 I=1,K1
IF(J.EQ.I)GO TO 209
IF(SH(J).LT.SH(I))GO TO 209
MSP=MSP+1
209 CONTINUE
SHNE(MSP)=SH(J)
M4NE(MSP)=M4(J)
IF(J.GE.K1)GO TO 229
J=J+1
GO TO 27
229 K1=K1+1
SHNE(K1)=Y(K)
YSF(K)=SHNE(1)
M4NE(K1)=M4NE(K1-1)
SH3=SHNE(K1)-SHNE(1)
IF((K1-1).LE.0)GO TO 641
DDH=0,
IF(SHNE(1).GT.UWL)W(K)=RW*(SHNE(1)-UWL)
DO 208 I=2,K1
M5=M4NE(I-1)
WL1=WL(K)-0.0001
IF(SHNE(I-1).LT.WL1)PDNO=PDW(M5)
IF(SHNE(I-1).GE.WL1)PDNO=PDS(M5)
DH=SHNE(I)-SHNE(I-1)
DWK=DH*PDNO
W(K)=W(K)+DWK
42 CT=CT+DH*PC(M5)
IF(PC(M5).LT.0.00001.AND.PF(M5).LT.0.00001)GO TO 589
FT=FT+DH*PF(M5)
GO TO 208
589 DDH=DDH+DH
DHR=DHR+DH
208 CONTINUE
DHRP=SH3-DHR
M7=M4NE(K1)
IF(MF.NE.0)M7=IC3(MF,3)
MLA(K)=M7
PP=0
XU(K)=SHNE(1)
HSL(K)=SH3
IF(IPH.EQ.0)GO TO 736
IF(Y(K).LT.WL(K))GOTO 590
WLM=WL(K)
IF(K.EQ.1.OR.K.EQ.N)GO TO 590
IF(LRU(M7).EQ.0.AND.SHNE(1).LE.UWL)PP=(Y(K)-WLM)*RW
IF(LRU(M7).EQ.0.AND.SHNE(1).GT.UWL)PP=(Y(K)-UWL)*RW
590 CONTINUE
WT1=Y(K)-UWL
IF(WT1.LT.0)WT1=0
W(K)=W(K)-WT1*RW
PP=PP-WT1*RW
736 IF(ABS(SH3).LT.0.001.OR.ABS(W(K)).LT.0.001)GO TO 641
350 土质边坡稳定分析?原理
方法
程序
IF(LRU(M7).GE.0)THEN
RU(K)=PP/W(K)
GOTO 961
ENDIF
IF(Y(K).GT.WL(K).OR.IPH.EQ.0)THEN
RU(K)=RUS(M7)
ELSE
RU(K)=0
ENDIF
961 BDA(K)=W(K)/SH3
CA(K)=CT/SH3
IF(ABS(SH3-DDH).LT.0.0001)GO TO 587
FAA(K)=FT/(SH3-DDH)
GO TO 586
587 FAA(K)=0,
586 CONTINUE
GOTO 55
641 RU(K)=0,
W(K)=0
BDA(K)=0,
CA(K)=0,
FAA(K)=0,
55 F(K)=PF(M7)
C(K)=PC(M7)
21 CONTINUE
211 CONTINUE
IF(ABS(W(N)).GT.0.001) GO TO 256
W(N)=0.001
100 FORMAT(I3/(7F11.6))
256 CONTINUE
C(1)=C(2)
F(1)=F(2)
265 WRITE(6,251)
252 FORMAT(1X,'NUMBER OF SLICES=',I5)
WRITE(6,252)N
WRITE(6,258)
251 FORMAT(1X,'DATA FOR EACH SLICE:')
258 FORMAT(T5,'NO,',T14,'X',T24,'Y',T34,'W',T44,'RU',T54,'C',
@T64,'F',T74,'BDA',T84,'WL')
WRITE(6,253)(I,MLA(I),X(I),Y(I),W(I),RU(I),C(I),F(I),BDA(I),WL(I)
%,I=1,N)
253 FORMAT(1X,2I5,8F10.3)
104 FORMAT(1X,'PORE PRESSURE',F10.3)
WRITE(6,254)
254 FORMAT(1X,'NOTE:W=WEIGHT PER UNIT WIDTH,RU=COEF,OF PORE PRES.')
WRITE(6,255)
255 FORMAT(6X,'C=COEFFICIENT OF COHESION,F=COEFFICIENT OF FRICTION')
WRITE(6,259)
259 FORMAT(6X,'BDA=AVERAGE UNIT WEIGHT,WL=Y VALUE OF PHREATIC LINE')
END
3,例题
(1) 数据文件PROFILE.DAT
21
9.8,0.0
9
1 2 3
2 3 1
第11章 程序设计 351
3 4 1
2 5 3
5 6 2
6 7 2
7 8 2
5 9 3
9 10 3
10
1,000 -25.000
2 -30.000 -25.000
3 -50.000 -35.000
4 -70.000 -35.000
5 -40.000 -27.000
6 -50.000 -29.000
7 -54.000 -31.000
8 -70.000 -31.000
9 -52.000 -24.000
10 -70.000 -24.000
1 -51.522050 -35.000000
2 -50.955680 -33.901810
3 -49.822940 -31.978910
4 -48.690200 -30.491210
5 -47.557460 -29.278430
6 -46.424730 -28.265840
7 -45.291990 -27.411070
8 -44.159260 -26.687410
9 -43.026520 -26.076840
10 -41.893780 -25.566700
11 -40.761050 -25.147810
12 -39.628310 -24.813490
13 -38.495570 -24.558860
14 -37.362830 -24.380470
15 -36.230100 -24.275990
16 -35.097360 -24.244130
17 -33.964630 -24.284500
18 -32.831890 -24.397580
19 -31.699150 -24.584810
20 -30.566420 -24.848610
21 -30.000050 -25.000020
0
3
1 38.000,000,000,000 19.500 19.500 0
2 23.000 5.300,000,000 19.500 19.500 0
3 20.000 7.200,000,000 19.500 19.500 0
(2) 计算成果PROFILE.C
DATA FOR EACH SLICE,
NUMBER OF SLICES= 21
NO,SOIL X Y W RU C F BDA WL
1 1 -51.522 -35.000,000,000,000 38.000,000 76.000
2 1 -50.956 -33.902 21.415,000,000 38.000 19.500 76.000
3 1 -49.823 -31.979 57.185,000,000 38.000 19.500 76.000
4 1 -48.690 -30.491 75.151,000,000 38.000 19.500 76.000
5 1 -47.557 -29.278 87.756,000,000 38.000 19.500 76.000
6 2 -46.425 -28.266 96.457,000 5.300 23.000 19.500 76.000
7 2 -45.292 -27.411 102.081,000 5.300 23.000 19.500 76.000
8 2 -44.159 -26.687 105.148,000 5.300 23.000 19.500 76.000
352 土质边坡稳定分析?原理
方法
程序
9 3 -43.027 -26.077 106.010,000 7.200 20.000 19.500 76.000
10 3 -41.894 -25.567 104.914,000 7.200 20.000 19.500 76.000
11 3 -40.761 -25.148 102.038,000 7.200 20.000 19.500 76.000
12 3 -39.628 -24.813 97.513,000 7.200 20.000 19.500 76.000
13 3 -38.496 -24.559 91.434,000 7.200 20.000 19.500 76.000
14 3 -37.363 -24.380 83.868,000 7.200 20.000 19.500 76.000
15 3 -36.230 -24.276 74.862,000 7.200 20.000 19.500 76.000
16 3 -35.097 -24.244 64.439,000 7.200 20.000 19.500 76.000
17 3 -33.965 -24.285 52.607,000 7.200 20.000 19.500 76.000
18 3 -32.832 -24.398 39.358,000 7.200 20.000 19.500 76.000
19 3 -31.699 -24.585 24.663,000 7.200 20.000 19.500 76.000
20 3 -30.566 -24.849 8.475,000 7.200 20.000 19.500 76.000
21 3 -30.000 -25.000,001,000 7.200 20.000,000 76.000
NOTE:W=WEIGHT PER UNIT WIDTH,RU=COEF,OF PORE PRES,
C=COEFFICIENT OF COHESION,F=COEFFICIENT OF FRICTION
BDA=AVERAGE UNIT WEIGHT,WL=Y VALUE OF PHREATIC LINE
11,3,5 计算安全系数的程序S.FOR
1,说明
本程序包括一个主程序和三个子程序子程序MP实现垂直条分法的计算过程包括陆军工程师团法子程序SAFE2实现瑞典法毕肖普法和简化法的计算过程这三个方法和相应计算公式已在第3章中曾作了详细介绍通过一个整型变量OPTION来控制计算方法通过变量ID?C来控制滑面形状(圆弧或非圆弧)对非圆弧滑裂面通过IWALL=0和=1
分别实现垂直条分法计算安全系数和土压力的功能本源程序中将水容重RW设为9.8也就是要求使用kN?m的量纲体系控制滑面几何形状的点有两套一套为界面在土条底的坐标用X2(I) Y2(I)代表土条总数为NS另一套为土条底中点的坐标用X(I) Y(I)来代表在输入W(I) RU(I)等一系列物理力学参数时都是相应土条中点的但是加上了左右两个端点的相应物理量因此读入了总数为N=NS+1行数据这些信息实际上就是程序PROFILE.FOR的输出成果在使用程序MP实现垂直条分法时程序首先使用子程序
SAFE2中简化法的计算安全系数作为初值然后进行第3节介绍的稳定分析方法的计算子程序SAFE2在进行Bishop法计算时同样用简化法计算安全系数初值读者还可从子程序MP中发现通过式(2.43)计算简化法1即λ初值的过程
2,源程序
REAL AB1
INTEGER N,NS,OPTION
COMMON/AS/ASP
COMMON/CIRCLE/ID_C,CX,CY,DS
C OPTION=0,SPENCER法
C OPTION=1,BISHOP法
C OPTION=2,瑞典法
C OPTION=3,工程师团法
C OPTION=4,简化法
C ID_C=0,任意形状滑裂面ID_C不为零圆弧滑裂面
REAL X2(80),Y2(80),W(80),X(80),Y(80),C(80),F(80),WL(80)
$,FXO(80),FX(80),BN(80)
REALBDA(80),FQUH(80),FDIS(80),YSF(80),RU(80),ALF(80),RW
WRITE(6,*)'输入数据文件名'
OPEN(5,FILE='',STATUS='UNKNOWN')
第11章 程序设计 353
WRITE(6,*)'输入计算成果文件名'
OPEN(6,FILE='',STATUS='UNKNOWN')
RW=9.8
AB1=0.0
WRITE(6,*)'边坡稳定分析'
READ(5,*)OPTION,ID_C
IF(ID_C.NE.0)READ(5,*)CX,CY,DS
READ(5,*)N!土条总数
NS=N-1
READ(5,*)(X(I),I=1,N)
READ(5,*)(Y(I),I=1,N)
READ(5,*)(ALF(I),I=1,N)
N=NS+1
DO302I=1,N
READ(5,*)W(I),RU(I),C(I),F(I),BDA(I),WL(I),FQUH(I),FDIS(I),YSF(I)
F(I)=TAN(F(I)*3.14159/180.)
302 CONTINUE
WRITE(6,*)(ALF(I),I=1,N)
!读入两端点和土条底中点的数据,其中涉及力的物理量均以单宽土条计
!X2,Y2=坐标,W=土条重量,RU=孔压系数,C=粘聚力=,F=摩擦系数,BDA=土容重,
!WL=浸润线Y坐标,FQUH,FDIS,YSF=水平地震力大小,Y坐标,垂直水平地震力大小
IF(OPTION.EQ.3)THEN
READ(5,*)ASP
WRITE(6,*)'陆军工程师团法平均坝坡=',ASP
ASP=ASP*3.14159/180,
ENDIF
WRITE(6,709)RW
709 FORMAT(1X,'水容重=',F10.3//)
WRITE(6,702)
702 FORMAT(T25,'DATAFORASSUMEDSIDEFORCEFUNCTION')
WRITE(6,703)
703 FORMAT(T25,'************************************'//)
!计算土条底中点的坐标和条底倾角
DO208I=2,NS
X2(I)=(X(I)+X(I+1))/2,
Y2(I)=(Y(I)+Y(I+1))/2,
208 CONTINUE
X2(N)=X(N)
Y2(N)=Y(N)
X2(1)=X(1)
Y2(1)=Y(1)
WRITE(6,716)
716 FORMAT(T5,'NO.',T14,'X',T24,'Y',T34,'W',T44,'RU',T54,'C',
# T64,'F',T74,'BDA',T84,'WL',T94,'QH',T104,'YE')
DO901I=1,N
WRITE(6,718)I,X(I),Y(I),W(I),RU(I),C(I),F(I),BDA(I),WL(I)
#,FQUH(I),FDIS(I)
901 CONTINUE
718 FORMAT(1X,I5,11F10.3)
WRITE(6,719)
719 FORMAT(1X,/'NOTE:')
WRITE(6,720)
720 FORMAT(5X,'W=土条的单宽重量')
WRITE(6,721)
721 FORMAT(5X,'RU=孔隙水压力系数')
WRITE(6,722)
722 FORMAT(5X,'C=粘聚力')
354 土质边坡稳定分析?原理
方法
程序
WRITE(6,723)
723 FORMAT(5X,'F=摩擦系数')
WRITE(6,724)
724 FORMAT(5X,'BDA=土条的平均容重')
WRITE(6,725)
725 FORMAT(5X,'WL=浸润线的Y坐标')
C WRITE(6,701)
C WRITE(6,*)'迭代过程'
SELECTCASE(OPTION)
CASE(0)
WRITE(6,701)
701 FORMAT(/10X,'SPENCERMETHOD')
CALLSAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,4)
IWALL=0
DO301I=1,80
FXO(I)=0
FX(I)=1
301 CONTINUE!TAN(BETA(I))=FXO(I)+ALAM*FX(I)
KXYX=0
CALLMP(N,KXYX,KS,AB1,DQ,DM,BN,FX,FXO,
& C,F,X2,Y2,X,Y,W,RU,FQUH,WL,FDIS,RW,DB,DF)
CASE(1)
WRITE(6,704)
704 FORMAT(/20X,'BISHOP法'/)
CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION)
CASE(2)
WRITE(6,705)
705 FORMAT(/20X,'瑞典法'/)
CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION)
CASE(3)
WRITE(6,731)
731 FORMAT(/10X,'工程师团法')
CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,4)
IWALL=0
DO311I=1,80
FXO(I)=0
FX(I)=1
311 CONTINUE!TAN(BETA(I))=FXO(I)+ALAM*FX(I)
KXYX=1
CALL MP(N,KXYX,KS,AB1,DQ,DM,BN,FX,FXO,
& C,F,X2,Y2,X,Y,W,RU,FQUH,WL,FDIS,RW,DB,DF)
CASE(4)
WRITE(6,789)
789 FORMAT(/20X,'简化法法'/)
CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION)
CASE DEFAULT
ENDSELECT
END
SUBROUTINE SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,ID)
REAL X2(80),F(80),FQUH(80),W(80),RU(80),ALF(80),X(80),C(80)
REAL AB1,FDIS(80)
REAL RO1,DM,RM,BX1,CDN,DG(80),G3,DX1,G,GX
REAL O,S1,FRIC,S2,S,DRM,P1,P2,DDM,DBX1,DCDN
COMMON/CIRCLE/ID_C,CX,CY,DS
INTEGER NS
IWALL=0
第11章 程序设计 355
RO1=1.0
NS=N-1
KL=2
DM=0,
RM=0,
BX1=0,
CDN=0,
DG(2)=0.
G3=0,
DO93I=2,NS
DX1=X2(I)-X2(I-1)
DG(I+1)=DG(I)+(F(I)-F(I+1))*ALF(I)
93G3=G3+(F(I)*ALF(I)+DG(I))*DX1
G3=G3/(X(N)-X(1))
DO202I=2,NS
IF(ID_C.EQ.0)THEN
RO1=1
ELSE
RO1=(FDIS(I)-CY)/(DS-CY)
ENDIF
G=ALF(I)
GX=G*F(I)-G3+DG(I)
DX1=X2(I)-X2(I-1)
O=FQUH(I)*SIN(G)
S1=W(I)*(COS(G)-RU(I)/COS(G))
FRIC=F(I)
S2=(S1-O)*FRIC
S=C(I)/COS(G)
DRM=S+S2
DRM=DRM*DX1
P1=W(I)*SIN(G)
P2=FQUH(I)*RO1
DDM=(P1+P2)*DX1
DM=DM+DDM
RM=RM+DRM
DBX1=-DRM-SIN(G)*GX*W(I)*DX1
BX1=BX1+DBX1
DCDN=DRM*GX
CDN=CDN+DCDN
202 CONTINUE
AB1=-BX1/DM+CDN/BX1!通过简化法求得AB1的迭代初值
IF(ID.EQ.4)WRITE(6,*)'简化法安全系数=',AB1
IF(ID.EQ.2)AB1=RM/DM
IF(ID.EQ.2)WRITE(6,*)'瑞典法安全系数=',AB1
IF(ID.NE.1)RETURN
WRITE(6,*)'迭代过程'
LM=1
41DT=0,
DM=0,
RM=0,
DO21I=2,NS
G=ALF(I)
DX1=X2(I)-X2(I-1)
S1=W(I)*(1-RU(I))*F(I)
S2=S1+C(I)
S=COS(G)*(1+F(I)*TAN(G)/AB1)
DRM=DX1*S2/S
RM=RM+DRM
356 土质边坡稳定分析?原理
方法
程序
P1=W(I)*SIN(G)
P2=FQUH(I)*RO1
DDM=(P1+P2)*DX1
DM=DM+DDM
DT=DT+DRM*TAN(G)*F(I)/(AB1**2+AB1*F(I)*TAN(G))
21 CONTINUE
AB2=RM/DM!BISHOP方法
IF(AB2.LT.0.OR.AB2.GE.15)THEN
WRITE(6,*)'F.SNOTREASONABLE=',AB2
AB2=9999.9
RETURN
ENDIF
WRITE(6,740)LM,AB2
740 FORMAT(1X,I5,F10.5)
IF(ABS(AB2-AB1).LT.0.0001)GOTO49
IF(LM.LT.10)GOTO9191
WRITE(6,9192)
9192 FORMAT(1X,'ITERATIONTROUBLESINBISHOPMETHOD')
GOTO49
9191 CK=1-DT/DM!牛顿迭代法
AB1=AB1-(AB1-AB2)/CK
LM=LM+1
GOTO41
49 WRITE(6,203)AB1
203 FORMAT(/5X,'BISHOP法安全系数=',F10.4)
RETURN
END
SUBROUTINE DETE(AB1,RW,DQ,DM,ALAM,KS,BN,FX,FXO,
& C,F,X2,Y2,W,RU,FQUH,YSF,WL,FDIS,N,DB,KXYX,DF)
REAL AB1,DQ,DT,DM,ALAM,DF
REAL C(80),F(80),X2(80),Y2(80),W(80),RU(80),FXO(80),FX(80)
REAL FQUH(80),YSF(80),WL(80),FDIS(80),BN(80)
INTEGER KS,N,KXYX
REAL X(80),ALF(80),G,Y(80)
REAL EXT,EXR,EYR,EMT,EMR,FMT,DMT,DEG,DX1
REAL DMG,BS,R1X,TX,QX,PT112,R1XT,POX
REAL POXT,POY,DPT112,DPX112,BSG,H3,YT(80),GF(80)
REAL GFX(80),GFY(80),BET(80),EFX(80)
REAL G1,TA,TA1,OP,FI,FI1,CE,OQ,OQ1,THO,CW
REAL SEGB,SEGB1,FI0,AXS11,OTEM,DDQ,PX112,DBS
REAL SX,SXG,DTX1,DTX2,DTX,TX1,DDP,DBS1,DBS2,EPX112
REAL AAKX1,DR1X1,DR1X2,DR1X,DR1XT
REAL DG,DDT,GAVE,BETA,DGX,OX,OY,HSL
REAL TAA,OQS2,OQS1,DG1,BT(80),EXR1,EXT1,EYR1
REAL EYT1,XCE,YCE,XRE,YRE,DR2X11,DR2X12,DR2X1
REAL GS,DPOY1,DPOY2,OY1,DDM,DDMT,PM,PMT,DDEG,DDMG
REAL OQ0,SEGB0,DGG,DMG1,DEG1,OXX
REAL TEMP1,TEMP2,GG,DDM2,DM1,EWALL,GWALL,PX
REAL DQM,DMM,DTM,DMTM,DEGM,DMGM,S1,HMW,DDM1,ETA
REAL GTENSION,HITE,WP,RW
INTEGER EYT,IWALL,NS
!本程序计算MORGENSTERN-PRICE,SPENCER和工程师团法
IWALL=0.!土压力的作用位置
HMW=0.
GTENSION=0,
HITE=0,
第11章 程序设计 357
GWALL=0.!土压力的大小
!GWALL=1计算土压力设定条件下的安全系数
!GWALL=2计算主动土压力
IF(KS.EQ.0)THEN!KS调用本程序段标志符,=0表示第1次调用
WRITE(6,705)
705 FORMAT(1X,/'VALUEOFFX0(I)INORDEROFI')
WRITE(6,706)(FXO(I),I=1,N)
706 FORMAT(1X,5F10.5)
WRITE(6,707)
707 FORMAT(1X,/'VALUEOFFX(I)INORDEROFI')
WRITE(6,708)(FX(I),I=1,N)
708 FORMAT(1X,5F10.5)
WRITE(6,726)
726 FORMAT(1X,/'迭代过程:'/)
WRITE(6,712)
712 FORMAT(T10,'NO.',T24,'DQ',T35,'DM',T55,'AB1',T65,'ALAM')
ENDIF
KS=1
DM=0,
DT=0.
EXT=GTENSION
EYT=0
EXR=0,
EYR=0,
EMT=-GTENSION*HITE/3,
EMR=0,
FMT=0,
DQ=0.
DMT=0,
DEG=0.
DMG=0,
BS=0,
R1X=0,
TX=0,
QX=0,
PT112=0,
PX112=0,
R1XT=0,
POX=0,
POXT=0.
POY=0,
DPT112=0,
DPX112=0,
BSG=0,
H3=HITE/3,
YT(1)=Y2(1)-H3
DQ=DQ-GTENSION
GF(1)=-DQ
DM=DM-GTENSION*H3!96.8.31
NS=N-1
DO303I=2,NS
X(I)=(X2(I)+X2(I-1))/2,
Y(I)=(Y2(I)+Y2(I-1))/2,
ALF(I)=ATAN((Y2(I)-Y2(I-1))/(X2(I)-X2(I-1)))
303 CONTINUE
ALF(1)=ALF(2)
ALF(N)=0,
358 土质边坡稳定分析?原理
方法
程序
X(N)=X2(NS)
Y(N)=Y2(NS)
X(1)=X2(1) !X,Y各土条的中点坐标(在滑裂面上)
Y(1)=Y2(1)
DO304I=1,N
BET(I)=ATAN(FXO(I)+ALAM*FX(I))
IF(C(I).LT.0.001.AND.F(I).LT.0.001)BET(I)=0
304 CONTINUE
401 CONTINUE
EWALL=BET(NS) !土压力的作用角度
GFX(1)=GF(1)*COS(BET(1))
GFY(1)=GF(1)*SIN(BET(1))
DO305I=2,NS
PX=X(I)-X(I-1)
DX1=X2(I)-X2(I-1)
G=ALF(I)
G1=ALF(1)
TA=BET(I)
TA1=BET(1)
OP=G-TA
FI=ATAN(F(I)/AB1)
FI1=ATAN(F(1)/AB1)
CE=C(I)/AB1
OQ=FI-G
OQ1=FI1-G1
THO=SIN(TA)*X(I)-COS(TA)*Y(I)-(SIN(TA)*X(1)-COS(TA)*Y(1))
CW=1
IF(I.GT.2)CW=(X2(I-1)-X2(I-2))/(X2(I)-X2(I-1))
IF((OQ+TA).GT.1.5.OR.(OQ1+TA1).GT.1.5.OR.TA.GT,
$ 1.5.OR.TA1.GT.1.5)GOTO317
SEGB=(1./COS(OQ+TA))**2
SEGB1=(1./COS(OQ1+TA1))**2
FMT=FMT+FQUH(I)*(Y(I)-FDIS(I))*DX1
FI0=ATAN(F(I+1)/AB1)
AXS11=ABS((FI-FI0))
OTEM=ABS(OQ+TA)*180./3.14259
IF(OTEM.LT.95.AND.OTEM.GT.85)RETURN
DDQ=(W(I)*SIN(OQ)+CE*COS(FI)/COS(G)-RU(I)*W(I)*SIN(FI)/COS(G)
& -FQUH(I)*COS(OQ))/COS(OQ+TA)!DDQ=P(X)SEC(FI-ALF-BETA)
IF(I.GT.2)GOTO307
PX112=-TAN(OQ1+TA1)*COS(TA1)**2.*(FX(1))
DBS1=TAN(OQ1+TA1)
307 DBS2=TAN(OQ+TA)
EPX112=-DBS2*COS(TA)**2.*(FX(I))
DBS=(DBS1*CW+DBS2)/(1+CW)
DBS1=DBS2
BS=BS+DBS*(BET(I)-BET(I-1))
BSG=BS+DBS*(BET(I)-BET(I-1))*(0.5*DX1)/PX
IF((-BS).GT.30.OR.BSG.GT.30)GOTO317
SX=2.71828**(-BS)!SX=S(X)*COS(FI-ALF+BETA)
SXG=2.71828**BSG
IF(I.GT.2)GOTO309
DTX1=(SIN(TA1)-TAN(G1)*COS(TA1))
309 DTX2=(SIN(TA)-TAN(G)*COS(TA))/SX
DTX=(DTX1*CW+DTX2)/(1+CW)
DTX1=DTX2
TX1=TX+DTX*PX*0.5
TX=TX+DTX*(X(I)-X(I-1))
第11章 程序设计 359
DDP=(-FQUH(I)*SIN(TA)*SIN(FI)+CE*COS(OP)*COS(FI)/
# COS(G)+W(I)*SIN(FI)*COS(TA)-RU(I)*W(I)*SIN(FI)*COS(OP)/
# COS(G))*COS(FI)/(-AB1*COS(OQ+TA)**2)
AAKX1=DDP/DDQ
IF(I.GT.2)GOTO310
DR1X1=(-SIN(FI1)*COS(FI1)*SEGB1/AB1)
310 DR1X2=(-SIN(FI)*COS(FI)*SEGB/AB1)
DR1X=(DR1X1*CW+DR1X2)/(1+CW)
DR1X1=DR1X2
DR1XT=DR1X*TX1
R1XT=R1XT+DR1XT*(BET(I)-BET(I-1))
R1X=R1X+DR1X*(BET(I)-BET(I-1))
DG=DDQ*SX*DX1
DQ=DQ+DG!DQ=DG的积分,DG=P(X)S(X),力的平衡
DDT=DG*(AAKX1-R1X)
DT=DT+DDT!DT=DDT的积分,DDT=P(X)S(X)T(X),力矩平衡
GF(I)=-DQ*SXG
GFX(I)=GF(I)*COS(BET(I))
GFY(I)=GF(I)*SIN(BET(I))
GAVE=(GFY(I)+GFY(I-1))/2,
BETA=(BET(I)+BET(I-1))/2,
DGX=GFX(I)-GFX(I-1)
OX=(GFX(I-1)+GFX(I))/2,
OY=(GFY(I-1)+GFY(I))/2,
IF(WL(I).GE.Y(I))THEN
WP=0,
ELSE
WP=0.5*(RU(I)*W(I))**2/RW
ENDIF
OXX=OX-WP
EFX(I)=OXX
IF(OXX.GT.0.)HSL=Y(I)-YSF(I)
TAA=COS(TA)*COS(TA)*FX(I)
OQS2=(BET(I+1)+BET(I))/2
OQS1=(BET(I-1)+BET(I))/2
DG1=GF(I)*SIN(G-OQS2)-GF(I-1)*SIN(G-OQS1)
BN(I)=W(I)*DX1*COS(G)+DG1-FQUH(I)*DX1*SIN(G)
BT(I)=(BN(I)-RU(I)*W(I)*DX1/COS(G))*TAN(FI)+CE*DX1/COS(G)
EXR1=BN(I)*SIN(G)-BT(I)*COS(G)
EXR=EXR+EXR1
EXT1=FQUH(I)*DX1
EXT=EXT+EXT1
EYR1=-BN(I)*COS(G)-BT(I)*SIN(G)
EYR=EYR+EYR1
EYT1=W(I)*DX1
EYT=EYT+EYT1
XCE=X(1)
YCE=Y(1)
XRE=X(I)-XCE
YRE=Y(I)-YCE
EMT=EMT+EXT1*(FDIS(I)-YCE)-EYT1*XRE!MOMENTOFBODYFORCE
EMR=EMR+EXR1*YRE-EYR1*XRE!MOMENTOFINTERNALFORCE
IF(ABS(DX1).LT.0.001)GOTO317
BN(I)=BN(I)*COS(G)/DX1
BT(I)=BT(I)*COS(G)/DX1
IF(KXYX.EQ.0)GOTO402
GOTO305
360 土质边坡稳定分析?原理
方法
程序
402 I F(I.GT.2)GOTO312
DR2X11=SEGB1*COS(TA1)**2.*(FX(1))
312 DR2X12=SEGB*COS(TA)**2.*(FX(I))
DR2X1=(DR2X11*CW+DR2X12)/(1+CW)
DR2X11=DR2X12
GS=ALF(I-1)
POX=POX-DR2X1*(ALF(I)-GS)
POXT=POXT-DR2X1*(ALF(I)-GS)*TX1
IF(I.GT.2)GOTO313
DPOY1=(COS(FI1)/(COS(G1)*COS(OQ1+TA1)))*COS(TA1)**2.*(FX(1))
313 DPOY2=(COS(FI)/(COS(G)*COS(OQ+TA)*SX))*COS(TA)**2.*(FX(I))
OY1=(DPOY1*CW+DPOY2)/(1+CW)
DPOY1=DPOY2
POY=POY+OY1*PX
DDM=DG*TX
DDMT=DG*(AAKX1*TX-R1XT)
PX112=PX112+DPX112
PT112=PT112+DPT112
DPX112=0
DPT112=0,
PM=POX-PX112
PMT=POXT+POY-PT112
DDEG=DG*PM
DDMG=DG*PMT
DM=DM+DDM
DMT=DMT+DDMT
DEG=DEG+DDEG
DMG=DMG+DDMG
IF(AXS11-0.005)305,305,315
G=ALF(I+1)
315 OQ0=FI0-G
IF((OQ0+TA).GT.1.5)GOTO317
SEGB0=(1./COS(OQ0+TA)**2)
DR2X11=SEGB0*COS(TA)**2.*(FX(I))
DPX112=-TAN(FI0+BET(I+1)-G)*COS(BET(I+1))**2*FX(I+1)
# +TAN(OQ+TA)*COS(TA)**2.*FX(I)
DPT112=DPX112*TX
305 CONTINUE
TEMP1=EPX112+PM
TEMP2=PMT-TX*TAN(OQ+TA)*TAA
IF(IWALL.EQ.2)THEN
GG=-DQ
DDM1=HMW*COS(EWALL)+TX*SX
DDM2=HMW*COS(EWALL)/SX+TX
DM1=DM
DM=DM1+GG/SX*DDM1-FMT
DDDM2=-HMW*COS(EWALL)/SX*TEMP1-TX*TEMP1+TEMP2
DGG=-DEG*(HMW*COS(EWALL)/SX+TX)
DMG1=DMG
DEG1=DEG
DMG=DMG1+GG*(-HMW*COS(EWALL)/SX*TEMP1-HMW/SX*SIN(EWALL)
# *FX(N)*COS(EWALL)*COS(EWALL)-TX*TEMP1+TEMP2)-
# DEG*(HMW*COS(EWALL)/SX+TX)
DB=-DM/DMG
ETA=1.0+DQ/GWALL/SX
EXT=EXT-GWALL*COS(EWALL)*(1-ETA)
EYT=EYT-GWALL*SIN(EWALL)*(1-ETA)
第11章 程序设计 361
EMT=EMT+GWALL*SIN(EWALL)*(X(N)-X(1))*(1-ETA)-
%GWALL*COS(EWALL)*(Y(N)-Y(1)-HMW)*(1-ETA)
GOTO317
ENDIF
IF(IWALL.EQ.1)THEN
DQM=SX*GWALL
DDM1=HMW*COS(EWALL)+TX*SX
DMM=GWALL*DDM1
DTM=-GWALL*SX*R1X
DMTM=-GWALL*SX*R1XT
DEGM=TEMP1*SX*GWALL
DMGM=TEMP2*SX*GWALL
DMGM=DMGM-GWALL*SIN(EWALL)*FX(N)*COS(EWALL)
% *COS(EWALL)*HMW
DQ=DQ+DQM !式(2.31)
DM=DM+DMM !式(2.32)
DT=DT+DTM !DT=DQ对AB1的偏导数
DMT=DMT+DMTM !DMT=DM对AB1的偏导数
DEG=DEG+DEGM !DEG=DQ对ALAM的偏导数
DMG=DMG+DMGM !DMG=DM对ALAM的偏导数
EXT=EXT-GWALL*COS(EWALL)
EYT=EYT-GWALL*SIN(EWALL)
EMT=EMT+GWALL*SIN(EWALL)*(X(N)-X(1))-GWALL*COS(EWALL)
$ *(Y(N)-Y(1)-HMW)
ENDIF
316 CONTINUE
IF(ABS(DM).GT.1.0E+20)GOTO317
IF(ABS(DMT).GT.1.0E+20)GOTO317
DM=DM-FMT !
S1=DEG*DMT-DT*DMG
IF(KXYX.EQ.1)THEN
DF=-DQ/DT
GOTO317
ENDIF
IF(ABS(S1).LE.0.01)GOTO317
DF=(DQ*DMG-DM*DEG)/S1 !式(2.33)
DB=(-DQ*DMT+DM*DT)/S1 !式(2.34)
317 RETURN
END
C
SUBROUTINE MP(N,KXYX,KS,AB1,DQ,DM,BN,FX,FXO,
& C,F,X2,Y2,X,Y,W,RU,FQUH,WL,FDIS,RW,DB,DF)
REAL AB1,DF,DQ,DM,ALAM,DB
INTEGER KS,N,KXYX
REAL X2(80),Y2(80),W(80),X(80),Y(80),C(80),F(80),WL(80)
REAL FQUH(80),FDIS(80),YSF(80),RU(80),RW
REAL GAV,PLA,FA,FB,PLA1
REAL BN(80),FX(80),FXO(80)
COMMON/AS/ASP
!KXYX力矩平衡控制符,=0时同时计算力和力矩的平衡
FA=0,
FB=0.
PLA=0,
IF(KXYX.EQ.0)THEN
GAV=(Y(N)-Y(1))/(X(N)-X(1))
DO277I=2,N
362 土质边坡稳定分析?原理
方法
程序
PLA=PLA+(X(I)-X(I-1))*(FX(I)+FX(I-1))/2,
277 CONTINUE
PLA1=(FA+FB)*(X(N)-X(1))/2,
ALAM=(GAV*(X(N)-X(1))-PLA1)/PLA !ALAM的初值,参见1.3.3节
ELSE
ALAM=ASP
ENDIF
MQ=1
15 CONTINUE
CALL DETE(AB1,RW,DQ,DM,ALAM,KS,BN,FX,FXO,
& C,F,X2,Y2,W,RU,FQUH,YSF,WL,FDIS,N,DB,KXYX,DF)
WRITE(6,*)MQ,DQ,DM,AB1,ALAM
AB1=AB1+DF
IF(KXYX.EQ.0)ALAM=ALAM+DB
IF(MQ.GE.20)GOTO44
MQ=MQ+1
IF(ABS(DF).GT.0.1E-4)GOTO15
IF(ABS(DB).GT.0.1E-3)GOTO15
44 CONTINUE
IF(MQ.GE.20)THEN
WRITE(6,*)'THEITERATIONFAILS'
RETURN
ENDIF
WRITE(6,727)
727 FORMAT(1X,//'THERESULTIS:')
WRITE(6,710)ALAM
710 FORMAT(1X,'LAMBDA=',F12.6)
WRITE(6,711)AB1
711 FORMAT(1X,'FACTOROFSAFETY=',F12.6/)
RETURN
END
3,例题
以11.3.1节的算例说明上节所示源程序的使用计算简图如图11.4所示数据文件分别为s0.dat,s1.dat,s2.dat,s3.dat,s4.dat分别相应于垂直条分法中Morgenstern?Price法
Bishop法瑞典法陆军工程师团法和简化法这些数据文件除控制符OPTION不同外其它均相同(陆军工程师团法在数据最后增加了一行输入平均坝坡数据)文件s0.dat列于后相应迭代过程如下
(1) 数据文件(s0.dat)
0 0
21
-51.522 -50.956 -49.823 -48.690 -47.557
-46.425 -45.292 -44.159 -43.027 -41.894
-40.761 -39.628 -38.496 -37.363 -36.230
-35.097 -33.965 -32.832 -31.699 -30.566
-30.000
-35.000 -33.902 -31.979 -30.491 -29.278
-28.266 -27.411 -26.687 -26.077 -25.567
-25.148 -24.813 -24.559 -24.380 -24.276
-24.244 -24.285 -24.398 -24.585 -24.849
-25.000
1.094647 1.094647,969015,863838,770531
第11章 程序设计 363
,685123,605365,529825,457513,387709
,319849,253490,188260,123833,059923
-.003743 -.067425 -.131382 -.195887 -.261233
-.261233
.000,000,000 38.000,000 76.000 0 0 0
21.415,000,000 38.000 19.500 76.000 0 0 0
57.185,000,000 38.000 19.500 76.000 0 0 0
75.151,000,000 38.000 19.500 76.000 0 0 0
87.756,000,000 38.000 19.500 76.000 0 0 0
96.457,000 5.300 23.000 19.500 76.000 0 0 0
102.081,000 5.300 23.000 19.500 76.000 0 0 0
105.148,000 5.300 23.000 19.500 76.000 0 0 0
106.010,000 7.200 20.000 19.500 76.000 0 0 0
104.914,000 7.200 20.000 19.500 76.000 0 0 0
102.038,000 7.200 20.000 19.500 76.000 0 0 0
97.513,000 7.200 20.000 19.500 76.000 0 0 0
91.434,000 7.200 20.000 19.500 76.000 0 0 0
83.868,000 7.200 20.000 19.500 76.000 0 0 0
74.862,000 7.200 20.000 19.500 76.000 0 0 0
64.439,000 7.200 20.000 19.500 76.000 0 0 0
52.607,000 7.200 20.000 19.500 76.000 0 0 0
39.358,000 7.200 20.000 19.500 76.000 0 0 0
24.663,000 7.200 20.000 19.500 76.000 0 0 0
8.475,000 7.200 20.000 19.500 76.000 0 0 0
.001,000 7.200 20.000,000 76.000 0 0 0
(2) 计算成果(S0.c)
条底倾角α值
1.094647 1.094647 9.690150E-01 8.638380E-01
7.705310E-01 6.851230E-01 6.053650E-01 5.298250E-01
4.575130E-01 3.877090E-01 3.198490E-01 2.534900E-01
1.882600E-01 1.238330E-01 5.992300E-02 -3.743000E-03
-6.742500E-02 -1.313820E-01 -1.958870E-01 -2.612330E-01
-2.612330E-01
水容重= 9.800
DATA FOR ASSUMED SIDE FORCE FUNCTION
************************************
NO,X Y W RU C F BDA WL QH YE
1 -51.522 -35.000,000,000,000,781,000 76.000,000,000
2 -50.956 -33.902 21.415,000,000,781 19.500 76.000,000,000
3 -49.823 -31.979 57.185,000,000,781 19.500 76.000,000,000
4 -48.690 -30.491 75.151,000,000,781 19.500 76.000,000,000
5 -47.557 -29.278 87.756,000,000,781 19.500 76.000,000,000
6 -46.425 -28.266 96.457,000 5.300,424 19.500 76.000,000,000
7 -45.292 -27.411 102.081,000 5.300,424 19.500 76.000,000,000
8 -44.159 -26.687 105.148,000 5.300,424 19.500 76.000,000,000
9 -43.027 -26.077 106.010,000 7.200,364 19.500 76.000,000,000
10 -41.894 -25.567 104.914,000 7.200,364 19.500 76.000,000,000
11 -40.761 -25.148 102.038,000 7.200,364 19.500 76.000,000,000
12 -39.628 -24.813 97.513,000 7.200,364 19.500 76.000,000,000
13 -38.496 -24.559 91.434,000 7.200,364 19.500 76.000,000,000
14 -37.363 -24.380 83.868,000 7.200,364 19.500 76.000,000,000
15 -36.230 -24.276 74.862,000 7.200,364 19.500 76.000,000,000
16 -35.097 -24.244 64.439,000 7.200,364 19.500 76.000,000,000
17 -33.965 -24.285 52.607,000 7.200,364 19.500 76.000,000,000
364 土质边坡稳定分析?原理
方法
程序
18 -32.832 -24.398 39.358,000 7.200,364 19.500 76.000,000,000
19 -31.699 -24.585 24.663,000 7.200,364 19.500 76.000,000,000
20 -30.566 -24.849 8.475,000 7.200,364 19.500 76.000,000,000
21 -30.000 -25.000,001,000 7.200,364,000 76.000,000,000
NOTE,
W=土条的单宽重量
RU=孔隙水压力系数
C=粘聚力
F=摩擦系数
BDA=土条的平均容重
WL=浸润线的Y坐标
SPENCER METHOD
简化法安全系数= 1.325056
迭代过程,
NO,DQ DM AB1 ALAM
1 20.968520 139.277800 1.325056 4.646409E-01
2 9.050338E-01 15.479030 1.344203 3.878648E-01
3 4.864426E-03 4.935455E-02 1.343819 3.800467E-01
THE RESULTS ARE,
LAMBDA=,380018
FACTOR OF SAFETY= 1.343821
上表中DQ和DM为按式(2.31)和式(2.32)计算的不平衡力和力矩AB1和ALAM分别为安全系数F和λ数值
相应s1.dat,s2.dat,s3.dat,d4.dat的安全系数分别为1.384 1.228 1.346和1.350
11,4 边坡稳定分析程序的考核
11,4,1 概述
随着计算机的普及推广几乎所有的岩土工程设计者都使用自己编制的或商用程序对边坡的稳定性进行分析计算鉴于计算成果对工程的经济和安全有着直接影响因此对稳定分析程序的可靠性进行评估和鉴定就成为广受关注的工作
Fredlund (1987)曾对加拿大常用的边坡稳定分析程序作了一次调查并以表格形式介绍各程序的功能和用途但未对这些程序进行考核他的调查发现获得考核的程序对同一算例的计算结果互相差异甚大1980年香港工程发展部岩土工程处 (Geotechnical Control
Office of the Engineering Development Department,Hong Kong) 曾进行过一次测验提出了29
道边坡稳定分析的考核题每道考核题都详细注明了几何断面土的力学强度指标地下水情况表面荷载和计算的滑裂面所有的考核题都是从香港的实际工程中提练出来的42
个单位或个人提交了54份答卷其中75是采用简化Jaubu法进行安全系数计算计算的成果互相接近但是由于使用的是简化Jambu法其成果都和严格方法(MorgesternPrice法
Sarma法)有较大的出入此外这次测验仅局限于对指定的滑裂面计算没有涉及最小安全系数计算Lumsdaine和Tang (1982)曾对这次测验作过详细报道香港工程发展部岩土工第11章 程序设计 365
程处随后建立了一个对各设计部门使用的边坡稳定分析程序进行审核的制度凡是提交该处审查的计算成果必须使用经审定的程序程序在报审时必须提交五道考核题的答卷包括全部的输入输出数据和使用手册
自1982年开始水利电力部规划设计院在潘家铮总工程师主持下曾对水电系统常用的土石坝边坡稳定分析和其它数值分析程序进行开发和鉴定工作原水利电力部昆明设计院根据天生桥水电站工程的实际情况设计了考核题
在国内外对边坡稳定分析程序进行考核的工作中澳大利亚计算机应用协会的一项调查工作细致并已正式记载于文献(Donald & Giam,1992)在下节中将着重介绍这一考核工作的成果
11,4,2 澳大利亚ACADS边坡稳定分析程序调查
1987年澳大利亚计算机应用协会(ACADS)对澳大利亚所使用的边坡稳定分析程序进行一次调查他们委托Monash大学的B,Donald教授和P,Giam博士主持此项工作此次调查设计的考核题五道总计十个小题向120个单位的发出了邀请有28个单位发回了计算答案东道主Donald教授还邀请了国际上在边坡稳定分析程序方面做过较多工作的研究者提供裁判程序(Referee Program)答案本书作者有幸与加拿大Saskatchwen大学Fredlund
教授和以色列工学院的Baker教授一起被邀请担任此项工作本研究成果已于1989年4月
19日在澳大利亚岩土工程协会维克多利亚分会的月会上公布,后正式发表(Donald & Giam,
1992)由于这次调查工作规模较大所获得的成果比较可靠本节详细介绍有关的考核题和答案这些内容可作为各单位和个人对自编程序考核的参考资料
1,考核题设计
(1) 考核题1本题包括四小题图11.5图11.6及表11.1和表11.2分别给出本题剖面的几何尺寸及物性指标对所有四种情况要求确定临界滑裂面(圆弧或非圆弧)以及最小安全系数
考核题1(a)为一均质边坡材料特性如表11.1所示
考核题1(b)为同样几何形状的均质边坡但材料特性如表11.2所示要求选择一个合适的拉力缝深度并假定缝内充水
考核题1(c)为非均质边坡材料特性如表11.3所示相信读者对考核题1(c)已很熟悉了因为在第2第4第10章中已多次使用
考核题1(d)情况与考核题1(c)情况相同只是土体作用一水平地震加速度系数为0.15g
366 土质边坡稳定分析?原理
方法
程序
图 11,5 ACADS考核题EX1(a)和EX1(b)
图 11,6 ACADS考核题EX1(c)和EX1(d)
表 11,1 考核题1(a)材料性质
c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
备注
3.0 19.6 20.0 1.0E4 0.25 0.65 1
#
土
表 11,2 考核题1 (b)材料性质
c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
备注
32.0 10.0 20.0 1.0E4 0.25 0.65 1
#
土
注 这一问题需设一合适的拉力缝
表 11,3 考核题1(c d)材料性质
土号 c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
1
#
土 0.0 38.0 19.5 1.0E4 0.25 0.65
2
#
土 5.3 23.0 19.5 1.0E4 0.25 0.65
3
#
土 7.2 20.0 19.5 1.0E4 0.25 0.65
(2) 考核题2本题选自Talbingo坝实例图11.7及表11.4~表11.6分别给出其剖面的第11章 程序设计 367
几何尺寸及物性指标材料1为堆石填料材料2为心墙料
考核题2(a)要求确定上游坡临界滑裂面及相应最小安全系数
考核题2(b)要求对指定圆弧计算其稳定安全系数圆弧圆心坐标和半径值参见表11.6
(3) 考核题3本题的边坡剖面几何特征和材料特性示于图11.8和表11.7地下水位假定位于软弱夹层底部地下水位以上孔压为零考核题3(a)要求确定临界滑裂面(圆弧或非圆弧)和相应最小安全系数考核题3(b)要求计算指定滑裂面(其坐标如表11.8所示)的安全系数
考核题3的数据如表11.7和表11.8
表 11,4 考核题2几何剖面数据
点号 X (m) Y (m) 点号 X (m) Y (m)
1 0 0 14 648 0
2 315.5 162 15 168.1 0
3 319.5 162 16 302.2 130.6
4 321.6 162 17 200.7 0
5 327.6 162 18 311.9 130.6
6 386.9 130.6 19 307.1 0
7 394.1 130.6 20 331.3 130.6
8 453.4 97.9 21 328.8 146.1
9 460.6 97.9 22 310.7 0
10 515 65.3 23 333.7 130.6
11 521.1 65.3 24 331.3 146.1
12 577.9 31.4 25 372.4 0
13 585.1 31.4 26 347 130.6
注 点号参考图11.7
表 11,5 考核题2材料性质
c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
堆石 0.0 45 20.4 5.0E4 0.30 0.65
过渡层 0.0 45 20.4 5.0E4 0.30 0.65
反滤层 0.0 45 20.4 4.0E4 0.30 0.65
心墙 85 23 18.1 1.0E4 0.44 0.65
注 堆石过渡层和反滤层几乎具有同样的材料性质
368 土质边坡稳定分析?原理
方法
程序
图 11,7 ACADS考核题EX2
表 11,6 考核题2(b)圆弧滑裂面信息
X (m) Y(m) 半径(m)
100.3 290.0 278.8
注 x y轴的原点位于上游堆石边坡和坝基的交点
表 11,7 考核题3材料性质
c(kN/m
2
) φ(°) γ(kN/m
3
) E(kN/m
2
) ν K
0
1
#
土 28.5 20.0 18.84 6.0E4 0.25 0.65
2
#
土 0.0 10.0 18.84 2.0E3 0.25 0.65
注 假设水面线与基础中的软弱夹层相重合(夹层和夹层上的孔压=0)不考虑水面线以上负的孔隙水压力
表 11,8 考核题3(b)滑裂面控制点坐标
X(m) Y(m)
41.85 27.75
44.00 26.50
63.50 27.00
73.31 40.00
图 11,8 ACADS考核题EX3
(4) 考核题4本题为一坡面上作用表面荷载坡内有浸润线的例子图11.9及表11.9
表11.10和表11.11分别给出本题的剖面几何尺寸物性指标浸润线及表面荷载要求确第11章 程序设计 369
定临界滑裂面和相应最小安全系数
表 11,9 考核题4材料性质
c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
1
#
土 28.5 20.0 18.84 6.0E4 0.25 0.65
2
#
土 0.0 10.0 18.84 2.0E3 0.25 0.65
表 11,10 考核题4中的外荷载
X(m) Y(m)
正应力
(kN/m
2
)
23.00 27.75 20.00
43.00 27.75 20.00
70.00 40.00 20.00
80.00 40.00 40.00
表 11,11 考核题4的浸润线数据
点号 X(m) Y(m)
1 20.0 27.75
2 43.0 27.75
3 49.0 29.8
4 60.0 34.0
5 66.0 35.8
6 74.0 37.6
7 80.0 38.4
9 84.0 38.4
注 点号参考图11.9
图 11,9 ACADS考核题EX4
(5) 考核题5考核题5如图11.10图11.17材料性质如表11.12要求计算最小安全系数
表 11,12 考核题5材料性质
c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
K
x
=K
y
(m/s)
11.0 28.0 20.00 1.5E4 0.25 0.65 5.0E-08
370 土质边坡稳定分析?原理
方法
程序
图 11,10 ACADS考核题EX5
2,计算成果汇总和分析
如前所述组织者向澳大利亚各咨询公司发出了问卷调整对各题目返回的答案数目在
7~33个之间表11.13总结了收到的答案
本次调查邀请了国外三个程序提供裁判答案Donald教授也使用了自编的GWEDGEM
和EMU程序分析了这些考核题综合各咨询公司提供的答案裁判答案和Donald自编程序的答案在表11.13汇总本次考核最终提出了推荐的裁判答案
现对各题的解答情况作一简要回顾
表 11,13 返回的问卷答案统计成果
安全系数
考核题号 分析方法
均值 标准差 min.F max.F
提交答案总数
(ALL) 0.991 0.031 0.940 1.080 32
BISHOP 0.993 0.015 0.960 1.030 18
1(a)
(1.00)
JANBU 0.978 0.041 0.940 1.043 7
(ALL) 1.381 0.067 1.230 1.520 30
BISHOP 1.406 0.039 1.350 1.500 16
1(c)
(1.39)
JANBU 1.325 0.057 1.264 1.400 6
(ALL) 0.972 0.054 0.900 1.100 14 1(d)
(1.00)
(ALL) 2.003 0.101 1.830 2.201 23 2(a)
(1.90) BISHOP 2.042 0.118 1.830 2.201 13
(ALL) 2.237 0.112 1.948 2.390 23 2(b)
(2.29) BISHOP 2.204 0.013 2.180 2.220 11
(ALL) 1.369 0.168 1.191 1.820 25
BISHOP 1.559 0.146 1.400 1.820 7
JANBU 1.307 0.138 1.191 1.564 8
3(a)
(1.26?)
NON-CIRCULAR 1.295 0.108 1.191 1.564 18
(ALL) 1.292 0.064 1.197 1.450 29 3(b)
(1.34) JANBU 1.252 0.060 1.197 1.333 10
(ALL) 0.897 0.254 0.410 1.360 25
BISHOP 1.151 0.175 0.850 1.360 6
JANBU 0.984 0.198 0.712 1.287 7
4
(0.78)
NON-CIRCULAR 0.817 0.223 0.410 1.287 19
5 (ALL) 1.462 0.150 0.910 1.600 20
第11章 程序设计 371
(1.35) BISHOP 1.431 0.192 0.910 1.600 10
(1) 考核题1本题是断面相对比较简单的土质边坡包括了四个小题主要考核用优化方法来确定的临界滑裂面和最小安全系数裁判答案作者均使用Bishop简化法从表11.14计算成果看对考核题1(a) 1(c)和1(d)
裁判答案各家给出的结果完全一致图11.11给出了对考核题1采用STAB
求得的临界滑裂面和相应的安全系数图11.12给出了对考核题1各家程序求得的临界滑裂面和相应的安全系数可以看出不同的程序得到的临界滑裂面都落在一个窄的区域里相应的安全系数也基本一致这说明对于土质边坡圆弧滑裂面目前水平已可以达到不同的程序均可获得一致成果的程度自动搜索临界滑裂面的功能也已经达到实用阶段
图 11,11 总结返回的考核题EX1计算成果的临界滑裂面
(a) EX1(a) (b) EX1(b) (c) EX1(c) (d) EX1(d)
表 11,14 裁判程序答案汇总(考核题1)
考核题 1(a) 1(b) 1(c) 1(d)
推荐的裁判答案Donald 1.00 1.65~1.70 1.39 1.00
SSA(Baker) 1.00 1.50 1.39
STAB(Chen) 0.991 1.618 1.385 1.007
GWEDGEM 1.00 1.65~1.70 1.39 1.00
EMU 1.00 1.65~1.70 1.39 1.00
Fredlund 0.99 1.619~1.588 1.406 1.015
372 土质边坡稳定分析?原理
方法
程序
图 11,12 各家程序对考核题EX1的计算成果(临界滑裂面位置和最小安全系数)
(a) EX1(c); (b) EX1(d)
图例SR? Sarma法B Bishop法JG Janbu法GW GWEDGEM法J简化扬布法MP Morgensten?Price法
考核题1(b)需要用户自行设计一个拉力缝这样由于设有的拉力缝高度不同安全系数有一定的出入但各家的差距也不大H
t
均在0.5~3.5m范围内图11.11(b)示出了笔者选用
H
t
=2.0 m时充水拉力缝的位置及临界滑裂面和相应的安全系数
因此考核题1可以作为一个经典的考核题考核自编的边坡稳定分析程序的准确性
(2) 考核题2本题分为两部分第一部分是求上游坝坡的安全系数由于坝壳料为无粘聚力的砾质材料故如使用线性强度指标将获得一个非常浅的临界滑裂面而其安全系数就是摩擦系数和坝的坡度之比即F
=
tanφ / tanα
=
1.947各家得到的结论是相同的考核第11章 程序设计 373
题还要求在厚心墙内寻找一个局部极值从表11.15中可见GWEDGEM和EMU给出的局部极值为1.90而Fredlund和陈祖煜的答案是2.14~2.15图11.13(a)为STAB程序获得考核题EX2的临界滑裂面和最小安全系数图11.14(a)列出了考核题2(a)各家得到的临界滑裂面和相应安全系数
第二部分是计算指定的圆弧滑裂面的安全系数,如图11.13(b)各家计算成果偏离不大
Baker和Donald的成果均为2.29陈祖煜和Fredlund的结果分别为2.220和2.206
表 11,15 裁判程序答案汇总(考核题2~考核题5)
2(a) 2(b) 3(a) 3(b) 4 5
推荐的裁判答案 (1.95)/1.90 2.29 1.26 1.34 0.78 1.53/1.53/1.59
SSA+(baker) 1.95/ 2.35 1.24 1.37 0.65 1.53
STAB+(Chen) 1.947~2.135 2.224 1.242 1.265 0.671 1.503
GWEDGEM 1.95/1.90 2.29 1.27 1.34 0.78 1.53/1.53/1.59
EMU 1.95/1.90 2.29 1.27 1.34 0.78 1.53/1.53/1.59
Fredland 2.138 2.206 1.261 0.630 1.460
图 11,13 STAB程序获得考核题EX2的临界滑裂面和最小安全系数
(a) 考核题 EX2(a) (b) 考核题 EX2(b)
(3) 考核题3本题为典型的沿软弱夹层进行非圆弧滑裂面的计算考核题考核题3(a)
要求计算临界滑裂面考核题3(b)要求计算指定滑裂面的安全系数STAB采用
Morgenstern?Price法图11.15所示STAB程序获得考核题EX3的临界滑裂面和最小安全系数对指定滑裂面算得安全系数为1.265[考核题3(b)]其值与Fredlund的解答1.261接近但Baker和Donald等给出了稍大的结果(1.34~1.37)最终各家搜索获得的最小安全系数比较接近在1.24~1.27之间图11.16为各家给出了本次调查考核题3(a)的临界滑裂面和相应的安全系数
374 土质边坡稳定分析?原理
方法
程序
图 11,14 各家程序对考核题EX2的计算成果(临界滑裂面位置和最小安全系数)
图 11,15 STAB程序对考核题EX3计算的滑裂面和安全系数
(a)考核题 EX3(a); (b) 考核题 EX3(b)
图 11,16 各家对考核题EX3计算的临界滑裂面和最小安全系数
(4) 考核题4本例在设计方面存在一定的问题软弱夹层的上部没有标明是否与边坡第11章 程序设计 375
面相交由于夹层的指标明显低于土层导致在搜索过程中临界滑裂面不断上移实际上直正的临界滑裂面可能要越出规定的图框临界滑裂面除趾部外全部通过软弱夹层另外
Bishop法在此例中不适用因控制滑裂面通过软弱夹层不可能是圆弧形的因此表11.15
提供的答案仅供参考
(5) 考核题5本文为一个地基快速开挖的例题考核要求采用三种不同的孔隙水压力处理方法把所给的地下水面作为测压管面把图11.17所给的流网转化为孔压网来计算根据渗流分析计算更精确的孔隙压力除外各家根据图11.17所示流网计算孔隙水压再用孔压计算安全系数表11.15分别示出这三种处理方法Baker和Donald的计算成果
STAB仅作了第一种情况如图11.18示各家提供的安全系数在1.46~1.53之间
图 11,17 考核题5渗流场
图 11,18 STAB程序对考核题EX5计算的临界滑裂面和最小安全系数
ACADS组织的这次考核的特点是详细地规定了计算原始数据和要求使不同的程序获得的成果具有可比性除第4题外其它考核题的设计都比较严密因此有条件使不同的程序算出基本相同的结果但由组织者提供的2(a)和3(b)的推荐的裁判答案与Fredlund
和陈祖煜提供的答案相差较大其权威性有待进一步论证综合参加的各程序提供的安全系数统计参数成果表11.13可见裁判答案互相比较接近而各单位自编程序提供的答案则比较分散早期Fredlund也曾经发现过这一问题这一事实说明在进行稳定分析时最好还是使用专门开发的经过考核的比较成熟的程序
边坡稳定分析是岩土工程设计和施工中经常要进行的工作在计算机迅速普及的今日
376 土质边坡稳定分析?原理
方法
程序
手算方式显然要被淘汰但确认各家自编程序的准确性确实成为了进入计算机时代面临的一个重要的问题已提到了议事日程这次程序考核工作做了一次有益的尝试本项工作为各单位自编的稳定分析程序的考核提供了十分有用的资料
11,5 本章附录
11,5,1 本章数据文件表11.6
表 11,16 本 章 数 据 文 件
有关章节 系列号 数据文件名 内容
11-01-01 CIRCLE.FOR 构筑圆弧滑裂面源程序
11-01-02 CIRCLE.DAT 构筑圆弧滑裂面例题
11.3.2
11-01-03 CIRCLE.C 构筑圆弧滑裂面例题计算结果
11-02-01 SPLINE.FOR 构筑任意形状滑裂面源程序
11-02-02 SPLINE-1.DAT 构筑光滑形状滑裂面例题
11-02-03 SPLINE-1.C 构筑光滑形状滑裂面例题计算结果
11-02-04 SPLINE-2.DAT 构筑组合形状滑裂面例题
11.3.3
11-02-05 SPLINE-2.C 构筑组合形状滑裂面例题计算结果
11-03-01 PROFILE.FOR 计算土条重量等特性数据的源程序
11-03-02 PROFILE.DAT 计算土条重量等特性数据的例题
11.3.4
11-03-03 PROFILE.C 计算土条重量等数据的例题计算结果
11-04-01 S.FOR 计算安全系数源程序
11-04-02 S0.DAT 计算通用条分法安全系数例题
11-04-03 S0.C 计算通用条分法安全系数例题计算结果
11-04-04 S1.DAT 计算Bishop法安全系数例题
11-04-05 S1.C 计算Bishop法安全系数例题计算结果
11-04-06 S2.DAT 计算瑞典法安全系数例题
11-04-07 S2.C 计算瑞典法法安全系数例题计算结果
11-04-08 S3.DAT 计算工程师团法安全系数例题
11-04-09 S3.C 计算工程师团法安全系数例题计算结果
11-04-10 S4.DAT 计算简化法安全系数例题
11.3.5
11-04-11 S4.C 计算简化法法安全系数例题计算结果
11.4.2 11-05-01 exp1a.dat ACADS考核题1(a)
11-05-02 exp1b.dat ACADS考核题1(b)
11-05-03 exp1c.dat ACADS考核题1(c)
11-05-04 exp1d.dat ACADS考核题1(d)
11-06-01 Exp2a.dat ACADS考核题2整体极值
11-06-02 Exp2b.dat ACADS考核题2心墙局部极值
11-06-03 Exp2c.dat ACADS考核题2指定滑裂面
11-07-01 Exp3a.dat ACADS考核题3整体极值
11-07-02 Exp3b.dat ACADS考核题3心墙局部极值
11-08-01 Exp4.dat ACADS考核题4
11-09-01 Exp5.dat ACADS考核题5
第11章 程序设计 377
参考文献
1 刘德贵费景高于泳江李广元,Fortran算法汇编第一分册,北京国防工业出版社
1980
2 Donald,I,B,And Giam,P,The ACADS slope stability programs review,Proc,6
Th
International Symposium
on Landslides,1992,Vol,3,1665-1670
3 Lumsdaine,R,W,and Tang,K,Y,A comparison of slope stability calculations1982,Proc,7
th
South East Asia
Geotechnical Conference,Hong Kong,31-38
4 Fredlund,D,G,Usage,requirements and features of slope stability computer software,Canadian Geotechnical
Journal,1987,15:83-95
5 ACADS Publication No,U255,Soil slope stability programs review,April,1989
11,1 前言
从本书介绍的内容可以发现随着近代数值分析技术的进展边坡稳定分析的极限平衡分析方法已经发展到十分成熟的程度但是随之而出现的问题是这一分析步骤必须通过计算机程序才能得以实现而现状是尽管计算机的硬件已经获得了飞速的发展尽管土力学的理论和实现这些理论的工具和手段均已十分先进但岩土工程许多应用软件仍然处于起步状态远不能满足工程实际的要求且大部分岩土工程师对边坡稳定分析计算的实际应用水平还不很高使用近代科学技术的成果来提高岩土工程设计的效率和水平是工程师的强烈愿望也是保证建筑物的安全和减少投资非常重要手段当前岩土工程软件的落后状况是有多方面因素造成的
首先开发岩土工程软件需要投入大量的人力和物力研究一个新的数值分析方法并编制相应的科研程序是一回事将其变为一个能够在工程设计中广泛应用的软件是另一回事其源代码的数量起码增加十倍在Windows的界面上编制一个直观方便的具有处理图形对话和纠错功能的软件不仅是源代码数量上的增加而且意味着对编程人员素质的新要求而这一要求是由所支付的高额报酬体现在程序价格中的
其次岩土工程软件的市场极为有限对于这种专业性的软件单靠市场回报来实现其自身的发展是不现实的而这一问题并未得到共识众所周知设计人员在图纸上多画一根预应力锚索多画一根抗滑桩则是十万甚至百万元的投资但对于开发这些有关软件的研究经费(尽管无法与施工费用相比)却难以得到批准
再次我们的软件管理还远未达到规范化市场化的水平盗版现象普遍存在因此营造一个保护软件的法律和道德约束环境是软件作为一个产业兴旺发达的最基本的条件
工程应用软件开发的另一个重要问题是确认其可靠性和准确性工程软件应用的对象是关系到人民生命财产安全的建筑物因此绝非儿戏对程序进行鉴定和考核自然是一个必须屡行的手续但笔者认为对一个包含成千上万行指命的又需要面临复杂多变的实际情况的工程应用软件单靠鉴定和考核是不可能真正发现其全部问题的这样的程序只有在长期的大量的实际应用过程中才能变得成熟成为可以信赖的设计工具
作者希望通过第11章和第12章的内容推动边坡稳定计算机程序的普及和推广本章将介绍编制一个边坡稳定分析程序的核心步骤和源程序将这几个子程序串联起来就是一个计算任意几何形状边坡包含第2章和第3章介绍的各种分析方法的程序本章还将讨论与程序推广使用相关的问题在第12章作者提供STAB程序的一个简化版使用这一版本读者可以进行本书介绍的部分稳定分析工作需要说明的是本书所列的所有源程序
336 土质边坡稳定分析?原理
方法
程序
只供科研和教学使用作者不对其承担法律责任他人也不得将其应用于商业目的
11,2 几何图形的识别和分析
11,2,1 边坡剖面图形的数字信息
通常坝体并非均一而是由若干个物理力学指标不同的土质区域组成图11.1所示坝体由6个土质不同的区域组成控制点17个边界线22条图中滑弧上阴影所示的土条中心线和3条边界线相交将该土条分成了三段其高度分别为H
1
H
2
和H
3
并分别位于
I区V区和VI区所有这些在图解法中依靠人的直观操作是很容易判断和量测得的编制程序的关键是如何让计算机代替人判断土条的中心线究竟与哪几个线条相交各段高度分别为多少我们的办法是让计算机对所有的边界线循环一遍逐个判断每一条线是否与土条中心线相交边界线用所压土层号的数组来表示输入这些数据如交上则记下该交点的X Y坐标值及该线段所压的土层号这样根据相邻两交点的Y坐标值的差值即可算得H
1
H
2
和H
3
根据土层编号即可查得相应的物理力学指标随即可算得土条重量地震力孔隙水压力等这一对几何图形进行处理的方法充分利用了计算机可进行大量重复的判断忌讳过多的逻辑判断的特点为编制一个效率高鲁棒性 (Robust) 强的边坡稳定分析程序创造了条件
图 11,1 由6个土质不同的区域组成坝体
11,2,2 构造圆弧滑裂面的步骤
设某一坝坡坡面为具有n个折点的折线如图11.2所示(其中起点1和尾点n在极远处)
建立坐标系统X轴水平指向滑坡方向为正Y坐标轴垂直向下为正某一半径为R
0
的滑弧圆心O′ 至坝坡各折点距离分别为R
1
,R
2
,…,R
n
该圆弧自左至右在连接i i+1两点的线段进入坝体上交点在连接j j+1线段逸出坝体下交点
不难证明某一半径为R
0
的圆弧与线段A
1
A
2
有一个且只有一个交点的必要充分条件是
R
2
R
0
R
1
其中R
1
R
2
为圆心至线段端点的距离如果R
0
同时大于R
1
R
2
则圆弧与A
1
A
2
没有交点如果R
0
同时小于R
1
R
2
则圆弧与A
1
A
2
不是没有交点就是有两个交点,见图第11章 程序设计 337
11.3
则根据上述原理可知i j应满足下面条件
图 11,3 具有n个折点的边坡线与圆弧滑裂面相交的关系
图 11,2 圆弧滑裂面与边坡相交的关系
(11.1)
ii
RRR <≤
+ 01
(11.2)
10 +
<≤
jj
RRR
找到了上下交点所在的线段号i j之后就可以根据该线段两端点的坐标值和圆弧的圆心坐标值x
0
y
0
和半径值R算得其交点的坐标值对于点m (m = i,j)计算公式为
(11.3)
22
0
2
0
)()( Ryyxx =?+?
(11.4)
sxxyy
mm
= )/()(
(11.5)
)/()(
11 mmmm
xxyys=
++
联立以上方程组求解得
])()1([
)1(
1
2
00
2
02
TsxRssTx
s
x+±+
+
=
(11.6)
(11.7)
0
yysxT
mm
+?=
此联立方程组有两个根计算下交点时取+号计算上交点时取?号
11,2,3 构造任意形状滑裂面的步骤
如果要将若干点用光滑的曲线相连可以采用样条函数予以实现在STAB程序中采用了第二种边界条件的三次样条函数
第二种边界条件三次样条函数曲线两端点的二阶导数值为已知值即)(
112
xyy ′′=
)(
2 nn
xyy ′′=在构造滑裂面时取0)(
112
=′′= xyy 0)(
2
=′′=
nn
xyy由此确定的样条称为
338 土质边坡稳定分析?原理
方法
程序
自然样条根据此边界条件计算n个节点上的一阶导数值m nj
j
,,2,1,L=相应的计算公式如下
(11.8) 5.0
1
=a
)(2
)()(
3
12
12
1
xx
xyxy
b
×= (11.9)
(11.10)
1,,2,1
1
=?=
+
njxxh
jjj
L
)(
11 jjjj
hhh +=
α
(11.11)
])())(1[(3
111 jjjjjjjjj
hyyhyy?+=
+
ααβ (11.12)
])1(2[
1?
+?=
jjjj
aa αα
(11.13)
1,,3,2])1(2[])1([
11
=?+=
njabb
j
j
jjjj
Lααβ (11.14)
)2/(])(3[
1111
+=
nnnnnn
abhyym (11.15)
(11.16)
1,,2,1
1
L=+=
+
nnjbmam
jjjj
于是对于,利用数值导数 (j )
1+
<≤
jj
xxx
j
m n,,2,1L=计算插值节点处函数值的公式如下
1
3
3
2
2
3
1
3
2
1
2
1
3
3
2
2
3
1
3
2
1
2
)(
2
)(
1
)(
1
)(
1
)(
2
)(
3
)(
2
)(
3
)(
+++
+++
+
+
=
jj
j
j
j
fjj
j
j
j
f
jj
j
j
j
jj
j
j
j
mxx
h
xx
h
hmxx
h
xx
h
h
yxx
h
xx
h
yxx
h
xx
h
xy
(11.17)
如果插值区间被等距划分为l段),(
1+jj
xx则位于第s (0<s<l) 段的被插值节点
jj
hlsxx )(+=处的函数值的计算公式如下
(11.18)
)]([)(
1234
κκκκ
sjsjsj
hhhxy +++=
其中
lshxxh
jjsj
=?= )(
(11.19)
(11.20)
jjjjj
hmmxyxy )()]()([2
111 ++
++?=κ
(11.21)
jjjjj
hmmxyxy )2()]()([3
112 ++
+=κ
(11.22)
jj
hm?=
3
κ
(11.23))(
4 j
xy=κ
第11章 程序设计 339
11,3 边坡稳定分析的源程序
11,3,1 概述
本章介绍边坡稳定分析程序STAB的核心步骤和程序程序 CIRCLE.FOR具体实施
11.2.2节的计算步骤对给定的边坡剖面信息构造一个圆弧滑裂面程序SPLINE.FOR具体实施11.2.3节的计算步骤对给定的滑面控制点信息构造一个任意形状滑裂面程序
PROFILE.FOR具体实施11.2.1节的计算步骤 计算土条重量和物理力学参数程序S.FOR
计算安全系数包含了第2章和第3章介绍的各种分析方法的源程序即极限平衡法的通用条分法瑞典法毕肖普法陆军工程师团法和简化法将上述几个子程序串联起来就是一个计算复杂几何形状边坡稳定的安全系数的程序
在介绍这些程序时将配以图11.4所示的一个例子这一例子是将在11.4节中介绍的澳大利亚ACADS程序考核的EX1(c)该例已在第2章和第4章中多次使用读者在移植这些程序时可通过这些考核题进行考核
图 11,4 滑裂面和条分
11,3,2 构造圆弧滑裂面的程序 CIRCLE.FOR
1,说明
本子程序对给定的边坡剖面信息构造一个圆弧滑裂面(图11.4)通过子程序SEAR找到圆弧滑裂面与边坡的上下交点所在线段子程序DIVI计算上下交点的坐标并划分土条
2,源程序
REAL YTENSION,DS,CX,CY,XN(99),YN(99),X(80),Y(80),ALF(80)
INTEGER IN1,LSL,N,LNUM(80),IC(80,3),IWR5,L,NN,NLOW,NUPP
OPEN(5,FILE='',STATUS='UNKNOWN')
OPEN(6,FILE='',STATUS='UNKNOWN')
NLOW=0
NUPP=0
READ(5,*)IWR5,YTENSION,DS,CX,CY
C IWR5=0,1分别相应顶部有无拉力缝YTENSION=拉力缝底Y坐标
C DS,CX,CY=滑弧深度圆心坐标
READ(5,*)N!土条总数
READ(5,*)IN1!线段总数
DO 302 I=1,IN1
READ(5,*)(IC(I,J),J=1,3)!线段信息
302 CONTINUE
340 土质边坡稳定分析?原理
方法
程序
READ(5,*)LSL!外边坡线总数
READ(5,*)(LNUM(I),I=1,LSL)!外边坡线编号
READ(5,*)NN!点总数
READ(5,*)(L,XN(I),YN(I),I=1,NN)!点坐标
WRITE(6,703)
703 FORMAT(T25,'******************'//)
WRITE(6,704)
704 FORMAT(10X,'THEABSCISSAVALUESOFTHEUPPERANDLOWER')
CALL SEAR(LNUM,IC,DS,CX,CY,XN,YN,LSL,NLOW,NUPP)!寻找上下交点所在线段
CALL DIVI(N,IC,IWR5,DS,CX,CY,XN,YN确定上下交点坐标条分
$,YTENSION,NLOW,NUPP,X,Y,ALF)
WRITE(6,711)
711 FORMAT(T5,'NO.',T17,'X',T32,'Y',T47,'ALF')
DO 309 I=1,N
WRITE(6,710)I,X(I),Y(I),ALF(I)
710 FORMAT(1X,I5,3F15.6)
309 CONTINUE
12 CALLCLOSEFL
END
SUBROUTINE CLOSEFL
CLOSE(5)
CLOSE(6)
RETURN
END
SUBROUTINE DIVI(N,IC,IWR5,DS,CX,CY,XN,YN
$,YTENSION,NLOW,NUPP,X,Y,ALF)
INTEGERN S,N,K,IC(80,3),K1,KTEMP,IWR5,J,J1,JTEMP,NLOW,NUPP
REAL RO,DS,CY,XN(80),YN(80),SK,T,XL,CX,XU1,YTENSION
REAL YU1,SJ,XU,ROO,YU,X2(80),Y2(80),X(80),Y(80),ALF(80)
S(E1,F1,E2,F2)=(F2-F1)/(E2-E1)
NS=N-1
RO=DS-CY
K=IC(NLOW,1)
K1=IC(NLOW,2)
IF(XN(K).GT.XN(K1))GOTO190
KTEMP=K
K=K1
K1=KTEMP
190 IF(ABS(XN(K)-XN(K1)).LT.0.00005)GOTO23
SK=S(XN(K),YN(K),XN(K1),YN(K1))
T=SK*XN(K)-YN(K)+CY
XL=(CX+SK*T+SQRT((1+SK**2)*RO**2-(SK*CX-T)**2))/(1+SK**2)
!按式(11.6)计算下交点的横坐标
GOTO28
23 XL=XN(K)
28 CONTINUE
IF(IWR5.NE.0)THEN
XU1=CX-SQRT(RO*RO-(YTENSION-CY)**2)
YU1=CY+SQRT(RO**2-ABS(XU1-CX)**2)
ENDIF
J=IC(NUPP,1)
J1=IC(NUPP,2)
IF(XN(J).GT.XN(J1))GOTO191
JTEMP=J
J=J1
J1=JTEMP
第11章 程序设计 341
191 F(ABS(XN(J)-XN(J1)).LT.0.00005)GOTO82
SJ=S(XN(J),YN(J),XN(J1),YN(J1))
T=SJ*XN(J)-YN(J)+CY
XU=(CX+SJ*T-SQRT((1+SJ**2)*RO**2-(SJ*CX-T)**2))/(1+SJ**2)
!XU上交点的横坐标
ROO=RO**2-ABS(XU-CX)**2
IF(ROO.LT.0)RETURN!无交点,需修改RO,重新计算
YU=CY+SQRT(RO**2-ABS(XU-CX)**2)
IF(CY.GT.YN(J).AND.CY.GT.YN(J1))RETURN!凸弧,否定,重算
IF(IWR5.NE.0.AND.YU.LT.YU1)XU=XU1
GOTO20
82 XU=XN(J)
20 X2(1)=XU
DO541I=1,NS!划分土条,求各点坐标
X2(I)=(I-1)*(XL-XU)/(NS-1)+XU!土条宽度B=(XL-XU)/(NS-1)
Y2(I)=CY+SQRT(RO**2-ABS(X2(I)-CX)**2)
541 CONTINUE
IF(ABS(Y2(NS)-Y2(1)).LT.0.01)THEN
WRITE(6,*)'THE SLIP SURFACE INTERCEPTS A HORIZONTAL SURFACE'
RETURN
ENDIF
N=NS+1
X2(N)=X2(NS)
Y2(N)=Y2(NS)
C 计算土条中点坐标
DO544I=2,NS
X(I)=(X2(I)+X2(I-1))/2,
Y(I)=(Y2(I)+Y2(I-1))/2,
ALF(I)=ATAN((Y2(I)-Y2(I-1))/(X2(I)-X2(I-1)))
544 CONTINUE
ALF(1)=ALF(2)
ALF(N)=ALF(NS)
X(N)=X2(N)
Y(N)=Y2(N)
X(1)=X2(1)
Y(1)=Y2(1)
RETURN
END
SUBROUTINE SEAR(LNUM,IC,DS,CX,CY,XN,YN,LSL,NLOW,NUPP)
REAL R,DS,CY,R2,CX,XN(80),YN(80),R1,D,XXL,SLO
INTEGER I,I1,K,LNUM(80),J1,J2,IC(80,3),IN,NLOW,NUPP,LSL
I=1 !NLOW=下交点所在直线段的线段号
I1=0 !NUPP=上交点所在直线段的线段号
12K=LNUM(I)
J1=IC(K,1)
J2=IC(K,2)
R=DS-CY!R=圆弧半径
R2=SQRT(ABS(CX-XN(J2))**2+ABS(CY-YN(J2))**2)
R1=SQRT(ABS(CX-XN(J1))**2+ABS(CY-YN(J1))**2)
!R1,R2为圆心至该线段两端点的距离
IF(R.EQ.R2.AND.I1.NE.0)GOTO11
IF(R.GT.R1.AND.R.GT.R2)GOTO11 !无交点
IF(R.LT.R1.AND.R.LT.R2)GOTO19 !有两个或零个交点
GOTO20
19 CALL DD(XN(J1),YN(J1),XN(J2),YN(J2),CX,CY,D,IN)
XXL=XN(J2)-XN(J1)
342 土质边坡稳定分析?原理
方法
程序
IF(ABS(XXL).LT.0.0001)THEN
SLO=1.0E-6
ELSE
SLO=(YN(J2)-YN(J1))/XXL !斜率
ENDIF
IF(IN.EQ.1.OR.SLO.LE.0)GOTO11
C 第一条外边坡线为水平,同时圆弧与其有两个交点,不接受,继续搜索
IF(R.LT.D.OR.R.EQ.D)GOTO11 !无交点
NLOW=LNUM(I)
NUPP=NLOW
RETURN
20 IF(I1.NE.0)GOTO14
NLOW=LNUM(I)
I1=I1+1
11 CONTINUE
IF(I.EQ.1.AND.R.LT.R1.AND.R.GT.R2)RETURN
I=I+1
IF(I.GT.LSL)GOTO14
GOTO12
14 NUPP=LNUM(I)
15 CONTINUE
RETURN
END
SUBROUTINE DD(X1,Y1,X2,Y2,XC,YC,D,IN)
REAL X1,Y1,X2,Y2,XC,YC,D
REAL XXL,AK,AKF,X,Y,XT,YT
INTEGER IN
D=9999,
IN=0
XXL=X2-X1
IF(ABS(XXL).LT.0.001)XXL=1.0E-6
IF(ABS((Y2-Y1)/XXL).LT.0.001)RETURN
AK=(Y2-Y1)/XXL
AKF=1./AK
X=(AK*X1-Y1+AKF*XC+YC)/(AK+AKF)
Y=AK*X-AK*X1+Y1 !从圆心作垂直与线段的线,获与线段的交点
XT=(X-X1)*(X-X2)
YT=(Y-Y1)*(Y-Y2)
IF(XT.LT.0.OR.YT.LT.0)GOTO12
IN=1 !交点不在线段内
RETURN
12 D=SQRT((X-XC)**2+(Y-YC)**2) !线段与圆心的垂直距离
RETURN
END
3,例题
数据文件CIRCLE.DAT
0 0,-24.235 -35.164 -42.046
21
9
1 2 3
2 3 1
3 4 1
2 5 3
5 6 2
6 7 2
第11章 程序设计 343
7 8 2
5 9 3
9 10 3
3
1 2 3
10
1,000 -25.000
2 -30.000 -25.000
3 -50.000 -35.000
4 -70.000 -35.000
5 -40.000 -27.000
6 -50.000 -29.000
7 -54.000 -31.000
8 -70.000 -31.000
9 -52.000 -24.000
10 -70.000 -24.000
计算成果CIRCLE.C
THE ABSCISSA VALUES OF THE UPPER AND LOWER
NO,X Y ALF
1 -51.522050 -35.000000 1.094647
2 -50.955680 -33.901810 1.094647
3 -49.822940 -31.978910,969015
4 -48.690200 -30.491210,863838
5 -47.557460 -29.278430,770531
6 -46.424730 -28.265840,685123
7 -45.291990 -27.411070,605365
8 -44.159260 -26.687410,529825
9 -43.026520 -26.076840,457513
10 -41.893780 -25.566700,387709
11 -40.761050 -25.147810,319849
12 -39.628310 -24.813490,253490
13 -38.495570 -24.558860,188260
14 -37.362830 -24.380470,123833
15 -36.230100 -24.275990,059923
16 -35.097360 -24.244130 -.003743
17 -33.964630 -24.284500 -.067425
18 -32.831890 -24.397580 -.131382
19 -31.699150 -24.584810 -.195887
20 -30.566420 -24.848610 -.261233
21 -30.000050 -25.000020 -.261233
11,3,3 构造任意形状滑裂面的程序SPLINE.FOR
1,说明
本程序刘德贵1980用来计算已知第二种边界条件(曲线两端点处的二阶导数值)
三次样条函数的一阶导数值并利用样条函数对一元函数进行插值进而构筑包含直线段的复合滑裂面对于滑裂面上控制点区间为直线段时则直接线性插值而对于非直线段区间则运用上述介绍的方法进行三次样条插值滑裂面被NS1个点分为NS1?1段从上交点向下交点编号为1,2,,NS1?1程序用LS代表此NS1?1段中为直线段的线段总数在数组
LSA(I)(I=1,2,,LS)中存入这些直线段的编号
LS,滑裂面上NS1?1条线段中为直线段的总数
LSA(LS)滑裂面上线段为直线段的编号
344 土质边坡稳定分析?原理
方法
程序
AI(I,4)滑裂面上控制点处四个与样条插值有关的数
1
,κ
2
,κ
3
,κ
4
κ
A(I)滑裂面上控制点处的一阶导数值
j
m
B(I)一维数组工作单元
2,源程序
C SUBROUTINE SP1(LS,LSA,X2,Y2)
DIMENSION KQ2(50),X1(50),Y1(50),AI(50,4),LSA(60),A(50),B(50)
DIMENSION X2(80),Y2(80)
COMMON /A23/NS1,KQ2,X1,Y1/A13/IDD,IWR1,IWR2,IWR3
COMMON/WALL/IWALL,GWALL,HMW,EWALL,ETA
INTEGER*4 OPTION(6),OP1(6)
COMMON/OPP/OPTION,OP1
OPEN(5,FILE=' ',STATUS='UNKNOWN')
OPEN(6,FILE=' ',STATUS='UNKNOWN')
READ(5,*)NS1
DO I=1,NS1
READ(5,*)KQ2(I),X1(I),Y1(I)
ENDDO
C 如果滑裂面上控制点总数NS1=2,则按直线段处理并计算滑裂面与土条侧面边界线交点的x,y值;
READ(5,*)LS
DO I=1,NS1-1
LSA(I)=0
ENDDO
IF(LS.NE.0)READ(5,*)(LSA(I),I=1,LS)
LSA(LS+1)=NS1
C 滑裂面上第NS1个控制点对应的土条侧面边界线的编号
C X1(I) Y1(I)滑裂面上第I个控制点的x y值
N1=0
DO 100 IL=1,LS+1
L1=N1+1
N1=LSA(IL)
N2=N1-1
IF(LS.EQ.0)GOTO 83 ! 滑裂面上没有直线段
IF(N1.GE.NS1) GOTO 40
AI(N1,1)=0,! 计算直线段左端点处κ
1
,κ
2
,κ
3
,κ
4
AI(N1,2)=0,
AI(N1,3)=Y1(N1+1)-Y1(N1)
AI(N1,4)=Y1(N1)
40 IF(N1-L1.EQ.0)GOTO 100 ! 滑裂面上没有直线段
83 B(L1)=1.5*(Y1(L1+1)-Y1(L1))/(X1(L1+1)-X1(L1))
A(L1)=-0.5 ! 插值区间左端点处的一阶导数值
CN=3.0*(Y1(N1)-Y1(N2))/(X1(N1)-X1(N2))
IF(N1-L1.EQ.1)GOTO 55 ! 插值区间只具有左右两个端点
L2=L1+1
DO 50 J=L2,N2
H1=X1(J)-X1(J-1)
H=X1(J+1)-X1(J)
AF=H1/(H+H1) ! 计算α值
BT=3.*((1.-AF)*(Y1(J)-Y1(J-1))/H1+AF*(Y1(J+1)-Y1(J))/H)
!计算β值
A(J)=-AF/(2.+(1.-AF)*A(J-1)) ! 计算插值区间内插值节点处的一阶导数值m
j
50 B(J)=(BT-(1.-AF)*B(J-1))/(2.+(1.-AF)*A(J-1))
55 A(N1)=(CN-B(N2))/(2.+A(N2))
60 A(N2)=A(N2)*A(N2+1)+B(N2) ! 计算插值区间右端点处的一阶导数值m
n
第11章 程序设计 345
N2=N2-1
IF(N2.GE.L1)GOTO 60
N2=N1-1
DO 80 I=L1,N2
H=X1(I+1)-X1(I)
AI(I,4)=Y1(I) ! 计算插值区间内插值节点处的κ
1
,κ
2
,κ
3
,κ
4
AI(I,3)=A(I)*H
AI(I,2)=3.*(Y1(I+1)-Y1(I))-(2.*A(I)+A(I+1))*H
80 AI(I,1)=2.*(Y1(I)-Y1(I+1))+(A(I)+A(I+1))*H
100 CONTINUE
DO 200 I=1,NS1-1
L1=KQ2(I)
N1=KQ2(I+1)
STEP=1./(N1-L1)
SX=X1(I+1)-X1(I)
S1=0,
DO 200 J=L1,N1
X2(J)=X1(I)+S1*SX ! 计算滑裂面与土条侧面边界线交点的X Y值
Y2(J)=AI(I,4)+S1*(AI(I,3)+S1*(AI(I,2)+S1*AI(I,1)))
200 S1=S1+STEP
NW=KQ2(NS1)
WRITE(6,*)'N= ',NW
WRITE(6,105)(M,X2(M),Y2(M),M=1,NW)
105 FORMAT(1X,I5,2F10.3)
STOP
END
3,光滑曲线滑裂面例题
(1) 数据文件Spline-1.dat
5
1 -51.522 -35.000
6 -46.314 -28.235
10 -41.685 -25.511
15 -35.898 -24.224
20 -30.000 -25.000
0
计算成果Spline-1.c
N= 20
1 -51.522 -35.000
2 -50.480 -33.466
3 -49.439 -31.978
4 -48.397 -30.580
5 -47.356 -29.317
6 -46.314 -28.235
7 -45.157 -27.280
8 -44.000 -26.545
9 -42.842 -25.975
10 -41.685 -25.511
11 -40.528 -25.107
12 -39.370 -24.765
13 -38.213 -24.495
14 -37.055 -24.311
15 -35.898 -24.224
16 -34.718 -24.243
17 -33.539 -24.352
18 -32.359 -24.530
19 -31.180 -24.754
20 -30.000 -25.000
346 土质边坡稳定分析?原理
方法
程序
4,曲线-直线组合滑裂面例题
(1) 数据文件Spline-2.dat
5
1 -51.522 -35.000
6 -46.314 -28.235
10 -41.685 -25.511
15 -35.898 -24.224
20 -30.000 -25.000
1
3
(2) 计算成果Spline-2.c
N= 20
1 -51.522 -35.000
2 -50.480 -33.459
3 -49.439 -31.965
4 -48.397 -30.565
5 -47.356 -29.306
6 -46.314 -28.235
7 -45.157 -27.300
8 -44.000 -26.583
9 -42.842 -26.011
10 -41.685 -25.511
11 -40.528 -25.511
12 -39.370 -25.511
13 -38.213 -25.511
14 -37.055 -25.511
15 -35.898 -25.511
16 -34.718 -25.409
17 -33.539 -25.307
18 -32.359 -25.204
19 -31.180 -25.102
20 -30.000 -25.000
11,3,4 计算土条重量和物理力学参数的程序 PROFILE.FOR
1,说明
本子程序根据11.2.1介绍的方法计算每个土条的重量条底的强度指标和孔隙水压力
2,源程序
DIMENSION SH(20),M4(20),SHNE(20),M4NE(20)
DIMENSION XU(80),HSL(80),X3(80),Y3(80)
COMMON/CRU/RUS(10),LINES(80),YSF(80)/B3/IKD
COMMON/A1/IN,IN1,PF(20),PC(20),PDW(20),PDS(20),IPH,ICPH(40,3)
COMMON/A2/NN,NJ(20),NHY,DHRP,WL(80),BDA(80),CA(80),FAA(80)
#,LRU(20)
COMMON/A3/RW,RU1,UWL/SWEAK/IWEAK,KQ3(50),LWK,IC3(40,3)
COMMON/A6/X(80),Y(80),XN(99),YN(99),X2(80),Y2(80),IC(80,3),N
COMMON/A7/W(80),RU(80),C(80),F(80)/A8/MYK,KXSX,MTE,KL
COMMON/A15/UWL1,NCOUN,MLA(80),PF1(20),PC1(20),MID(20),PNOR(20)
#,NHT,MRD
INTEGER*4 OPTION(6),OP1(6)
COMMON/OPP/OPTION,OP1
OPEN(5,FILE=' ',STATUS='UNKNOWN')
OPEN(6,FILE=' ',STATUS='UNKNOWN')
C IWR5=0,1 分别相应顶部有无拉力缝YTENSION=拉力缝底Y坐标
C DS,CX,CY=滑弧深度圆心坐标
第11章 程序设计 347
READ(5,*) N ! 土条总数
READ(5,*) RW,UWL ! 水容重坡外水位
READ(5,*) IN1 ! 线段总数
DO 302 I=1,IN1
READ(5,*) (IC(I,J),J=1,3) !线段信息
302 CONTINUE
READ(5,*)NN ! 点总数
READ(5,*) (L,XN(I),YN(I),I=1,NN) ! 点坐标
DO I=1,N
READ(5,*) I,X(I),Y(I)
ENDDO
READ(5,*)IPH
IF(IPH.EQ.0)GOTO 18
DO I=1,IPH
READ(5,*) (ICPH(I,J),J=1,3) !线段信息
ENDDO
18 READ(5,*)IN,(L,PF(I),PC(I),PF1(I),PC1(I),PDW(I),PDS(I),
%LRU(I),I=1,IN)
CONTINUE
C 对土条信息重新处理保证IC(I,1)的x坐标比IC(I,2)的小
DO 37 I=1,IN1
M1=IC(I,1)
M2=IC(I,2)
IF(XN(M1).LT.XN(M2))GO TO 37
IC(I,1)=IC(I,2)
IC(I,2)=M1
37 CONTINUE
IF(IPH.EQ.0)GO TO 672
C REARRANGE ICPH(I,J) IN THE SAME WAY AS IC(I,J)
DO 38 I=1,IPH
M1=ICPH(I,1)
M2=ICPH(I,2)
IF(XN(M1).LT.XN(M2))GO TO 38
ICPH(I,1)=ICPH(I,2)
ICPH(I,2)=M1
38 CONTINUE
672 CONTINUE
GTENSION=0
C FIND THE LINES WHICH INTERCEPT THE CONSIDERED SLICE
MX=IC(1,1)
YYY=YN(MX)
DO 813 I=1,IN1
DO 813 II=1,2
IY=IC(I,II)
IF(YN(IY).GT.YYY)YYY=YN(IY)
813 CONTINUE
YYY=YYY+100,
X3(1)=X2(1)
Y3(1)=Y2(1)
DO 211 K=1,N
MF=0
QN=0,
W(K)=0,
DHR=0
WL(K)=YYY
348 土质边坡稳定分析?原理
方法
程序
LPHS=IC(IN1,3)
IF(IPH.EQ.0)GO TO 739
DO 735 I1=1,IPH
M1=ICPH(I1,1)
M2=ICPH(I1,2)
IF(X(K).LT.XN(M1))GO TO 735
IF(X(K).GE.XN(M2))GO TO 735
IF(ABS(XN(M2)-XN(M1)).LT.0.001)GOTO 735
WL(K)=YN(M1)+(X(K)-XN(M1))*(YN(M2)-YN(M1))/(XN(M2)-XN(M1))
LPHS=ICPH(I1,3)
735 CONTINUE
739 CONTINUE
PP=0,
CT=0,
FT=0,
K1=0
DO 207 I1=1,IN1
M1=IC(I1,1)
M2=IC(I1,2)
M3=I1
IF(X(K).LT.XN(M1))GO TO 207
IF(X(K).GE.XN(M2))GO TO 207
SH2=YN(M1)+(X(K)-XN(M1))*(YN(M2)-YN(M1))/(XN(M2)-XN(M1))
152 FORMAT(1X,'**',I5,2F15.6)
DSH2=ABS((Y(K)-SH2))
DSH1=SH2-Y(K)
C ASSURE THAT THE STRENGTH PARAMETERS OF THE UNDERLYING SOIL
IF(DSH2.GT.0.00005.AND.DSH1.GT.0.001)GO TO 207
K1=K1+1
SH(K1)=SH2
M4(K1)=M3
207 CONTINUE
IF(K1.NE.1)GO TO 707
LX=M4(1)
GO TO 708
707 LX=M4(1)
HN=SH(1)
DO 45 L1=2,K1
IF(SH(L1).GT.HN)GO TO 45
LX=M4(L1)
HN=SH(L1)
45 CONTINUE
708 M1=IC(LX,1)
M2=IC(LX,2)
DELX=XN(M2)-XN(M1)
IF(ABS(DELX).GT.0.001)BBK=(YN(M2)-YN(M1))/(XN(M2)-XN(M1))
X3(K)=X2(K)
Y3(K)=YN(M1)+(X3(K)-XN(M1))*BBK
139 FORMAT(1X,'K1,LX,M2,M1,XM2,XM1',4I5,2F10.3)
LINES(K)=LX !96.12.29
47 DO 46 L1=1,K1
LSB=M4(L1)
M4(L1)=IC(LSB,3)
46 CONTINUE
IF(WL(K).GT.Y(K))GO TO 734
DO 530 ISP=1,K1
SDIF=ABS(WL(K)-SH(ISP))
IF(SDIF.LT.0.00001)GO TO 734
第11章 程序设计 349
530 CONTINUE
K1=K1+1
SH(K1)=WL(K)
M4(K1)=LPHS
734 J=1
27 MSP=1
DO 209 I=1,K1
IF(J.EQ.I)GO TO 209
IF(SH(J).LT.SH(I))GO TO 209
MSP=MSP+1
209 CONTINUE
SHNE(MSP)=SH(J)
M4NE(MSP)=M4(J)
IF(J.GE.K1)GO TO 229
J=J+1
GO TO 27
229 K1=K1+1
SHNE(K1)=Y(K)
YSF(K)=SHNE(1)
M4NE(K1)=M4NE(K1-1)
SH3=SHNE(K1)-SHNE(1)
IF((K1-1).LE.0)GO TO 641
DDH=0,
IF(SHNE(1).GT.UWL)W(K)=RW*(SHNE(1)-UWL)
DO 208 I=2,K1
M5=M4NE(I-1)
WL1=WL(K)-0.0001
IF(SHNE(I-1).LT.WL1)PDNO=PDW(M5)
IF(SHNE(I-1).GE.WL1)PDNO=PDS(M5)
DH=SHNE(I)-SHNE(I-1)
DWK=DH*PDNO
W(K)=W(K)+DWK
42 CT=CT+DH*PC(M5)
IF(PC(M5).LT.0.00001.AND.PF(M5).LT.0.00001)GO TO 589
FT=FT+DH*PF(M5)
GO TO 208
589 DDH=DDH+DH
DHR=DHR+DH
208 CONTINUE
DHRP=SH3-DHR
M7=M4NE(K1)
IF(MF.NE.0)M7=IC3(MF,3)
MLA(K)=M7
PP=0
XU(K)=SHNE(1)
HSL(K)=SH3
IF(IPH.EQ.0)GO TO 736
IF(Y(K).LT.WL(K))GOTO 590
WLM=WL(K)
IF(K.EQ.1.OR.K.EQ.N)GO TO 590
IF(LRU(M7).EQ.0.AND.SHNE(1).LE.UWL)PP=(Y(K)-WLM)*RW
IF(LRU(M7).EQ.0.AND.SHNE(1).GT.UWL)PP=(Y(K)-UWL)*RW
590 CONTINUE
WT1=Y(K)-UWL
IF(WT1.LT.0)WT1=0
W(K)=W(K)-WT1*RW
PP=PP-WT1*RW
736 IF(ABS(SH3).LT.0.001.OR.ABS(W(K)).LT.0.001)GO TO 641
350 土质边坡稳定分析?原理
方法
程序
IF(LRU(M7).GE.0)THEN
RU(K)=PP/W(K)
GOTO 961
ENDIF
IF(Y(K).GT.WL(K).OR.IPH.EQ.0)THEN
RU(K)=RUS(M7)
ELSE
RU(K)=0
ENDIF
961 BDA(K)=W(K)/SH3
CA(K)=CT/SH3
IF(ABS(SH3-DDH).LT.0.0001)GO TO 587
FAA(K)=FT/(SH3-DDH)
GO TO 586
587 FAA(K)=0,
586 CONTINUE
GOTO 55
641 RU(K)=0,
W(K)=0
BDA(K)=0,
CA(K)=0,
FAA(K)=0,
55 F(K)=PF(M7)
C(K)=PC(M7)
21 CONTINUE
211 CONTINUE
IF(ABS(W(N)).GT.0.001) GO TO 256
W(N)=0.001
100 FORMAT(I3/(7F11.6))
256 CONTINUE
C(1)=C(2)
F(1)=F(2)
265 WRITE(6,251)
252 FORMAT(1X,'NUMBER OF SLICES=',I5)
WRITE(6,252)N
WRITE(6,258)
251 FORMAT(1X,'DATA FOR EACH SLICE:')
258 FORMAT(T5,'NO,',T14,'X',T24,'Y',T34,'W',T44,'RU',T54,'C',
@T64,'F',T74,'BDA',T84,'WL')
WRITE(6,253)(I,MLA(I),X(I),Y(I),W(I),RU(I),C(I),F(I),BDA(I),WL(I)
%,I=1,N)
253 FORMAT(1X,2I5,8F10.3)
104 FORMAT(1X,'PORE PRESSURE',F10.3)
WRITE(6,254)
254 FORMAT(1X,'NOTE:W=WEIGHT PER UNIT WIDTH,RU=COEF,OF PORE PRES.')
WRITE(6,255)
255 FORMAT(6X,'C=COEFFICIENT OF COHESION,F=COEFFICIENT OF FRICTION')
WRITE(6,259)
259 FORMAT(6X,'BDA=AVERAGE UNIT WEIGHT,WL=Y VALUE OF PHREATIC LINE')
END
3,例题
(1) 数据文件PROFILE.DAT
21
9.8,0.0
9
1 2 3
2 3 1
第11章 程序设计 351
3 4 1
2 5 3
5 6 2
6 7 2
7 8 2
5 9 3
9 10 3
10
1,000 -25.000
2 -30.000 -25.000
3 -50.000 -35.000
4 -70.000 -35.000
5 -40.000 -27.000
6 -50.000 -29.000
7 -54.000 -31.000
8 -70.000 -31.000
9 -52.000 -24.000
10 -70.000 -24.000
1 -51.522050 -35.000000
2 -50.955680 -33.901810
3 -49.822940 -31.978910
4 -48.690200 -30.491210
5 -47.557460 -29.278430
6 -46.424730 -28.265840
7 -45.291990 -27.411070
8 -44.159260 -26.687410
9 -43.026520 -26.076840
10 -41.893780 -25.566700
11 -40.761050 -25.147810
12 -39.628310 -24.813490
13 -38.495570 -24.558860
14 -37.362830 -24.380470
15 -36.230100 -24.275990
16 -35.097360 -24.244130
17 -33.964630 -24.284500
18 -32.831890 -24.397580
19 -31.699150 -24.584810
20 -30.566420 -24.848610
21 -30.000050 -25.000020
0
3
1 38.000,000,000,000 19.500 19.500 0
2 23.000 5.300,000,000 19.500 19.500 0
3 20.000 7.200,000,000 19.500 19.500 0
(2) 计算成果PROFILE.C
DATA FOR EACH SLICE,
NUMBER OF SLICES= 21
NO,SOIL X Y W RU C F BDA WL
1 1 -51.522 -35.000,000,000,000 38.000,000 76.000
2 1 -50.956 -33.902 21.415,000,000 38.000 19.500 76.000
3 1 -49.823 -31.979 57.185,000,000 38.000 19.500 76.000
4 1 -48.690 -30.491 75.151,000,000 38.000 19.500 76.000
5 1 -47.557 -29.278 87.756,000,000 38.000 19.500 76.000
6 2 -46.425 -28.266 96.457,000 5.300 23.000 19.500 76.000
7 2 -45.292 -27.411 102.081,000 5.300 23.000 19.500 76.000
8 2 -44.159 -26.687 105.148,000 5.300 23.000 19.500 76.000
352 土质边坡稳定分析?原理
方法
程序
9 3 -43.027 -26.077 106.010,000 7.200 20.000 19.500 76.000
10 3 -41.894 -25.567 104.914,000 7.200 20.000 19.500 76.000
11 3 -40.761 -25.148 102.038,000 7.200 20.000 19.500 76.000
12 3 -39.628 -24.813 97.513,000 7.200 20.000 19.500 76.000
13 3 -38.496 -24.559 91.434,000 7.200 20.000 19.500 76.000
14 3 -37.363 -24.380 83.868,000 7.200 20.000 19.500 76.000
15 3 -36.230 -24.276 74.862,000 7.200 20.000 19.500 76.000
16 3 -35.097 -24.244 64.439,000 7.200 20.000 19.500 76.000
17 3 -33.965 -24.285 52.607,000 7.200 20.000 19.500 76.000
18 3 -32.832 -24.398 39.358,000 7.200 20.000 19.500 76.000
19 3 -31.699 -24.585 24.663,000 7.200 20.000 19.500 76.000
20 3 -30.566 -24.849 8.475,000 7.200 20.000 19.500 76.000
21 3 -30.000 -25.000,001,000 7.200 20.000,000 76.000
NOTE:W=WEIGHT PER UNIT WIDTH,RU=COEF,OF PORE PRES,
C=COEFFICIENT OF COHESION,F=COEFFICIENT OF FRICTION
BDA=AVERAGE UNIT WEIGHT,WL=Y VALUE OF PHREATIC LINE
11,3,5 计算安全系数的程序S.FOR
1,说明
本程序包括一个主程序和三个子程序子程序MP实现垂直条分法的计算过程包括陆军工程师团法子程序SAFE2实现瑞典法毕肖普法和简化法的计算过程这三个方法和相应计算公式已在第3章中曾作了详细介绍通过一个整型变量OPTION来控制计算方法通过变量ID?C来控制滑面形状(圆弧或非圆弧)对非圆弧滑裂面通过IWALL=0和=1
分别实现垂直条分法计算安全系数和土压力的功能本源程序中将水容重RW设为9.8也就是要求使用kN?m的量纲体系控制滑面几何形状的点有两套一套为界面在土条底的坐标用X2(I) Y2(I)代表土条总数为NS另一套为土条底中点的坐标用X(I) Y(I)来代表在输入W(I) RU(I)等一系列物理力学参数时都是相应土条中点的但是加上了左右两个端点的相应物理量因此读入了总数为N=NS+1行数据这些信息实际上就是程序PROFILE.FOR的输出成果在使用程序MP实现垂直条分法时程序首先使用子程序
SAFE2中简化法的计算安全系数作为初值然后进行第3节介绍的稳定分析方法的计算子程序SAFE2在进行Bishop法计算时同样用简化法计算安全系数初值读者还可从子程序MP中发现通过式(2.43)计算简化法1即λ初值的过程
2,源程序
REAL AB1
INTEGER N,NS,OPTION
COMMON/AS/ASP
COMMON/CIRCLE/ID_C,CX,CY,DS
C OPTION=0,SPENCER法
C OPTION=1,BISHOP法
C OPTION=2,瑞典法
C OPTION=3,工程师团法
C OPTION=4,简化法
C ID_C=0,任意形状滑裂面ID_C不为零圆弧滑裂面
REAL X2(80),Y2(80),W(80),X(80),Y(80),C(80),F(80),WL(80)
$,FXO(80),FX(80),BN(80)
REALBDA(80),FQUH(80),FDIS(80),YSF(80),RU(80),ALF(80),RW
WRITE(6,*)'输入数据文件名'
OPEN(5,FILE='',STATUS='UNKNOWN')
第11章 程序设计 353
WRITE(6,*)'输入计算成果文件名'
OPEN(6,FILE='',STATUS='UNKNOWN')
RW=9.8
AB1=0.0
WRITE(6,*)'边坡稳定分析'
READ(5,*)OPTION,ID_C
IF(ID_C.NE.0)READ(5,*)CX,CY,DS
READ(5,*)N!土条总数
NS=N-1
READ(5,*)(X(I),I=1,N)
READ(5,*)(Y(I),I=1,N)
READ(5,*)(ALF(I),I=1,N)
N=NS+1
DO302I=1,N
READ(5,*)W(I),RU(I),C(I),F(I),BDA(I),WL(I),FQUH(I),FDIS(I),YSF(I)
F(I)=TAN(F(I)*3.14159/180.)
302 CONTINUE
WRITE(6,*)(ALF(I),I=1,N)
!读入两端点和土条底中点的数据,其中涉及力的物理量均以单宽土条计
!X2,Y2=坐标,W=土条重量,RU=孔压系数,C=粘聚力=,F=摩擦系数,BDA=土容重,
!WL=浸润线Y坐标,FQUH,FDIS,YSF=水平地震力大小,Y坐标,垂直水平地震力大小
IF(OPTION.EQ.3)THEN
READ(5,*)ASP
WRITE(6,*)'陆军工程师团法平均坝坡=',ASP
ASP=ASP*3.14159/180,
ENDIF
WRITE(6,709)RW
709 FORMAT(1X,'水容重=',F10.3//)
WRITE(6,702)
702 FORMAT(T25,'DATAFORASSUMEDSIDEFORCEFUNCTION')
WRITE(6,703)
703 FORMAT(T25,'************************************'//)
!计算土条底中点的坐标和条底倾角
DO208I=2,NS
X2(I)=(X(I)+X(I+1))/2,
Y2(I)=(Y(I)+Y(I+1))/2,
208 CONTINUE
X2(N)=X(N)
Y2(N)=Y(N)
X2(1)=X(1)
Y2(1)=Y(1)
WRITE(6,716)
716 FORMAT(T5,'NO.',T14,'X',T24,'Y',T34,'W',T44,'RU',T54,'C',
# T64,'F',T74,'BDA',T84,'WL',T94,'QH',T104,'YE')
DO901I=1,N
WRITE(6,718)I,X(I),Y(I),W(I),RU(I),C(I),F(I),BDA(I),WL(I)
#,FQUH(I),FDIS(I)
901 CONTINUE
718 FORMAT(1X,I5,11F10.3)
WRITE(6,719)
719 FORMAT(1X,/'NOTE:')
WRITE(6,720)
720 FORMAT(5X,'W=土条的单宽重量')
WRITE(6,721)
721 FORMAT(5X,'RU=孔隙水压力系数')
WRITE(6,722)
722 FORMAT(5X,'C=粘聚力')
354 土质边坡稳定分析?原理
方法
程序
WRITE(6,723)
723 FORMAT(5X,'F=摩擦系数')
WRITE(6,724)
724 FORMAT(5X,'BDA=土条的平均容重')
WRITE(6,725)
725 FORMAT(5X,'WL=浸润线的Y坐标')
C WRITE(6,701)
C WRITE(6,*)'迭代过程'
SELECTCASE(OPTION)
CASE(0)
WRITE(6,701)
701 FORMAT(/10X,'SPENCERMETHOD')
CALLSAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,4)
IWALL=0
DO301I=1,80
FXO(I)=0
FX(I)=1
301 CONTINUE!TAN(BETA(I))=FXO(I)+ALAM*FX(I)
KXYX=0
CALLMP(N,KXYX,KS,AB1,DQ,DM,BN,FX,FXO,
& C,F,X2,Y2,X,Y,W,RU,FQUH,WL,FDIS,RW,DB,DF)
CASE(1)
WRITE(6,704)
704 FORMAT(/20X,'BISHOP法'/)
CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION)
CASE(2)
WRITE(6,705)
705 FORMAT(/20X,'瑞典法'/)
CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION)
CASE(3)
WRITE(6,731)
731 FORMAT(/10X,'工程师团法')
CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,4)
IWALL=0
DO311I=1,80
FXO(I)=0
FX(I)=1
311 CONTINUE!TAN(BETA(I))=FXO(I)+ALAM*FX(I)
KXYX=1
CALL MP(N,KXYX,KS,AB1,DQ,DM,BN,FX,FXO,
& C,F,X2,Y2,X,Y,W,RU,FQUH,WL,FDIS,RW,DB,DF)
CASE(4)
WRITE(6,789)
789 FORMAT(/20X,'简化法法'/)
CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION)
CASE DEFAULT
ENDSELECT
END
SUBROUTINE SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,ID)
REAL X2(80),F(80),FQUH(80),W(80),RU(80),ALF(80),X(80),C(80)
REAL AB1,FDIS(80)
REAL RO1,DM,RM,BX1,CDN,DG(80),G3,DX1,G,GX
REAL O,S1,FRIC,S2,S,DRM,P1,P2,DDM,DBX1,DCDN
COMMON/CIRCLE/ID_C,CX,CY,DS
INTEGER NS
IWALL=0
第11章 程序设计 355
RO1=1.0
NS=N-1
KL=2
DM=0,
RM=0,
BX1=0,
CDN=0,
DG(2)=0.
G3=0,
DO93I=2,NS
DX1=X2(I)-X2(I-1)
DG(I+1)=DG(I)+(F(I)-F(I+1))*ALF(I)
93G3=G3+(F(I)*ALF(I)+DG(I))*DX1
G3=G3/(X(N)-X(1))
DO202I=2,NS
IF(ID_C.EQ.0)THEN
RO1=1
ELSE
RO1=(FDIS(I)-CY)/(DS-CY)
ENDIF
G=ALF(I)
GX=G*F(I)-G3+DG(I)
DX1=X2(I)-X2(I-1)
O=FQUH(I)*SIN(G)
S1=W(I)*(COS(G)-RU(I)/COS(G))
FRIC=F(I)
S2=(S1-O)*FRIC
S=C(I)/COS(G)
DRM=S+S2
DRM=DRM*DX1
P1=W(I)*SIN(G)
P2=FQUH(I)*RO1
DDM=(P1+P2)*DX1
DM=DM+DDM
RM=RM+DRM
DBX1=-DRM-SIN(G)*GX*W(I)*DX1
BX1=BX1+DBX1
DCDN=DRM*GX
CDN=CDN+DCDN
202 CONTINUE
AB1=-BX1/DM+CDN/BX1!通过简化法求得AB1的迭代初值
IF(ID.EQ.4)WRITE(6,*)'简化法安全系数=',AB1
IF(ID.EQ.2)AB1=RM/DM
IF(ID.EQ.2)WRITE(6,*)'瑞典法安全系数=',AB1
IF(ID.NE.1)RETURN
WRITE(6,*)'迭代过程'
LM=1
41DT=0,
DM=0,
RM=0,
DO21I=2,NS
G=ALF(I)
DX1=X2(I)-X2(I-1)
S1=W(I)*(1-RU(I))*F(I)
S2=S1+C(I)
S=COS(G)*(1+F(I)*TAN(G)/AB1)
DRM=DX1*S2/S
RM=RM+DRM
356 土质边坡稳定分析?原理
方法
程序
P1=W(I)*SIN(G)
P2=FQUH(I)*RO1
DDM=(P1+P2)*DX1
DM=DM+DDM
DT=DT+DRM*TAN(G)*F(I)/(AB1**2+AB1*F(I)*TAN(G))
21 CONTINUE
AB2=RM/DM!BISHOP方法
IF(AB2.LT.0.OR.AB2.GE.15)THEN
WRITE(6,*)'F.SNOTREASONABLE=',AB2
AB2=9999.9
RETURN
ENDIF
WRITE(6,740)LM,AB2
740 FORMAT(1X,I5,F10.5)
IF(ABS(AB2-AB1).LT.0.0001)GOTO49
IF(LM.LT.10)GOTO9191
WRITE(6,9192)
9192 FORMAT(1X,'ITERATIONTROUBLESINBISHOPMETHOD')
GOTO49
9191 CK=1-DT/DM!牛顿迭代法
AB1=AB1-(AB1-AB2)/CK
LM=LM+1
GOTO41
49 WRITE(6,203)AB1
203 FORMAT(/5X,'BISHOP法安全系数=',F10.4)
RETURN
END
SUBROUTINE DETE(AB1,RW,DQ,DM,ALAM,KS,BN,FX,FXO,
& C,F,X2,Y2,W,RU,FQUH,YSF,WL,FDIS,N,DB,KXYX,DF)
REAL AB1,DQ,DT,DM,ALAM,DF
REAL C(80),F(80),X2(80),Y2(80),W(80),RU(80),FXO(80),FX(80)
REAL FQUH(80),YSF(80),WL(80),FDIS(80),BN(80)
INTEGER KS,N,KXYX
REAL X(80),ALF(80),G,Y(80)
REAL EXT,EXR,EYR,EMT,EMR,FMT,DMT,DEG,DX1
REAL DMG,BS,R1X,TX,QX,PT112,R1XT,POX
REAL POXT,POY,DPT112,DPX112,BSG,H3,YT(80),GF(80)
REAL GFX(80),GFY(80),BET(80),EFX(80)
REAL G1,TA,TA1,OP,FI,FI1,CE,OQ,OQ1,THO,CW
REAL SEGB,SEGB1,FI0,AXS11,OTEM,DDQ,PX112,DBS
REAL SX,SXG,DTX1,DTX2,DTX,TX1,DDP,DBS1,DBS2,EPX112
REAL AAKX1,DR1X1,DR1X2,DR1X,DR1XT
REAL DG,DDT,GAVE,BETA,DGX,OX,OY,HSL
REAL TAA,OQS2,OQS1,DG1,BT(80),EXR1,EXT1,EYR1
REAL EYT1,XCE,YCE,XRE,YRE,DR2X11,DR2X12,DR2X1
REAL GS,DPOY1,DPOY2,OY1,DDM,DDMT,PM,PMT,DDEG,DDMG
REAL OQ0,SEGB0,DGG,DMG1,DEG1,OXX
REAL TEMP1,TEMP2,GG,DDM2,DM1,EWALL,GWALL,PX
REAL DQM,DMM,DTM,DMTM,DEGM,DMGM,S1,HMW,DDM1,ETA
REAL GTENSION,HITE,WP,RW
INTEGER EYT,IWALL,NS
!本程序计算MORGENSTERN-PRICE,SPENCER和工程师团法
IWALL=0.!土压力的作用位置
HMW=0.
GTENSION=0,
HITE=0,
第11章 程序设计 357
GWALL=0.!土压力的大小
!GWALL=1计算土压力设定条件下的安全系数
!GWALL=2计算主动土压力
IF(KS.EQ.0)THEN!KS调用本程序段标志符,=0表示第1次调用
WRITE(6,705)
705 FORMAT(1X,/'VALUEOFFX0(I)INORDEROFI')
WRITE(6,706)(FXO(I),I=1,N)
706 FORMAT(1X,5F10.5)
WRITE(6,707)
707 FORMAT(1X,/'VALUEOFFX(I)INORDEROFI')
WRITE(6,708)(FX(I),I=1,N)
708 FORMAT(1X,5F10.5)
WRITE(6,726)
726 FORMAT(1X,/'迭代过程:'/)
WRITE(6,712)
712 FORMAT(T10,'NO.',T24,'DQ',T35,'DM',T55,'AB1',T65,'ALAM')
ENDIF
KS=1
DM=0,
DT=0.
EXT=GTENSION
EYT=0
EXR=0,
EYR=0,
EMT=-GTENSION*HITE/3,
EMR=0,
FMT=0,
DQ=0.
DMT=0,
DEG=0.
DMG=0,
BS=0,
R1X=0,
TX=0,
QX=0,
PT112=0,
PX112=0,
R1XT=0,
POX=0,
POXT=0.
POY=0,
DPT112=0,
DPX112=0,
BSG=0,
H3=HITE/3,
YT(1)=Y2(1)-H3
DQ=DQ-GTENSION
GF(1)=-DQ
DM=DM-GTENSION*H3!96.8.31
NS=N-1
DO303I=2,NS
X(I)=(X2(I)+X2(I-1))/2,
Y(I)=(Y2(I)+Y2(I-1))/2,
ALF(I)=ATAN((Y2(I)-Y2(I-1))/(X2(I)-X2(I-1)))
303 CONTINUE
ALF(1)=ALF(2)
ALF(N)=0,
358 土质边坡稳定分析?原理
方法
程序
X(N)=X2(NS)
Y(N)=Y2(NS)
X(1)=X2(1) !X,Y各土条的中点坐标(在滑裂面上)
Y(1)=Y2(1)
DO304I=1,N
BET(I)=ATAN(FXO(I)+ALAM*FX(I))
IF(C(I).LT.0.001.AND.F(I).LT.0.001)BET(I)=0
304 CONTINUE
401 CONTINUE
EWALL=BET(NS) !土压力的作用角度
GFX(1)=GF(1)*COS(BET(1))
GFY(1)=GF(1)*SIN(BET(1))
DO305I=2,NS
PX=X(I)-X(I-1)
DX1=X2(I)-X2(I-1)
G=ALF(I)
G1=ALF(1)
TA=BET(I)
TA1=BET(1)
OP=G-TA
FI=ATAN(F(I)/AB1)
FI1=ATAN(F(1)/AB1)
CE=C(I)/AB1
OQ=FI-G
OQ1=FI1-G1
THO=SIN(TA)*X(I)-COS(TA)*Y(I)-(SIN(TA)*X(1)-COS(TA)*Y(1))
CW=1
IF(I.GT.2)CW=(X2(I-1)-X2(I-2))/(X2(I)-X2(I-1))
IF((OQ+TA).GT.1.5.OR.(OQ1+TA1).GT.1.5.OR.TA.GT,
$ 1.5.OR.TA1.GT.1.5)GOTO317
SEGB=(1./COS(OQ+TA))**2
SEGB1=(1./COS(OQ1+TA1))**2
FMT=FMT+FQUH(I)*(Y(I)-FDIS(I))*DX1
FI0=ATAN(F(I+1)/AB1)
AXS11=ABS((FI-FI0))
OTEM=ABS(OQ+TA)*180./3.14259
IF(OTEM.LT.95.AND.OTEM.GT.85)RETURN
DDQ=(W(I)*SIN(OQ)+CE*COS(FI)/COS(G)-RU(I)*W(I)*SIN(FI)/COS(G)
& -FQUH(I)*COS(OQ))/COS(OQ+TA)!DDQ=P(X)SEC(FI-ALF-BETA)
IF(I.GT.2)GOTO307
PX112=-TAN(OQ1+TA1)*COS(TA1)**2.*(FX(1))
DBS1=TAN(OQ1+TA1)
307 DBS2=TAN(OQ+TA)
EPX112=-DBS2*COS(TA)**2.*(FX(I))
DBS=(DBS1*CW+DBS2)/(1+CW)
DBS1=DBS2
BS=BS+DBS*(BET(I)-BET(I-1))
BSG=BS+DBS*(BET(I)-BET(I-1))*(0.5*DX1)/PX
IF((-BS).GT.30.OR.BSG.GT.30)GOTO317
SX=2.71828**(-BS)!SX=S(X)*COS(FI-ALF+BETA)
SXG=2.71828**BSG
IF(I.GT.2)GOTO309
DTX1=(SIN(TA1)-TAN(G1)*COS(TA1))
309 DTX2=(SIN(TA)-TAN(G)*COS(TA))/SX
DTX=(DTX1*CW+DTX2)/(1+CW)
DTX1=DTX2
TX1=TX+DTX*PX*0.5
TX=TX+DTX*(X(I)-X(I-1))
第11章 程序设计 359
DDP=(-FQUH(I)*SIN(TA)*SIN(FI)+CE*COS(OP)*COS(FI)/
# COS(G)+W(I)*SIN(FI)*COS(TA)-RU(I)*W(I)*SIN(FI)*COS(OP)/
# COS(G))*COS(FI)/(-AB1*COS(OQ+TA)**2)
AAKX1=DDP/DDQ
IF(I.GT.2)GOTO310
DR1X1=(-SIN(FI1)*COS(FI1)*SEGB1/AB1)
310 DR1X2=(-SIN(FI)*COS(FI)*SEGB/AB1)
DR1X=(DR1X1*CW+DR1X2)/(1+CW)
DR1X1=DR1X2
DR1XT=DR1X*TX1
R1XT=R1XT+DR1XT*(BET(I)-BET(I-1))
R1X=R1X+DR1X*(BET(I)-BET(I-1))
DG=DDQ*SX*DX1
DQ=DQ+DG!DQ=DG的积分,DG=P(X)S(X),力的平衡
DDT=DG*(AAKX1-R1X)
DT=DT+DDT!DT=DDT的积分,DDT=P(X)S(X)T(X),力矩平衡
GF(I)=-DQ*SXG
GFX(I)=GF(I)*COS(BET(I))
GFY(I)=GF(I)*SIN(BET(I))
GAVE=(GFY(I)+GFY(I-1))/2,
BETA=(BET(I)+BET(I-1))/2,
DGX=GFX(I)-GFX(I-1)
OX=(GFX(I-1)+GFX(I))/2,
OY=(GFY(I-1)+GFY(I))/2,
IF(WL(I).GE.Y(I))THEN
WP=0,
ELSE
WP=0.5*(RU(I)*W(I))**2/RW
ENDIF
OXX=OX-WP
EFX(I)=OXX
IF(OXX.GT.0.)HSL=Y(I)-YSF(I)
TAA=COS(TA)*COS(TA)*FX(I)
OQS2=(BET(I+1)+BET(I))/2
OQS1=(BET(I-1)+BET(I))/2
DG1=GF(I)*SIN(G-OQS2)-GF(I-1)*SIN(G-OQS1)
BN(I)=W(I)*DX1*COS(G)+DG1-FQUH(I)*DX1*SIN(G)
BT(I)=(BN(I)-RU(I)*W(I)*DX1/COS(G))*TAN(FI)+CE*DX1/COS(G)
EXR1=BN(I)*SIN(G)-BT(I)*COS(G)
EXR=EXR+EXR1
EXT1=FQUH(I)*DX1
EXT=EXT+EXT1
EYR1=-BN(I)*COS(G)-BT(I)*SIN(G)
EYR=EYR+EYR1
EYT1=W(I)*DX1
EYT=EYT+EYT1
XCE=X(1)
YCE=Y(1)
XRE=X(I)-XCE
YRE=Y(I)-YCE
EMT=EMT+EXT1*(FDIS(I)-YCE)-EYT1*XRE!MOMENTOFBODYFORCE
EMR=EMR+EXR1*YRE-EYR1*XRE!MOMENTOFINTERNALFORCE
IF(ABS(DX1).LT.0.001)GOTO317
BN(I)=BN(I)*COS(G)/DX1
BT(I)=BT(I)*COS(G)/DX1
IF(KXYX.EQ.0)GOTO402
GOTO305
360 土质边坡稳定分析?原理
方法
程序
402 I F(I.GT.2)GOTO312
DR2X11=SEGB1*COS(TA1)**2.*(FX(1))
312 DR2X12=SEGB*COS(TA)**2.*(FX(I))
DR2X1=(DR2X11*CW+DR2X12)/(1+CW)
DR2X11=DR2X12
GS=ALF(I-1)
POX=POX-DR2X1*(ALF(I)-GS)
POXT=POXT-DR2X1*(ALF(I)-GS)*TX1
IF(I.GT.2)GOTO313
DPOY1=(COS(FI1)/(COS(G1)*COS(OQ1+TA1)))*COS(TA1)**2.*(FX(1))
313 DPOY2=(COS(FI)/(COS(G)*COS(OQ+TA)*SX))*COS(TA)**2.*(FX(I))
OY1=(DPOY1*CW+DPOY2)/(1+CW)
DPOY1=DPOY2
POY=POY+OY1*PX
DDM=DG*TX
DDMT=DG*(AAKX1*TX-R1XT)
PX112=PX112+DPX112
PT112=PT112+DPT112
DPX112=0
DPT112=0,
PM=POX-PX112
PMT=POXT+POY-PT112
DDEG=DG*PM
DDMG=DG*PMT
DM=DM+DDM
DMT=DMT+DDMT
DEG=DEG+DDEG
DMG=DMG+DDMG
IF(AXS11-0.005)305,305,315
G=ALF(I+1)
315 OQ0=FI0-G
IF((OQ0+TA).GT.1.5)GOTO317
SEGB0=(1./COS(OQ0+TA)**2)
DR2X11=SEGB0*COS(TA)**2.*(FX(I))
DPX112=-TAN(FI0+BET(I+1)-G)*COS(BET(I+1))**2*FX(I+1)
# +TAN(OQ+TA)*COS(TA)**2.*FX(I)
DPT112=DPX112*TX
305 CONTINUE
TEMP1=EPX112+PM
TEMP2=PMT-TX*TAN(OQ+TA)*TAA
IF(IWALL.EQ.2)THEN
GG=-DQ
DDM1=HMW*COS(EWALL)+TX*SX
DDM2=HMW*COS(EWALL)/SX+TX
DM1=DM
DM=DM1+GG/SX*DDM1-FMT
DDDM2=-HMW*COS(EWALL)/SX*TEMP1-TX*TEMP1+TEMP2
DGG=-DEG*(HMW*COS(EWALL)/SX+TX)
DMG1=DMG
DEG1=DEG
DMG=DMG1+GG*(-HMW*COS(EWALL)/SX*TEMP1-HMW/SX*SIN(EWALL)
# *FX(N)*COS(EWALL)*COS(EWALL)-TX*TEMP1+TEMP2)-
# DEG*(HMW*COS(EWALL)/SX+TX)
DB=-DM/DMG
ETA=1.0+DQ/GWALL/SX
EXT=EXT-GWALL*COS(EWALL)*(1-ETA)
EYT=EYT-GWALL*SIN(EWALL)*(1-ETA)
第11章 程序设计 361
EMT=EMT+GWALL*SIN(EWALL)*(X(N)-X(1))*(1-ETA)-
%GWALL*COS(EWALL)*(Y(N)-Y(1)-HMW)*(1-ETA)
GOTO317
ENDIF
IF(IWALL.EQ.1)THEN
DQM=SX*GWALL
DDM1=HMW*COS(EWALL)+TX*SX
DMM=GWALL*DDM1
DTM=-GWALL*SX*R1X
DMTM=-GWALL*SX*R1XT
DEGM=TEMP1*SX*GWALL
DMGM=TEMP2*SX*GWALL
DMGM=DMGM-GWALL*SIN(EWALL)*FX(N)*COS(EWALL)
% *COS(EWALL)*HMW
DQ=DQ+DQM !式(2.31)
DM=DM+DMM !式(2.32)
DT=DT+DTM !DT=DQ对AB1的偏导数
DMT=DMT+DMTM !DMT=DM对AB1的偏导数
DEG=DEG+DEGM !DEG=DQ对ALAM的偏导数
DMG=DMG+DMGM !DMG=DM对ALAM的偏导数
EXT=EXT-GWALL*COS(EWALL)
EYT=EYT-GWALL*SIN(EWALL)
EMT=EMT+GWALL*SIN(EWALL)*(X(N)-X(1))-GWALL*COS(EWALL)
$ *(Y(N)-Y(1)-HMW)
ENDIF
316 CONTINUE
IF(ABS(DM).GT.1.0E+20)GOTO317
IF(ABS(DMT).GT.1.0E+20)GOTO317
DM=DM-FMT !
S1=DEG*DMT-DT*DMG
IF(KXYX.EQ.1)THEN
DF=-DQ/DT
GOTO317
ENDIF
IF(ABS(S1).LE.0.01)GOTO317
DF=(DQ*DMG-DM*DEG)/S1 !式(2.33)
DB=(-DQ*DMT+DM*DT)/S1 !式(2.34)
317 RETURN
END
C
SUBROUTINE MP(N,KXYX,KS,AB1,DQ,DM,BN,FX,FXO,
& C,F,X2,Y2,X,Y,W,RU,FQUH,WL,FDIS,RW,DB,DF)
REAL AB1,DF,DQ,DM,ALAM,DB
INTEGER KS,N,KXYX
REAL X2(80),Y2(80),W(80),X(80),Y(80),C(80),F(80),WL(80)
REAL FQUH(80),FDIS(80),YSF(80),RU(80),RW
REAL GAV,PLA,FA,FB,PLA1
REAL BN(80),FX(80),FXO(80)
COMMON/AS/ASP
!KXYX力矩平衡控制符,=0时同时计算力和力矩的平衡
FA=0,
FB=0.
PLA=0,
IF(KXYX.EQ.0)THEN
GAV=(Y(N)-Y(1))/(X(N)-X(1))
DO277I=2,N
362 土质边坡稳定分析?原理
方法
程序
PLA=PLA+(X(I)-X(I-1))*(FX(I)+FX(I-1))/2,
277 CONTINUE
PLA1=(FA+FB)*(X(N)-X(1))/2,
ALAM=(GAV*(X(N)-X(1))-PLA1)/PLA !ALAM的初值,参见1.3.3节
ELSE
ALAM=ASP
ENDIF
MQ=1
15 CONTINUE
CALL DETE(AB1,RW,DQ,DM,ALAM,KS,BN,FX,FXO,
& C,F,X2,Y2,W,RU,FQUH,YSF,WL,FDIS,N,DB,KXYX,DF)
WRITE(6,*)MQ,DQ,DM,AB1,ALAM
AB1=AB1+DF
IF(KXYX.EQ.0)ALAM=ALAM+DB
IF(MQ.GE.20)GOTO44
MQ=MQ+1
IF(ABS(DF).GT.0.1E-4)GOTO15
IF(ABS(DB).GT.0.1E-3)GOTO15
44 CONTINUE
IF(MQ.GE.20)THEN
WRITE(6,*)'THEITERATIONFAILS'
RETURN
ENDIF
WRITE(6,727)
727 FORMAT(1X,//'THERESULTIS:')
WRITE(6,710)ALAM
710 FORMAT(1X,'LAMBDA=',F12.6)
WRITE(6,711)AB1
711 FORMAT(1X,'FACTOROFSAFETY=',F12.6/)
RETURN
END
3,例题
以11.3.1节的算例说明上节所示源程序的使用计算简图如图11.4所示数据文件分别为s0.dat,s1.dat,s2.dat,s3.dat,s4.dat分别相应于垂直条分法中Morgenstern?Price法
Bishop法瑞典法陆军工程师团法和简化法这些数据文件除控制符OPTION不同外其它均相同(陆军工程师团法在数据最后增加了一行输入平均坝坡数据)文件s0.dat列于后相应迭代过程如下
(1) 数据文件(s0.dat)
0 0
21
-51.522 -50.956 -49.823 -48.690 -47.557
-46.425 -45.292 -44.159 -43.027 -41.894
-40.761 -39.628 -38.496 -37.363 -36.230
-35.097 -33.965 -32.832 -31.699 -30.566
-30.000
-35.000 -33.902 -31.979 -30.491 -29.278
-28.266 -27.411 -26.687 -26.077 -25.567
-25.148 -24.813 -24.559 -24.380 -24.276
-24.244 -24.285 -24.398 -24.585 -24.849
-25.000
1.094647 1.094647,969015,863838,770531
第11章 程序设计 363
,685123,605365,529825,457513,387709
,319849,253490,188260,123833,059923
-.003743 -.067425 -.131382 -.195887 -.261233
-.261233
.000,000,000 38.000,000 76.000 0 0 0
21.415,000,000 38.000 19.500 76.000 0 0 0
57.185,000,000 38.000 19.500 76.000 0 0 0
75.151,000,000 38.000 19.500 76.000 0 0 0
87.756,000,000 38.000 19.500 76.000 0 0 0
96.457,000 5.300 23.000 19.500 76.000 0 0 0
102.081,000 5.300 23.000 19.500 76.000 0 0 0
105.148,000 5.300 23.000 19.500 76.000 0 0 0
106.010,000 7.200 20.000 19.500 76.000 0 0 0
104.914,000 7.200 20.000 19.500 76.000 0 0 0
102.038,000 7.200 20.000 19.500 76.000 0 0 0
97.513,000 7.200 20.000 19.500 76.000 0 0 0
91.434,000 7.200 20.000 19.500 76.000 0 0 0
83.868,000 7.200 20.000 19.500 76.000 0 0 0
74.862,000 7.200 20.000 19.500 76.000 0 0 0
64.439,000 7.200 20.000 19.500 76.000 0 0 0
52.607,000 7.200 20.000 19.500 76.000 0 0 0
39.358,000 7.200 20.000 19.500 76.000 0 0 0
24.663,000 7.200 20.000 19.500 76.000 0 0 0
8.475,000 7.200 20.000 19.500 76.000 0 0 0
.001,000 7.200 20.000,000 76.000 0 0 0
(2) 计算成果(S0.c)
条底倾角α值
1.094647 1.094647 9.690150E-01 8.638380E-01
7.705310E-01 6.851230E-01 6.053650E-01 5.298250E-01
4.575130E-01 3.877090E-01 3.198490E-01 2.534900E-01
1.882600E-01 1.238330E-01 5.992300E-02 -3.743000E-03
-6.742500E-02 -1.313820E-01 -1.958870E-01 -2.612330E-01
-2.612330E-01
水容重= 9.800
DATA FOR ASSUMED SIDE FORCE FUNCTION
************************************
NO,X Y W RU C F BDA WL QH YE
1 -51.522 -35.000,000,000,000,781,000 76.000,000,000
2 -50.956 -33.902 21.415,000,000,781 19.500 76.000,000,000
3 -49.823 -31.979 57.185,000,000,781 19.500 76.000,000,000
4 -48.690 -30.491 75.151,000,000,781 19.500 76.000,000,000
5 -47.557 -29.278 87.756,000,000,781 19.500 76.000,000,000
6 -46.425 -28.266 96.457,000 5.300,424 19.500 76.000,000,000
7 -45.292 -27.411 102.081,000 5.300,424 19.500 76.000,000,000
8 -44.159 -26.687 105.148,000 5.300,424 19.500 76.000,000,000
9 -43.027 -26.077 106.010,000 7.200,364 19.500 76.000,000,000
10 -41.894 -25.567 104.914,000 7.200,364 19.500 76.000,000,000
11 -40.761 -25.148 102.038,000 7.200,364 19.500 76.000,000,000
12 -39.628 -24.813 97.513,000 7.200,364 19.500 76.000,000,000
13 -38.496 -24.559 91.434,000 7.200,364 19.500 76.000,000,000
14 -37.363 -24.380 83.868,000 7.200,364 19.500 76.000,000,000
15 -36.230 -24.276 74.862,000 7.200,364 19.500 76.000,000,000
16 -35.097 -24.244 64.439,000 7.200,364 19.500 76.000,000,000
17 -33.965 -24.285 52.607,000 7.200,364 19.500 76.000,000,000
364 土质边坡稳定分析?原理
方法
程序
18 -32.832 -24.398 39.358,000 7.200,364 19.500 76.000,000,000
19 -31.699 -24.585 24.663,000 7.200,364 19.500 76.000,000,000
20 -30.566 -24.849 8.475,000 7.200,364 19.500 76.000,000,000
21 -30.000 -25.000,001,000 7.200,364,000 76.000,000,000
NOTE,
W=土条的单宽重量
RU=孔隙水压力系数
C=粘聚力
F=摩擦系数
BDA=土条的平均容重
WL=浸润线的Y坐标
SPENCER METHOD
简化法安全系数= 1.325056
迭代过程,
NO,DQ DM AB1 ALAM
1 20.968520 139.277800 1.325056 4.646409E-01
2 9.050338E-01 15.479030 1.344203 3.878648E-01
3 4.864426E-03 4.935455E-02 1.343819 3.800467E-01
THE RESULTS ARE,
LAMBDA=,380018
FACTOR OF SAFETY= 1.343821
上表中DQ和DM为按式(2.31)和式(2.32)计算的不平衡力和力矩AB1和ALAM分别为安全系数F和λ数值
相应s1.dat,s2.dat,s3.dat,d4.dat的安全系数分别为1.384 1.228 1.346和1.350
11,4 边坡稳定分析程序的考核
11,4,1 概述
随着计算机的普及推广几乎所有的岩土工程设计者都使用自己编制的或商用程序对边坡的稳定性进行分析计算鉴于计算成果对工程的经济和安全有着直接影响因此对稳定分析程序的可靠性进行评估和鉴定就成为广受关注的工作
Fredlund (1987)曾对加拿大常用的边坡稳定分析程序作了一次调查并以表格形式介绍各程序的功能和用途但未对这些程序进行考核他的调查发现获得考核的程序对同一算例的计算结果互相差异甚大1980年香港工程发展部岩土工程处 (Geotechnical Control
Office of the Engineering Development Department,Hong Kong) 曾进行过一次测验提出了29
道边坡稳定分析的考核题每道考核题都详细注明了几何断面土的力学强度指标地下水情况表面荷载和计算的滑裂面所有的考核题都是从香港的实际工程中提练出来的42
个单位或个人提交了54份答卷其中75是采用简化Jaubu法进行安全系数计算计算的成果互相接近但是由于使用的是简化Jambu法其成果都和严格方法(MorgesternPrice法
Sarma法)有较大的出入此外这次测验仅局限于对指定的滑裂面计算没有涉及最小安全系数计算Lumsdaine和Tang (1982)曾对这次测验作过详细报道香港工程发展部岩土工第11章 程序设计 365
程处随后建立了一个对各设计部门使用的边坡稳定分析程序进行审核的制度凡是提交该处审查的计算成果必须使用经审定的程序程序在报审时必须提交五道考核题的答卷包括全部的输入输出数据和使用手册
自1982年开始水利电力部规划设计院在潘家铮总工程师主持下曾对水电系统常用的土石坝边坡稳定分析和其它数值分析程序进行开发和鉴定工作原水利电力部昆明设计院根据天生桥水电站工程的实际情况设计了考核题
在国内外对边坡稳定分析程序进行考核的工作中澳大利亚计算机应用协会的一项调查工作细致并已正式记载于文献(Donald & Giam,1992)在下节中将着重介绍这一考核工作的成果
11,4,2 澳大利亚ACADS边坡稳定分析程序调查
1987年澳大利亚计算机应用协会(ACADS)对澳大利亚所使用的边坡稳定分析程序进行一次调查他们委托Monash大学的B,Donald教授和P,Giam博士主持此项工作此次调查设计的考核题五道总计十个小题向120个单位的发出了邀请有28个单位发回了计算答案东道主Donald教授还邀请了国际上在边坡稳定分析程序方面做过较多工作的研究者提供裁判程序(Referee Program)答案本书作者有幸与加拿大Saskatchwen大学Fredlund
教授和以色列工学院的Baker教授一起被邀请担任此项工作本研究成果已于1989年4月
19日在澳大利亚岩土工程协会维克多利亚分会的月会上公布,后正式发表(Donald & Giam,
1992)由于这次调查工作规模较大所获得的成果比较可靠本节详细介绍有关的考核题和答案这些内容可作为各单位和个人对自编程序考核的参考资料
1,考核题设计
(1) 考核题1本题包括四小题图11.5图11.6及表11.1和表11.2分别给出本题剖面的几何尺寸及物性指标对所有四种情况要求确定临界滑裂面(圆弧或非圆弧)以及最小安全系数
考核题1(a)为一均质边坡材料特性如表11.1所示
考核题1(b)为同样几何形状的均质边坡但材料特性如表11.2所示要求选择一个合适的拉力缝深度并假定缝内充水
考核题1(c)为非均质边坡材料特性如表11.3所示相信读者对考核题1(c)已很熟悉了因为在第2第4第10章中已多次使用
考核题1(d)情况与考核题1(c)情况相同只是土体作用一水平地震加速度系数为0.15g
366 土质边坡稳定分析?原理
方法
程序
图 11,5 ACADS考核题EX1(a)和EX1(b)
图 11,6 ACADS考核题EX1(c)和EX1(d)
表 11,1 考核题1(a)材料性质
c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
备注
3.0 19.6 20.0 1.0E4 0.25 0.65 1
#
土
表 11,2 考核题1 (b)材料性质
c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
备注
32.0 10.0 20.0 1.0E4 0.25 0.65 1
#
土
注 这一问题需设一合适的拉力缝
表 11,3 考核题1(c d)材料性质
土号 c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
1
#
土 0.0 38.0 19.5 1.0E4 0.25 0.65
2
#
土 5.3 23.0 19.5 1.0E4 0.25 0.65
3
#
土 7.2 20.0 19.5 1.0E4 0.25 0.65
(2) 考核题2本题选自Talbingo坝实例图11.7及表11.4~表11.6分别给出其剖面的第11章 程序设计 367
几何尺寸及物性指标材料1为堆石填料材料2为心墙料
考核题2(a)要求确定上游坡临界滑裂面及相应最小安全系数
考核题2(b)要求对指定圆弧计算其稳定安全系数圆弧圆心坐标和半径值参见表11.6
(3) 考核题3本题的边坡剖面几何特征和材料特性示于图11.8和表11.7地下水位假定位于软弱夹层底部地下水位以上孔压为零考核题3(a)要求确定临界滑裂面(圆弧或非圆弧)和相应最小安全系数考核题3(b)要求计算指定滑裂面(其坐标如表11.8所示)的安全系数
考核题3的数据如表11.7和表11.8
表 11,4 考核题2几何剖面数据
点号 X (m) Y (m) 点号 X (m) Y (m)
1 0 0 14 648 0
2 315.5 162 15 168.1 0
3 319.5 162 16 302.2 130.6
4 321.6 162 17 200.7 0
5 327.6 162 18 311.9 130.6
6 386.9 130.6 19 307.1 0
7 394.1 130.6 20 331.3 130.6
8 453.4 97.9 21 328.8 146.1
9 460.6 97.9 22 310.7 0
10 515 65.3 23 333.7 130.6
11 521.1 65.3 24 331.3 146.1
12 577.9 31.4 25 372.4 0
13 585.1 31.4 26 347 130.6
注 点号参考图11.7
表 11,5 考核题2材料性质
c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
堆石 0.0 45 20.4 5.0E4 0.30 0.65
过渡层 0.0 45 20.4 5.0E4 0.30 0.65
反滤层 0.0 45 20.4 4.0E4 0.30 0.65
心墙 85 23 18.1 1.0E4 0.44 0.65
注 堆石过渡层和反滤层几乎具有同样的材料性质
368 土质边坡稳定分析?原理
方法
程序
图 11,7 ACADS考核题EX2
表 11,6 考核题2(b)圆弧滑裂面信息
X (m) Y(m) 半径(m)
100.3 290.0 278.8
注 x y轴的原点位于上游堆石边坡和坝基的交点
表 11,7 考核题3材料性质
c(kN/m
2
) φ(°) γ(kN/m
3
) E(kN/m
2
) ν K
0
1
#
土 28.5 20.0 18.84 6.0E4 0.25 0.65
2
#
土 0.0 10.0 18.84 2.0E3 0.25 0.65
注 假设水面线与基础中的软弱夹层相重合(夹层和夹层上的孔压=0)不考虑水面线以上负的孔隙水压力
表 11,8 考核题3(b)滑裂面控制点坐标
X(m) Y(m)
41.85 27.75
44.00 26.50
63.50 27.00
73.31 40.00
图 11,8 ACADS考核题EX3
(4) 考核题4本题为一坡面上作用表面荷载坡内有浸润线的例子图11.9及表11.9
表11.10和表11.11分别给出本题的剖面几何尺寸物性指标浸润线及表面荷载要求确第11章 程序设计 369
定临界滑裂面和相应最小安全系数
表 11,9 考核题4材料性质
c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
1
#
土 28.5 20.0 18.84 6.0E4 0.25 0.65
2
#
土 0.0 10.0 18.84 2.0E3 0.25 0.65
表 11,10 考核题4中的外荷载
X(m) Y(m)
正应力
(kN/m
2
)
23.00 27.75 20.00
43.00 27.75 20.00
70.00 40.00 20.00
80.00 40.00 40.00
表 11,11 考核题4的浸润线数据
点号 X(m) Y(m)
1 20.0 27.75
2 43.0 27.75
3 49.0 29.8
4 60.0 34.0
5 66.0 35.8
6 74.0 37.6
7 80.0 38.4
9 84.0 38.4
注 点号参考图11.9
图 11,9 ACADS考核题EX4
(5) 考核题5考核题5如图11.10图11.17材料性质如表11.12要求计算最小安全系数
表 11,12 考核题5材料性质
c (kN/m
2
) φ (°) γ (kN/m
3
) E (kN/m
2
) ν K
0
K
x
=K
y
(m/s)
11.0 28.0 20.00 1.5E4 0.25 0.65 5.0E-08
370 土质边坡稳定分析?原理
方法
程序
图 11,10 ACADS考核题EX5
2,计算成果汇总和分析
如前所述组织者向澳大利亚各咨询公司发出了问卷调整对各题目返回的答案数目在
7~33个之间表11.13总结了收到的答案
本次调查邀请了国外三个程序提供裁判答案Donald教授也使用了自编的GWEDGEM
和EMU程序分析了这些考核题综合各咨询公司提供的答案裁判答案和Donald自编程序的答案在表11.13汇总本次考核最终提出了推荐的裁判答案
现对各题的解答情况作一简要回顾
表 11,13 返回的问卷答案统计成果
安全系数
考核题号 分析方法
均值 标准差 min.F max.F
提交答案总数
(ALL) 0.991 0.031 0.940 1.080 32
BISHOP 0.993 0.015 0.960 1.030 18
1(a)
(1.00)
JANBU 0.978 0.041 0.940 1.043 7
(ALL) 1.381 0.067 1.230 1.520 30
BISHOP 1.406 0.039 1.350 1.500 16
1(c)
(1.39)
JANBU 1.325 0.057 1.264 1.400 6
(ALL) 0.972 0.054 0.900 1.100 14 1(d)
(1.00)
(ALL) 2.003 0.101 1.830 2.201 23 2(a)
(1.90) BISHOP 2.042 0.118 1.830 2.201 13
(ALL) 2.237 0.112 1.948 2.390 23 2(b)
(2.29) BISHOP 2.204 0.013 2.180 2.220 11
(ALL) 1.369 0.168 1.191 1.820 25
BISHOP 1.559 0.146 1.400 1.820 7
JANBU 1.307 0.138 1.191 1.564 8
3(a)
(1.26?)
NON-CIRCULAR 1.295 0.108 1.191 1.564 18
(ALL) 1.292 0.064 1.197 1.450 29 3(b)
(1.34) JANBU 1.252 0.060 1.197 1.333 10
(ALL) 0.897 0.254 0.410 1.360 25
BISHOP 1.151 0.175 0.850 1.360 6
JANBU 0.984 0.198 0.712 1.287 7
4
(0.78)
NON-CIRCULAR 0.817 0.223 0.410 1.287 19
5 (ALL) 1.462 0.150 0.910 1.600 20
第11章 程序设计 371
(1.35) BISHOP 1.431 0.192 0.910 1.600 10
(1) 考核题1本题是断面相对比较简单的土质边坡包括了四个小题主要考核用优化方法来确定的临界滑裂面和最小安全系数裁判答案作者均使用Bishop简化法从表11.14计算成果看对考核题1(a) 1(c)和1(d)
裁判答案各家给出的结果完全一致图11.11给出了对考核题1采用STAB
求得的临界滑裂面和相应的安全系数图11.12给出了对考核题1各家程序求得的临界滑裂面和相应的安全系数可以看出不同的程序得到的临界滑裂面都落在一个窄的区域里相应的安全系数也基本一致这说明对于土质边坡圆弧滑裂面目前水平已可以达到不同的程序均可获得一致成果的程度自动搜索临界滑裂面的功能也已经达到实用阶段
图 11,11 总结返回的考核题EX1计算成果的临界滑裂面
(a) EX1(a) (b) EX1(b) (c) EX1(c) (d) EX1(d)
表 11,14 裁判程序答案汇总(考核题1)
考核题 1(a) 1(b) 1(c) 1(d)
推荐的裁判答案Donald 1.00 1.65~1.70 1.39 1.00
SSA(Baker) 1.00 1.50 1.39
STAB(Chen) 0.991 1.618 1.385 1.007
GWEDGEM 1.00 1.65~1.70 1.39 1.00
EMU 1.00 1.65~1.70 1.39 1.00
Fredlund 0.99 1.619~1.588 1.406 1.015
372 土质边坡稳定分析?原理
方法
程序
图 11,12 各家程序对考核题EX1的计算成果(临界滑裂面位置和最小安全系数)
(a) EX1(c); (b) EX1(d)
图例SR? Sarma法B Bishop法JG Janbu法GW GWEDGEM法J简化扬布法MP Morgensten?Price法
考核题1(b)需要用户自行设计一个拉力缝这样由于设有的拉力缝高度不同安全系数有一定的出入但各家的差距也不大H
t
均在0.5~3.5m范围内图11.11(b)示出了笔者选用
H
t
=2.0 m时充水拉力缝的位置及临界滑裂面和相应的安全系数
因此考核题1可以作为一个经典的考核题考核自编的边坡稳定分析程序的准确性
(2) 考核题2本题分为两部分第一部分是求上游坝坡的安全系数由于坝壳料为无粘聚力的砾质材料故如使用线性强度指标将获得一个非常浅的临界滑裂面而其安全系数就是摩擦系数和坝的坡度之比即F
=
tanφ / tanα
=
1.947各家得到的结论是相同的考核第11章 程序设计 373
题还要求在厚心墙内寻找一个局部极值从表11.15中可见GWEDGEM和EMU给出的局部极值为1.90而Fredlund和陈祖煜的答案是2.14~2.15图11.13(a)为STAB程序获得考核题EX2的临界滑裂面和最小安全系数图11.14(a)列出了考核题2(a)各家得到的临界滑裂面和相应安全系数
第二部分是计算指定的圆弧滑裂面的安全系数,如图11.13(b)各家计算成果偏离不大
Baker和Donald的成果均为2.29陈祖煜和Fredlund的结果分别为2.220和2.206
表 11,15 裁判程序答案汇总(考核题2~考核题5)
2(a) 2(b) 3(a) 3(b) 4 5
推荐的裁判答案 (1.95)/1.90 2.29 1.26 1.34 0.78 1.53/1.53/1.59
SSA+(baker) 1.95/ 2.35 1.24 1.37 0.65 1.53
STAB+(Chen) 1.947~2.135 2.224 1.242 1.265 0.671 1.503
GWEDGEM 1.95/1.90 2.29 1.27 1.34 0.78 1.53/1.53/1.59
EMU 1.95/1.90 2.29 1.27 1.34 0.78 1.53/1.53/1.59
Fredland 2.138 2.206 1.261 0.630 1.460
图 11,13 STAB程序获得考核题EX2的临界滑裂面和最小安全系数
(a) 考核题 EX2(a) (b) 考核题 EX2(b)
(3) 考核题3本题为典型的沿软弱夹层进行非圆弧滑裂面的计算考核题考核题3(a)
要求计算临界滑裂面考核题3(b)要求计算指定滑裂面的安全系数STAB采用
Morgenstern?Price法图11.15所示STAB程序获得考核题EX3的临界滑裂面和最小安全系数对指定滑裂面算得安全系数为1.265[考核题3(b)]其值与Fredlund的解答1.261接近但Baker和Donald等给出了稍大的结果(1.34~1.37)最终各家搜索获得的最小安全系数比较接近在1.24~1.27之间图11.16为各家给出了本次调查考核题3(a)的临界滑裂面和相应的安全系数
374 土质边坡稳定分析?原理
方法
程序
图 11,14 各家程序对考核题EX2的计算成果(临界滑裂面位置和最小安全系数)
图 11,15 STAB程序对考核题EX3计算的滑裂面和安全系数
(a)考核题 EX3(a); (b) 考核题 EX3(b)
图 11,16 各家对考核题EX3计算的临界滑裂面和最小安全系数
(4) 考核题4本例在设计方面存在一定的问题软弱夹层的上部没有标明是否与边坡第11章 程序设计 375
面相交由于夹层的指标明显低于土层导致在搜索过程中临界滑裂面不断上移实际上直正的临界滑裂面可能要越出规定的图框临界滑裂面除趾部外全部通过软弱夹层另外
Bishop法在此例中不适用因控制滑裂面通过软弱夹层不可能是圆弧形的因此表11.15
提供的答案仅供参考
(5) 考核题5本文为一个地基快速开挖的例题考核要求采用三种不同的孔隙水压力处理方法把所给的地下水面作为测压管面把图11.17所给的流网转化为孔压网来计算根据渗流分析计算更精确的孔隙压力除外各家根据图11.17所示流网计算孔隙水压再用孔压计算安全系数表11.15分别示出这三种处理方法Baker和Donald的计算成果
STAB仅作了第一种情况如图11.18示各家提供的安全系数在1.46~1.53之间
图 11,17 考核题5渗流场
图 11,18 STAB程序对考核题EX5计算的临界滑裂面和最小安全系数
ACADS组织的这次考核的特点是详细地规定了计算原始数据和要求使不同的程序获得的成果具有可比性除第4题外其它考核题的设计都比较严密因此有条件使不同的程序算出基本相同的结果但由组织者提供的2(a)和3(b)的推荐的裁判答案与Fredlund
和陈祖煜提供的答案相差较大其权威性有待进一步论证综合参加的各程序提供的安全系数统计参数成果表11.13可见裁判答案互相比较接近而各单位自编程序提供的答案则比较分散早期Fredlund也曾经发现过这一问题这一事实说明在进行稳定分析时最好还是使用专门开发的经过考核的比较成熟的程序
边坡稳定分析是岩土工程设计和施工中经常要进行的工作在计算机迅速普及的今日
376 土质边坡稳定分析?原理
方法
程序
手算方式显然要被淘汰但确认各家自编程序的准确性确实成为了进入计算机时代面临的一个重要的问题已提到了议事日程这次程序考核工作做了一次有益的尝试本项工作为各单位自编的稳定分析程序的考核提供了十分有用的资料
11,5 本章附录
11,5,1 本章数据文件表11.6
表 11,16 本 章 数 据 文 件
有关章节 系列号 数据文件名 内容
11-01-01 CIRCLE.FOR 构筑圆弧滑裂面源程序
11-01-02 CIRCLE.DAT 构筑圆弧滑裂面例题
11.3.2
11-01-03 CIRCLE.C 构筑圆弧滑裂面例题计算结果
11-02-01 SPLINE.FOR 构筑任意形状滑裂面源程序
11-02-02 SPLINE-1.DAT 构筑光滑形状滑裂面例题
11-02-03 SPLINE-1.C 构筑光滑形状滑裂面例题计算结果
11-02-04 SPLINE-2.DAT 构筑组合形状滑裂面例题
11.3.3
11-02-05 SPLINE-2.C 构筑组合形状滑裂面例题计算结果
11-03-01 PROFILE.FOR 计算土条重量等特性数据的源程序
11-03-02 PROFILE.DAT 计算土条重量等特性数据的例题
11.3.4
11-03-03 PROFILE.C 计算土条重量等数据的例题计算结果
11-04-01 S.FOR 计算安全系数源程序
11-04-02 S0.DAT 计算通用条分法安全系数例题
11-04-03 S0.C 计算通用条分法安全系数例题计算结果
11-04-04 S1.DAT 计算Bishop法安全系数例题
11-04-05 S1.C 计算Bishop法安全系数例题计算结果
11-04-06 S2.DAT 计算瑞典法安全系数例题
11-04-07 S2.C 计算瑞典法法安全系数例题计算结果
11-04-08 S3.DAT 计算工程师团法安全系数例题
11-04-09 S3.C 计算工程师团法安全系数例题计算结果
11-04-10 S4.DAT 计算简化法安全系数例题
11.3.5
11-04-11 S4.C 计算简化法法安全系数例题计算结果
11.4.2 11-05-01 exp1a.dat ACADS考核题1(a)
11-05-02 exp1b.dat ACADS考核题1(b)
11-05-03 exp1c.dat ACADS考核题1(c)
11-05-04 exp1d.dat ACADS考核题1(d)
11-06-01 Exp2a.dat ACADS考核题2整体极值
11-06-02 Exp2b.dat ACADS考核题2心墙局部极值
11-06-03 Exp2c.dat ACADS考核题2指定滑裂面
11-07-01 Exp3a.dat ACADS考核题3整体极值
11-07-02 Exp3b.dat ACADS考核题3心墙局部极值
11-08-01 Exp4.dat ACADS考核题4
11-09-01 Exp5.dat ACADS考核题5
第11章 程序设计 377
参考文献
1 刘德贵费景高于泳江李广元,Fortran算法汇编第一分册,北京国防工业出版社
1980
2 Donald,I,B,And Giam,P,The ACADS slope stability programs review,Proc,6
Th
International Symposium
on Landslides,1992,Vol,3,1665-1670
3 Lumsdaine,R,W,and Tang,K,Y,A comparison of slope stability calculations1982,Proc,7
th
South East Asia
Geotechnical Conference,Hong Kong,31-38
4 Fredlund,D,G,Usage,requirements and features of slope stability computer software,Canadian Geotechnical
Journal,1987,15:83-95
5 ACADS Publication No,U255,Soil slope stability programs review,April,1989
