PFC模拟砂土滑坡

1年前浏览5028

滑坡、隧道、基坑算是岩土的三大问题了，这篇文章叙述一下砂土滑坡离散元模拟的标准过程。这里也是我自己的经验总结，不一定和传统的做法完全一致。

先说一下砂土滑坡的模拟步骤为：成样、预压、自重、切坡

一、成样

这部分也是经常做的事情，不多讲，注意这里的摩擦系数最好不要加，原因在后面说。


new def par    width=2.0    height=width*0.5     rdmax=0.009    rdmin=0.006     poro=0.08end@pardomain extent [-height*5] [height*5]set random 10001wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5ball distribute porosity @poro radius @rdmin @rdmax box [-width*0.5] [width*0.5] ...                                                        [-height*0.5] [height*0.5]cmat default model linear method deformability emod 100e6 kratio 1.5
ball attribute density 2.7e3 damp 0.7cycle 2000 calm 50solve save sample

二、预压

关于在边值问题（非单元试验）中是否加入预压过程，我也犹豫过，其实采用砂雨法成样应该才可以真实的模拟现实的过程。但是毕竟是模拟，这里采用预压也是为了释放试样中过大的内应力，这样在自重的时候，试样不至于飞出。



restore sample[txx=-1e6][tyy=-1e6][sevro_factor=0.5][do_xSevro=true][do_ySevro=true]
[sevro_freq=100][timestepNow=global.step-1]def sevro_walls    compute_stress   if timestepNow<global.step then        get_g(sevro_factor)        timestepNow =sevro_freq    endif    if do_xSevro=true then        Xvel=gx*(wxss-txx)        wall.vel.x(wpRight)=-Xvel        wall.vel.x(wpLeft)=Xvel    endif     if do_ySevro=true then        Yvel=gy*(wyss-tyy)        wall.vel.y(wpUp)=-Yvel        wall.vel.y(wpDown)=Yvel    endifenddef wp_ini    wpDown=wall.find(1)    wpRight=wall.find(2)    wpUp=wall.find(3)    wpLeft=wall.find(4)end@wp_inidef computer_chiCun    wlx=wall.pos.x(wpRight)-wall.pos.x(wpLeft)    wly=wall.pos.y(wpUp)-wall.pos.y(wpDown)enddef compute_stress    computer_chiCun    wxss=-(wall.force.contact.x(wpRight)-wall.force.contact.x(wpLeft))*0.5/wly    wyss=-(wall.force.contact.y(wpUp)-wall.force.contact.y(wpDown))*0.5/wlxenddef get_g(fac)    gx=0    gy=0    zongKNX=100e6*2*10    zongKNY=100e6*2*10    loop foreach ct wall.contactmap(wpLeft)        zongKNX =contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpRight)        zongKNX =contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpUp)        zongKNY =contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpDown)        zongKNY =contact.prop(ct,"kn")    endloop    gx=fac*wly/(zongKNX*global.timestep)    gy=fac*wlx/(zongKNY*global.timestep) endset fish callback -1.0 @sevro_wallshistory id 1 @wxsshistory id 2 @wysshistory id 3 @gxhistory id 4 @gycycle 1solve
save yuya

三、加自重

这里对应离心机原理，将重力加速度放大n倍，对应现实80m宽40m高的模型。

这里土体参数对应单元试验，也是在预压后加入，我这里没有变形参数，如果改变了变形参数的话，建议自重前再平衡一下。



restore yuyawall delete walls range id 3ball property fric 0.5set gravity [9.8*80/wlx]set fish callback -1.0 remove @sevro_wallscycle 1solvesave zizhong

平衡后，力链图如下所示：

四、切坡

这里采用删除颗粒模拟边坡开挖的过程，这里定义了边坡底部的位置：kengjiaoX与kengjiaoY。和边坡的角度 pojiao。这里设定了截止时间是算10s，这里也可以采用前文说的方法定时保存sav文件。



restore zizhongball attribute displacement multiply 0set mech age 0[kengjiaoX=wlx*0.5*0.2][kengjiaoY=0][pojiao=70]ball delete range plane origin @kengjiaoX @kengjiaoY dip @pojiao above plane origin @kengjiaoX @kengjiaoY dip 0  abovecycle 1solve time 10save qiepo