引入Flac计算结果-浅析PFC6.0降雨入渗模拟框架（赠项目包）

8小时前浏览11

导读：近日，某粉丝私信我：“你仿真秀官网的PFC5.0教程即将实现500万播放量，PFC6.0教程也已经突破100万播放量”，这确实让人惊讶！

近期我正在写一些新的PFC6.0案例，计划策划全新的PFC视频课程，大家敬请期待。

PFC本身模拟水的能力并不是很强大，通常依赖于第三方的计算软件进行流体计算后进行数据交互。降雨入渗算是在岩土行业里面很广泛的一个验算场景了，而PFC在此方面还略有不足。本文参考高品红等于今年5月份在《岩石力学》期刊发表的《降雨作用下花岗岩残积土边坡模型试验及离散元分析》一文。利用flac3d的渗流模块计算降雨入渗下的饱和度，然后根据饱和度更新接触上的强度参数，达到失稳计算的状态。

本文模拟步骤是：

step-1：首先生成离散元边坡模型；

step-2：生成与之相对应的有限元模型，设置到流体边界和参数；

step-3：开始流体计算，并获取饱和度，更新到颗粒和接触参数上。

一、成样

这里先生成指定尺寸的立方体式样，注意高度需要设置的略高一点，比如我这里目标高度其实是80，但是设置了120，方便后面做地基削平。

model newdef BianpoPar    width=100    length=30    height=120    rdmin=0.8    rdmax=1.0end@BianpoParmodel domain extent -1 [width+1] -1 [length+1] -1 [height+1]wall generate box 0 [width] 0 [length] 0 [height]ball distribute radius @rdmin @rdmax porosity 0.3  box 0 [width] ...    0 [length] 0 [height]ball attribute density 3e3 damp 0.7contact cmat default type ball-facet model linear method deform emod 1e7 kratio 1.5 property fric 0.3model cycle 2000 calm 50ball delete range position-x 0 [width] notball delete range position-y 0 [length] notball delete range position-z 0 [height] notmodel solve model save "sample"

二、自重

在地应力下自重平衡

model restore "sample"model gravity 9.8model cycle 10model solvemodel save "zizhong"

三、地基削平

这里进行4次地基表面削平的操作

model restore "zizhong"def diji    loop n(1,4)        command            ball delete range position-z 80 1000            model cycle 1            model solve        endcommand    endloopend@dijimodel save "diji"

四、加参数

胶结在这一步施加结束

model restore "diji"contact cmat default type ball-facet model linear method deform emod 1e7 kratio 1.5 property fric 0contact cmat default type ball-ball model linearpbond method deform emod 1e7 kratio 1.5 ...    pb_deform emod 1e7 kratio 1.5 property ...    pb_ten 1e6 pb_coh 1e6 pb_fa 30 fric 0.3cmat applymodel cycle 1model solvecontact method bond gap [rdmin]model cycle 1model solvemodel save "canshu"

五、削坡

到这一步为止，实际上离散元的边坡我们已经做好了。

model restore  "canshu"ball delete range  plane origin [20] 0 80 dip-direction 90 dip [math.atan2(60.0,40.0)*180/math.pi] above ...        position-z  [20.0] 100model cycle 1model solveball delete range  plane origin [20] 0 80 dip-direction 90 dip [math.atan2(60.0,40.0)*180/math.pi] above ...        position-z  [20.0] 100model cycle 1model solvemodel save "xuepo"

六、有限元网格生成

这里引入流体计算，边坡用两个brick和一个uniform-wedge组成。webge是三棱柱的建模方法，加不加uniform是两种网格划分方式，用户可以自己比对一下。

这里加了一个scale参数，方便对网格进行细分。基于本工况的平面模型特征，y向不做细分。

然后对入渗边界上的zone划入分组。下面就是设定一些常规参数，比如渗透率，给到常规黏土的数值，还有液体的体积模量，土体孔隙率。这里设定渗透边界的流体通量为1.39e-3 m3/s-1，是一般大降雨的数值。我们这里不关心zone的力学计算，力学计算虽然开着，但是把zone的节点全部固定住。

model restore "xuepo"model configure fluid[sclae=2][singleHoudu=4/float(sclae)]zone create brick size [40*sclae] 1 [5*sclae]  point 0 (0,0,0) point 1 (100.0,0,0) point 2 (0,30,0) point 3 (0,0,20)zone create brick size [8*sclae] 1 [15*sclae] point 0 (0,0,20) point 1 (20.0,0,20) point 2 (0,30,20) point 3 (0,0,80)zone create uniform-wedge size [16*sclae] 1 [15*sclae] point 0 (20,0,20) point 1 (60,0,20) point 2 (20,30,20) point 3 (20,0,80)zone group "None"zone  group "upBianjie" range position-z  [80-singleHoudu] 80.0 zone  group "upBianjie" range position-x 60 100 position-z  [20-singleHoudu] 20.0 zone  group "upBianjie" range plane origin [20-singleHoudu] 0 80 dip-direction 90 dip [math.atan2(60.0,40.0)*180/math.pi] above ...        position-z  [20.0-singleHoudu] 100zone fluid cmodel assign isotropiczone fluid property permeability 3.22e-7  porosity 0.3zone fluid biot offzone gridpoint initialize fluid-modulus 2.15e2zone gridpoint initialize saturation 0zone  apply well 1.39e-3 range group "upBianjie"zone initialize density 3e3zone gridpoint fix velocity-x zone gridpoint fix velocity-yzone gridpoint fix velocity-zzone initialize fluid-density 1e3zone gridpoint initialize fluid-tension 0.0model mechanical active onmodel fluid active onmodel save "liuti"

七、渗流计算

下面就是需要耦合流体计算和离散元颗粒的部分。这里有两个update开头的函数，updateBallBaohe函数主要是对每一个ball，找到离它最近的一个gridpoint，获取上面的saturation数值，认为是当前颗粒的饱和度。第二个updateConatctPar是将胶结两端颗粒的饱和度的平均值，作为接触上的饱和度，并且认为饱和度越高，胶结半径越小。

注意，这里的饱和度和胶结半径的关系是我随意提的，具体衰减规律应当结合实际的物理实验得到，但这个不是阻塞当前技术的概念了。

再往下就是跳跃结构进行update,也进行文件的保存，不再赘述。

本工况工计算20分钟，每10秒保存一个sav文件

model restore "liuti"ball attribute displacement multiply 0zone cmodel elasticzone property young 10e6 poisson 0.3model mechanical time-total 0model fluid time-total 0def delete_float    loop foreach bp ball.list        ct_num=0        loop foreach ct ball.contactmap(bp)            if contact.model(ct)=="linearpbond" then                if contact.prop(ct,"pb_state")==3 then                    ct_num+=1                endif            endif        endloop        if ct_num==0 then            ball.delete(bp)        endif    endloopend@delete_floatdef updateBallBaohe    loop foreach bp ball.list        ball_pos=ball.pos(bp)        gp=gp.near(ball_pos)        s=gp.sat(gp)        ball.extra(bp,1)=s    endloopenddef updateContactPar    loop foreach ct contact.list        if contact.model(ct)#"linearpbond" then            continue        endif        s1=ball.extra(contact.end1(ct),1)        s2=ball.extra(contact.end2(ct),1)        contact.prop(ct,"pb_rmul")=1-(s1+s2)*0.5    endloopend[update_freq=1000][update_time_record=-1000]def updateLiutiCanshu    whilestepping    time_mech=mech.time.total    if global.step<update_time_record then        exit    endif    update_time_record=global.step+update_freq    updateBallBaohe    updateContactParend

[fluid_timerecord=-100][fluid_save_freq=10][count=0]def savefile      whilestepping    time=fluid.time.total    if time-fluid_timerecord >= fluid_save_freq then        filename=string.build("fluid_jieguo_%1",count)        command            model save @filename        endcommand         fluid_timerecord=time        count +=1    endif    end


model solve time [20*60]        
model save "result"