工程地质学程序.docx
《工程地质学程序.docx》由会员分享，可在线阅读，更多相关《工程地质学程序.docx（24页珍藏版）》请在冰豆网上搜索。
工程地质学程序
READ（5,*）NS!
土条总数
READ（5,*）（X2（I）,I=1,NS）
READ（5,*）（Y2（I）,I=1,NS）
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）
302CONTINUE
!
读入两端点和土条底中点的数据,其中涉及力的物理量均以单宽土条计
!
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
709FORMAT（1X,'水容重=',F10.3//）
WRITE（6,702）
702FORMAT（T25,'DATAFORASSUMEDSIDE
READ（5,*）NS!
土条总数
READ（5,*）（X2（I）,I=1,NS）
READ（5,*）（Y2（I）,I=1,NS）
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）
302CONTINUE
!
读入两端点和土条底中点的数据,其中涉及力的物理量均以单宽土条计
!
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
709FORMAT（1X,'水容重=',F10.3//）
WRITE（6,702）DO208I=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）））
208CONTINUE
ALF
（1）=ALF
（2）
ALF（N）=0
X（N）=X2（NS）
Y（N）=Y2（NS）
X
（1）=X2
（1）
Y
（1）=Y2
（1）
WRITE（6,716）
716FORMAT（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）
901CONTINUE
718FORMAT（1X,I5,11F10.3）
WRITE（6,719）
719FORMAT（1X,/'NOTE:
'）
WRITE（6,720）
720FORMAT（5X,'W=土条的单宽重量'）
WRITE（6,721）
721FORMAT（5X,'RU=孔隙水压力系数'）
WRITE（6,722）
722FORMAT（5X,'C=凝聚力'）
WRITE（6,723）
723FORMAT（5X,'F=摩擦系数'）
WRITE（6,724）
724FORMAT（5X,'BDA=土条的平均容重'）
WRITE（6,725）
FORMAT（5X,'WL=浸润线的Y坐标'）
SELECTCASE（OPTION）
CASE（0）
WRITE（6,701）
701FORMAT（/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
301CONTINUE!
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）
704FORMAT（/20X,'BISHOP法'/）
CALLSAFE2（X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION）
CASE
（2）
WRITE（6,705）
705FORMAT（/20X,'瑞典法'/）
CALLSAFE2（X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION）
CASE（3）
WRITE（6,731）
731FORMAT（/10X,'工程师团法'）
CALLSAFE2（X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,4）
IWALL=0
DO311I=1,80
FXO（I）=0
FX（I）=1
311CONTINUE!
tan（BETA（I））=FXO（I）+ALAM*FX（I）
KXYX=1
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（4）
WRITE（6,789）
789FORMAT（/20X,'简化法法'/）
CALLSAFE2（X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION）
CASEDEFAULT
ENDSELECT
END
SUBROUTINESAFE2（X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,ID）
REALX2（80）,F（80）,FQUH（80）,W（80）,RU（80）,ALF（80）,X（80）,C（80）
REALAB1,FDIS（80）
CINTEGERN
REALRO1,DM,RM,BX1,CDN,DG（80）,G3,DX1,G,GX
REALO,S1,FRIC,S2,S,DRM,P1,P2,DDM,DBX1,DCDN
CALLSAFE2（X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION）
CASE
（2）
WRITE（6,705）
705FORMAT（/20X,'瑞典法'/）
CALLSAFE2（X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION）
CASE（3）
WRITE（6,731）
731FORMAT（/10X,'工程师团法'）
CALLSAFE2（X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,4）
IWALL=0
DO311I=1,80
FXO（I）=0
FX（I）=1
311CONTINUE!
tan（BETA（I））=FXO（I）+ALAM*FX（I）
KXYX=1
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（4）
WRITE（6,789）
789FORMAT（/20X,'简化法法'/）
CALLSAFE2（X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION）
CASEDEFAULT
ENDSELECT
END
SUBROUTINESAFE2（X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,ID）
REALX2（80）,F（80）,FQUH（80）,W（80）,RU（80）,ALF（80）,X（80）,C（80）
REALAB1,FDIS（80）
CINTEGERN
REALRO1,DM,RM,BX1,CDN,DG（80）,G3,DX1,G,GX
REALO,S1,FRIC,S2,S,DRM,P1,P2,DDM,DBX1,DCDN
C79CONTINUE
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
C79CONTINUE
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
202CONTINUE
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
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））
21CONTINUE
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
740FORMAT（1X,I5,F10.5）
IF（ABS（AB2-AB1）.LT.0.0001）GOTO49
IF（LM.LT.10）GOTO9191
WRITE（6,9192）
9192FORMAT（1X,'ITERATIONTROUBLESINBISHOPMETHOD'）
GOTO49
9191CK=1-DT/DM!
牛顿迭代法
AB1=AB1-（AB1-AB2）/CK
LM=LM+1
GOTO41
49WRITE（6,203）AB1
203FORMAT（/5X,'BISHOP法安全系数=',F10.4）
RETURN
END
SUBROUTINEDETE（AB1,RW,DQ,DM,ALAM,KS,BN,FX,FXO,
&C,F,X2,Y2,W,RU,FQUH,YSF,WL,FDIS,N,DB,KXYX,DF）
REALAB1,DQ,DT,DM,ALAM,DF
REALC（80）,F（80）,X2（80）,Y2（80）,W（80）,RU（80）,FXO（80）,FX（80）
REALFQUH（80）,YSF（80）,WL（80）,FDIS（80）,BN（80）
INTEGERKS,N,KXYX
REALX（80）,ALF（80）,G,Y（80）
REALEXT,EXR,EYR,EMT,EMR,FMT,DMT,DEG,DX1
REALDMG,BS,R1X,TX,QX,PT112,R1XT,POX
REALPOXT,POY,DPT112,DPX112,BSG,H3,YT（80）,GF（80）
REALGFX（80）,GFY（80）,BET（80）,EFX（80）
REALG1,TA,TA1,OP,FI,FI1,CE,OQ,OQ1,THO,CW
REALSEGB,SEGB1,FI0,AXS11,OTEM,DDQ,PX112,DBS
REALSX,SXG,DTX1,DTX2,DTX,TX1,DDP,DBS1,DBS2,EPX112
REALAAKX1,DR1X1,DR1X2,DR1X,DR1XT
REALDG,DDT,GAVE,BETA,DGX,OX,OY,HSL
REALTAA,OQS2,OQS1,DG1,BT（80）,EXR1,EXT1,EYR1
REALEYT1,XCE,YCE,XRE,YRE,DR2X11,DR2X12,DR2X1
REALGS,DPOY1,DPOY2,OY1,DDM,DDMT,PM,PMT,DDEG,DDMG
REALOQ0,SEGB0,DGG,DMG1,DEG1,OXX
REALTEMP1,TEMP2,GG,DDM2,DM1,EWALL,GWALL,PX
REALDQM,DMM,DTM,DMTM,DEGM,DMGM,S1,HMW,DDM1,ETA
REALGTENSION,HITE,WP,RW
INTEGEREYT,IWALL,NS
!
本程序计算Morgenstern-Price,Spencer和工程师团法
IWALL=0.!
土压力的作用位置
HMW=0.
GTENSION=0.
HITE=0.
GWALL=0.!
土压力的大小
!
GWALL=1计算土压力设定条件下的安全系数
!
GWALL=2计算主动土压力
IF（KS.EQ.0）THEN!
KS调用本程序段标志符,=0表示第1次调用
WRITE（6,705）
705FORMAT（1X,/'VALUEOFFX0（I）INORDEROFI'）
WRITE（6,706）（FXO（I）,I=1,N）
706FORMAT（1X,5F10.5）
WRITE（6,707）
707FORMAT（1X,/'VALUEOFFX（I）INORDEROFI'）
WRITE（6,708）（FX（I）,I=1,N）
708FORMAT（1X,5F10.5）
WRITE（6,726）
726FORMAT（1X,/'迭代过程:
'/）
WRITE（6,712）
712FORMAT（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）））
303CONTINUE
ALF
（1）=ALF
（2）
ALF（N）=0.
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
304CONTINUE
401CONTINUE
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
47
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）
307DBS2=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））
309DTX2=（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））
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）
310DR1X2=（-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
48
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）*C
升级会员 