PFC-FLAC耦合技术在解决这种大尺度变形上就有很高的价值了,将非大变形区域的颗粒用zone来模拟,而我们将颗粒数花费在大变形破坏区域,这样就可以提高模拟的价值了。
一、二维的耦合测试
下方是耦合的测试代码,FLAC2D的关键词并没有在PFC的手册中,各位可以参考ITASCA官网的关键词手册(Itasca Software Guide — Itasca Software 9.0 documentation (itascacg.com))。这里我们用quad关键词生成一个矩形的地基,地基首先进行地应力沉降。注意边界设置,底部固定y向,左右固定x向即可。
沉降完成后,生成一个密度较高的颗粒落在地基上,非常简单的案例,测试足够。
model new
domain extent -2 2
zone create quad size 50 5 point 0 0 0 ...
point 1 1 0 ...
point 2 0 0.2
zone gridpoint fix velocity-y range pos-y -0.001 0.001
zone gridpoint fix velocity-x range pos-x -0.001 0.001
zone gridpoint fix velocity-x range pos-x [1-0.001] [1+0.001]
zone initialize density 3e3
zone cmodel assign elastic
zone property young 10e6 poisson 0.3
model gravity 9.8
model cycle 1
model solve
model save "diji"
model restore "diji"
wall-zone create
ball create position 0.5 0.4 radius 0.02
contact cmat default model linear property kn 1e9 ks 1e9 fric 0.5
ball attribute density 3e6 damp 0.7
model cycle 2000
可以看到耦合效果还是可以的,在落点位置出现较大的位移。
这个测试完成后,就尝试进行一个略复杂的砂土边坡模拟了。
二、颗粒部分生成
我们先生成边坡的颗粒部分,这也是我一般的做法,离散-有限耦合的时候,先完成离散部分的成样,然后再扩展有限元部分的元素。我们将坡脚定义在原点位置,然后定义好边坡的宽和高。里面还有个数值控制着离散元左边界与边坡左边界的距离 width*0.1 ,这个数值越大,离散的部分就越小,这个数值在设置的时候应当充分考虑边坡的大变形范围。
还有一点是distribute成样是先在domain范围内生成颗粒,然后删除range以外的部分,所以domain应当足够小。
model new
def par
width=2.0
height=width*0.5
rdmin=0.001
rdmax=0.002
poro=0.08
end
@par
model domain extent [-width*0.5] [width*0.1] [-width*0.1] [height*0.6]
wall generate box [-width*0.5+width*0.1] [0] [0] [height*0.5]
ball distribute porosity @poro radius @rdmin @rdmax box ...
[0] [0] [height*0.5]
ball attribute density 3e3 damp 0.7
cmat default model linear method deform emod 10e6 kratio 1.5 property fric 0.5
model cycle 2000 calm 20
model solve
model save "sample_keli"
三、有限元部分生成
离散有限部分耦合的关键点是要保证材料的一致性,这个其实难度很大,对参数需要足够多的时间去调整。对于这种考虑重力的工况而言,除了常规的力学参数,主要是变形参数保持基本一致外,密度是需要调整的。PFC的密度是给到的ball,同样的式样体积下,zone部分的固体体积和式样体积一致,但是ball部分由于孔隙率和重叠量的存在,总的质量是有浮动的,所以需要调整密度,使得zone部分的地应力或者变形场基本水平即可。
model restore "sample_keli"
model mechanical time-total 0
ball attribute displacement multiply 0
geometry set "zone_range"
ball attribute density 2.8e3 damp 0.7
zone create quad size 2 10 point 0 [-width*0.5] [-height*0.5] ...
point 1 [-width*0.5+width*0.1] [-height*0.5] ...
point 2 [-width*0.5] [height*0.5]
zone create quad size 10 5 point 0 [-width*0.5+width*0.1] [-height*0.5] ...
point 1 [0] [-height*0.5] ...
point 2 [-width*0.5+width*0.1] [0]
zone create quad size 10 5 point 0 [0] [-height*0.5] ...
point 1 [width*0.5] [-height*0.5] ...
point 2 [0] [0]
zone create quad size 10 5 point 0 [0] [0] ...
point 1 [width*0.5] [0] ...
point 2 [0] [height*0.5]
zone initialize density 3e3
zone cmodel assign elastic
zone property young 10e6 poisson 0.3
cmat default type ball-facet model linear method deform emod 10e6 kratio 1.5 property fric 0 lin_mode 1
cmat default type ball-ball model linear method deform emod 10e6 kratio 1.5 property fric 0 lin_mode 1
cmat apply
model clean
contact property lin_force 0 0
ball attribute force-contact 0 0
wall delete
wall create id 100 vertices [0] [0] [0] [height*0.6]
wall create id 101 vertices [-width*0.5+width*0.1] [0] [-width*0.5+width*0.1] [height*0.6]
wall-zone create range pos-x [-width*0.5+width*0.1-0.0001] [0.0001] ...
pos-y [-0.00001] [height*0.5+0.0001]
zone gridpoint fix velocity-y range pos-y [-height*0.5-0.0001] [-height*0.5+0.0001]
zone gridpoint fix velocity-x range pos-x [-width*0.5-0.0001] [-width*0.5+0.0001]
zone gridpoint fix velocity-x range pos-x [width*0.5-0.0001] [width*0.5+0.0001]
model gravity [9.8*40]
model cycle 1000
model solve ratio-average 1e-3
wall delete walls range id 100 101
model solve
model save "zizhong"
将色条设为一致,可以发现这里的密度还不是一个非常好的状态,颗粒底部的位移明显小了点,说明密度还需要增大。而颗粒内部的连续并不需要过于在意,由于颗粒洒落的时候并不是连续变形,会导致初始时候中间颗粒的位移较大。
而应力也会有相应的变化,一般而言并不会完完全全的呈现水平带状分布,调整参数到差不多的状态就可以了。
四、开挖
这里开挖主要删除两个部分,第一是右上方的zone,第二十根据坡度产生的斜半平面以上的颗粒。常规有限元的生死单元和设置模型是null的方法,在这里只需要删除颗粒就可以了,这个还是比较好的。
注意一下初始位移量清零,还有一点是ball-ball接触参数应和耦合部分的ball-facet参数保持一致。
model restore "zizhong"
model mechanical time-total 0
model domain extent [-width*0.5] [width*0.6] [-height*0.6] [height*0.6]
ball attribute displacement multiply 0
zone gridpoint initialize displacement 0 0
zone delete range pos-x 0 [width*0.5] pos-y 0 [height*0.5]
wall-zone create range pos-x [0] [width*0.6] ...
pos-y [-0.2] [2]
wall property "fric" 0.5
ball property "fric" 0.5
ball delete range plane origin 0 0 normal ...
[math.cos(math.pi/6.0)] [math.sin(math.pi/6.0)]
[final_time=2]
[baocunpinlv=final_time/20.0]
[time_record_savefile=-1]
[count=0]
def savefile
time=mech.time.total
if time-time_record_savefile >= baocunpinlv then
filename=string.build("jieguo%1",count)
command
model save @filename
endcommand
time_record_savefile=time
count +=1
endif
end
fish callback add @savefile -0.5
model cycle 1
model solve time @final_time
测试的时候把收敛时间设置为1s,各位可以根据计算情况继续往下面算。下面为1s后的y向位移场,在分析的时候应当将色卡设置一致,可以看到总的来说,变形的连续性还尚可,调整变形参数应当会有更好的结果。
做成动图效果也是非常好看的。
七、Lobby的PFC6.0入门学习讲义
在此向读者朋友推荐我今年新上线的新课程《Lobby的PFC6.0入门学习讲义》适合新手入门。
本课程为PFC6.0软件入门,常用的案例技术解析;全字幕的课程,由简入深的讲解,用户可以获得多达二十种案例,丰富的技术点讲解,完全从零开始的案例实操,覆盖clump/cluster/rblock/离散连续耦合/流固耦合/柔性膜等内容,让用户掌握基本的后处理技能。
以下是我的课程安排
扫码立即试看
(完)