陈育明和徐鼎平 FLAC/FLAC3D 基础与工程实例中P267 有一个二分法强度折减法的例子。我有几个问题想咨询下，程序的前半部分都看懂了，看到这里就不明白了， 请有经验的同行指教。
fosfile0='_fos'+'.sav'
command
save fosfile0
endcommand
end
;程序执行及结果显示
SSR
pr ks
此外这个程序在FLAC3D 3.0 上能正常运行，可在FLAC3D 4.0 测试版中运行不了，请问需要做改动吗？请有经验的同行指教。程序如下，
1.
new
;=====================================
;建立网格模型
gen zone brick &
p0 0 0 0 p1 2 0 0 p2 0 0.5 0 p3 0 0 3 size 3 1 3
gen zone brick &
p0 2 0 0 p1 20 0 0 p2 2 0.5 0 p3 2 0 3 size 17 1 3 &
ratio 1.03 1 1
gen zone brick &
p0 2 0 3 p1 20 0 3 p2 2 0.5 3 p3 12 0 13 &
p4 20 0.5 3 p5 12 0.5 13 p6 20 0 13 &
p7 20 0.5 13 size 17 1 17 ratio 1.03 1 1
;**********************************************
;自定义强度折减法
def SSR
;=====================================
;定义有关参数及循环终止条件
ait1=0.02
k11=1.0
k12=2.0
ks=(k11+k12)/2
loop while (k12-k11)>ait1
coh1=12380/ks
fri1=(atan((tan(20*pi/180))/ks))*180/pi
dila1=20.0
ten1=1e6
grav0=-10
dens1=2000
K1=1e8
G1=3e7
;=====================================
;折减的实现过程
command
model null
;初始应力场的生成
model elastic
pro bulk 1e10 she 3e9 dens dens1
fix x y z range z -0.1 0.1
fix x range x 19.9 20.1
fix x range x -0.1 0.1
fix y
set grav 0 0 grav0
solve
ini xdisp 0 ydisp 0 zdisp 0
ini xvel 0 yvel 0 zvel 0
;塑性阶段求解
model mohr
pro bulk K1 she G1 dens dens1 coh coh1 &
friction fri1 dil dila1 tens ten1
set mech ratio 9.8e-6
solve step 10000
endcommand
;二分法的实现过程
if mech_ratio<1.0e-5
k11=ks
k12=k12
else
k12=ks
k11=k11
endif
ks=(k11+k12)/2
endloop
;=====================================
;计算结果的保存
fosfile0='_fos'+'.sav'
command
save fosfile0
endcommand
end
;**********************************************
;程序执行及结果显示
SSR
pr ks
