PFC接触分组赋值
今天将接触分组这部分思考了一下，总结了一下接触分组的做法，并提出一个做法可以实现接触层面赋值的方法。
首先我们这里提出一个工况，有两层土，底层是砾石，上层是砂土，砾石的参数为：kn=ks=5e7，砂土为：kn=ks=2e7，组间接触为：kn=ks=3e7，颗粒与墙的接触定义为：kn=ks=1e8。
下面的demo是我经常用的做法，在cmat entry中将除了组间接触之外的所有接触都定义好，id 1对应的是砾石内部，2对应的是砂土内部，3对应的是颗粒与墙的。那组间接触在cmat entry中找不到自己的range时，就会返回到default中。
new
domain extent -1 1
wall generate box -0.5 0.5
ball distribute porosity 0.2 radius 0.02 0.03 group lishi box -0.5 0.5 -0.5 0
ball distribute porosity 0.2 radius 0.006 0.01 group shatu box -0.5 0.5 0 0.5
ball attribute density 3.8e3 damp 0.7 range group lishi
ball attribute density 2.6e3 damp 0.7 range group shatu
contact groupbehavior and
cmat default model linear property kn 3e7 ks 3e7
cmat add 1 model linear property kn 5e7 ks 5e7 range group lishi
cmat add 2 model linear property kn 2e7 ks 2e7 range group shatu
cmat add 3 model linear property kn 1e8 ks 1e8 range contact type ball-facet
cycle 2000 calm 10
set gravity 9.8
solve
运行的结果为：
这个是达到我们预期的结果的，但是我们也会发现这个问题的不足，当出现多个组时，这个方法就不太适用了。
受到裂纹函数的启发，这里利用contact_activated事件进行接触的定义，demo如下：
new
domain extent -1 1
wall generate box -0.5 0.5
ball distribute porosity 0.2 radius 0.02 0.03 group lishi box -0.5 0.5 -0.5 0
ball distribute porosity 0.2 radius 0.006 0.01 group shatu box -0.5 0.5 0 0.5
ball attribute density 3.8e3 damp 0.7 range group lishi
ball attribute density 2.6e3 damp 0.7 range group shatu
contact groupbehavior and
cmat default model linear property kn 1e9 ks 1e9
cmat add 1 model linear property kn 5e7 ks 5e7 range group lishi
cmat add 2 model linear property kn 2e7 ks 2e7 range group shatu
cmat add 3 model linear property kn 1e8 ks 1e8 range contact type ball-facet
;
def set_contact(entry)
local ct=entry(1)
if type.pointer.name(contact.end1(ct)) # "Facet" then
if type.pointer.name(contact.end2(ct)) # "Facet" then
if ball.group(contact.end1(ct))="lishi" then
if ball.group(contact.end2(ct))="shatu" then
contact.model(ct)="linear"
contact.prop(ct,"kn")=3e7
contact.prop(ct,"ks")=3e7
endif
endif
if ball.group(contact.end2(ct))="lishi" then
if ball.group(contact.end1(ct))="shatu" then
contact.model(ct)="linear"
contact.prop(ct,"kn")=3e7
contact.prop(ct,"ks")=3e7
endif
endif
endif
endif
end
set fish callback contact_activated @set_contact
cycle 2000 calm 10
set gravity 9.8
solve
这个是用set_contact函数进行接触的定义，在其中筛选出组间的接触，然后对其进行定义，为了区分之前的，这里将cmat default中的参数设置为1e9，我们可以看一下运行结果：
这个也是符合我们需求的。
下面测试一下三个分组的：
这里增加了底层为块石，定义其接触参数为：kn=ks=8e7，砾石和块石间的接触参数为：kn=ks=7e7。
将demo修改为：
new
domain extent -1 1
wall generate box -0.5 0.5
ball distribute porosity 0.2 radius 0.03 0.05 group kuaishi box -0.5 0.5 -0.5 -0.3
ball distribute porosity 0.2 radius 0.02 0.03 group lishi box -0.5 0.5 -0.3 0.1
ball distribute porosity 0.2 radius 0.006 0.01 group shatu box -0.5 0.5 0.1 0.5
ball attribute density 3.8e3 damp 0.7 range group kuaishi
ball attribute density 3.8e3 damp 0.7 range group lishi
ball attribute density 2.6e3 damp 0.7 range group shatu
contact groupbehavior and
cmat default model linear property kn 1e9 ks 1e9
cmat add 1 model linear property kn 5e7 ks 5e7 range group lishi
cmat add 2 model linear property kn 2e7 ks 2e7 range group shatu
cmat add 3 model linear property kn 1e8 ks 1e8 range contact type ball-facet
cmat add 1 model linear property kn 8e7 ks 8e7 range group kuaishi
;
def set_contact(entry)
local ct=entry(1)
if type.pointer.name(contact.end1(ct)) # "Facet" then
if type.pointer.name(contact.end2(ct)) # "Facet" then
if ball.group(contact.end1(ct))="lishi" then
if ball.group(contact.end2(ct))="shatu" then
contact.model(ct)="linear"
contact.prop(ct,"kn")=3e7
contact.prop(ct,"ks")=3e7
endif
endif
if ball.group(contact.end2(ct))="lishi" then
if ball.group(contact.end1(ct))="shatu" then
contact.model(ct)="linear"
contact.prop(ct,"kn")=3e7
contact.prop(ct,"ks")=3e7
endif
endif
if ball.group(contact.end1(ct))="lishi" then
if ball.group(contact.end2(ct))="kuaishi" then
contact.model(ct)="linear"
contact.prop(ct,"kn")=3e7
contact.prop(ct,"ks")=3e7
endif
endif
if ball.group(contact.end2(ct))="lishi" then
if ball.group(contact.end1(ct))="kuaishi" then
contact.model(ct)="linear"
contact.prop(ct,"kn")=7e7
contact.prop(ct,"ks")=7e7
endif
endif
endif
endif
end
set fish callback contact_activated @set_contact
cycle 2000 calm 10
set gravity 9.8
solve
运行的结果为：
可以看到也是达到预期了。
这个方法有个缺点就是速度会稍微慢点，如果两个分组的话，采用第一个方法好点。
查看更多评论 >