隧道开挖地层损失产生的土拱效应

1年前浏览2769

这个案例并没有做出理想的位移状态，可能和颗粒数或者别的什么东西影响，欢迎读者进行讨论。但是在土拱效应方面做得比较好。

首先讲一下模拟步骤：

成样 - 预压 - 加自重 - 染色 - 开挖

前面四个步骤前面都讲过了，这里讲一下隧道开挖的模拟。我们采用刚性墙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 0def 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"    endloopend@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    endloopend@delete_ball

wall generate circle position 0 [wly*0.5-suidaoshendu] radius [suidaobanjing]cycle 1solve 

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

测量圆布置为：

微信截图_20220718155956.png

十字架图：

应力十字架的旋转还是显而易见的。

下面给出更新后的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    endloopend





def get_poro_x    tb_poro=table.create("poro")    loop foreach mp measure.list        table(tb_poro,measure.pos.x(mp))=measure.porosity(mp)    endloopend

def get_poro_y    tb_poro=table.create("poro")    loop foreach mp measure.list        table(tb_poro,measure.pos.y(mp))=measure.porosity(mp)    endloopend

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)    endloopend

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    endloopend

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    endloopend