这个案例并没有做出理想的位移状态,可能和颗粒数或者别的什么东西影响,欢迎读者进行讨论。但是在土拱效应方面做得比较好。
首先讲一下模拟步骤:
成样 - 预压 - 加自重 - 染色 - 开挖
前面四个步骤前面都讲过了,这里讲一下隧道开挖的模拟。我们采用刚性墙wall来模拟隧道,wall不可以因为力和位移或者变形。所以这里就必须对地层损失的形态进行改变,如下图所示为常用的地层损失的模拟形态:
假定衬砌已经因为重力下沉,开挖面和衬砌相切。地层损失率的计算为:
[suidaobanjing=6.6*0.5*wlx/40.0]
[suidaoshendu=10*wlx/40.0]
[sunshilv=14e-2]
[waibuR=math.sqrt((1 sunshilv)*suidaobanjing^2)]
[g=waibuR-suidaobanjing]
需要自己定义隧道半径,隧道深度,地层损失率,然后自动计算后面的外部半径和g。这里为了摆脱颗粒数的影响,半径很大,深度比较小,地层损失率也比较大。
完整开挖代码如下,前面的部分就不放了,用户可以认真学一下边坡那个文章(PFC模拟砂土滑坡),注意这里模型尺寸为40*20。
restore ranse
ball attribute displacement multiply 0
def set_cedian
cedianfanwei=wlx*0.8
cediangeshu=21
loop n(1,cediangeshu)
pos_x=-cedianfanwei*0.5 n*cedianfanwei/float(cediangeshu-1)
pos_y=wly
bp=ball.near(pos_x,pos_y)
ball.group(bp,3)="cedian"
endloop
end
@set_cedian
[suidaobanjing=6.6*0.5*wlx/40.0]
[suidaoshendu=10*wlx/40.0]
[sunshilv=14e-2]
[waibuR=math.sqrt((1 sunshilv)*suidaobanjing^2)]
[g=waibuR-suidaobanjing]
def delete_ball
loop foreach bp ball.list
local dist=math.sqrt((ball.pos.x(bp))^2 (ball.pos.y(bp)-(wly*0.5-suidaoshendu g))^2)-ball.radius(bp)
if dist<waibuR then
ball.delete(bp)
endif
endloop
end
@delete_ball
wall generate circle position 0 [wly*0.5-suidaoshendu] radius [suidaobanjing]
cycle 1
solve
save kaiwa
计算完成后的位移图为:
可以看出由于砂土颗粒的拱效应,地层损失并没有对地面产生很大的影响,这里和peck曲线有点出入,也有可能是砂土的强度太大了,后面再继续研究。
拱效应无疑是隧道开挖研究的一个重点了,通过力链图可以非常直观的看出来:
当然布置测量圆绘制十字架效果也很好:
十字架绘制代码为:
restore kaiwa
call measure_tools
@create_measure([-wlx*0.5*0.5],[wlx*0.5*0.5],[wly*0.5*0.35],[wly*0.5*0.9],10,5,[wlx*0.03])
@create_measure([-wlx*0.5*0.5],[-wlx*0.5*0.2],[-wly*0.5*0.3],[wly*0.5*0.35],4,5,[wlx*0.03])
@create_measure([wlx*0.5*0.2],[wlx*0.5*0.5],[-wly*0.5*0.3],[wly*0.5*0.35],4,5,[wlx*0.03])
@outtensor
测量圆布置为:
十字架图:
应力十字架的旋转还是显而易见的。
下面给出更新后的measure_tools文件:
def create_measure(x_min,x_max,y_min,y_max,n_x,n_y,rad)
if n_x=1 then
split_x=0
else
split_x=(x_max-x_min-2*rad)/float(n_x-1)
endif
if n_y=1 then
split_y=0
else
split_y=(y_max-y_min-2*rad)/float(n_y-1)
endif
loop local n(1,n_x)
x_pos=x_min rad split_x*(n-1)
loop local m(1,n_y)
y_pos=y_min rad split_y*(m-1)
command
measure create position [x_pos] [y_pos] radius [rad]
endcommand
endloop
endloop
end
def get_poro_x
tb_poro=table.create("poro")
loop foreach mp measure.list
table(tb_poro,measure.pos.x(mp))=measure.porosity(mp)
endloop
end
def get_poro_y
tb_poro=table.create("poro")
loop foreach mp measure.list
table(tb_poro,measure.pos.y(mp))=measure.porosity(mp)
endloop
end
def get_poro_xy
tb_poro=table.create("poro")
tb_x_pos=table.create("x_pos")
tb_y_pos=table.create("y_pos")
loop foreach mp measure.list
table(tb_poro,measure.id(mp))=measure.porosity(mp)
endloop
loop foreach mp measure.list
table(tb_x_pos,measure.id(mp))=measure.pos.x(mp)
endloop
loop foreach mp measure.list
table(tb_y_pos,measure.id(mp))=measure.pos.y(mp)
endloop
end
def outtensor
array yinglitensor(10000)
weizhi_count=1
loop foreach local mplocal measure.list
p_sc = user.tensor.create(measure.pos.x(mplocal),measure.pos.y(mplocal))
user.tensor.value(p_sc)=measur.stress(mplocal)
yinglitensor(weizhi_count)=p_sc
weizhi_count =1
endloop
end
def outyinglibi(a,b)
array yinglibiarr(10000)
weizhi_count=1
loop foreach local mplocal measure.list
soil_stress_xx=measure.stress.xx(mplocal)
soil_stress_yy=measure.stress.yy(mplocal)
soil_stress_xy=measure.stress.xz(mplocal)
p=math.abs(soil_stress_xx soil_stress_yy)*0.5
q=math.sqrt((soil_stress_xx-p)^2 soil_stress_xy^2)
qf=a*p b
n=q/qf
p_sc = user.scalar.create(measure.pos.x(mplocal),measure.pos.y(mplocal))
user.scalar.value(p_sc)=n
yinglibiarr(weizhi_count)=p_sc
weizhi_count =1
endloop
end