PFC教程_pb模型

孤亭远见

已于 2023-08-06 21:33:07 修改

阅读量6k

点赞数 16

分类专栏： PFC 文章标签：开发语言

于 2023-04-30 16:41:03 首次发布

本文链接：https://blog.youkuaiyun.com/weixin_52546971/article/details/130448410

版权

PFC 专栏收录该内容

30 篇文章

订阅专栏

pb模型

模型简介

linear contact bond model ：模拟胶结物和被胶结物分界不是很清楚

linear parallel bond model ：岩石（线性弹簧+）

控制变形:pb_kn、pb_ks

控制强度：pb_ten、pb_coh、pb_fa

胶结模型可以承受一定拉力。

图3：平行结合力和力矩的力-位移规律：（a）法向力与平行结合面间隙；（b）剪切力与相对剪切位移；（c）扭曲力矩与相对扭曲旋转；（d）弯曲力矩与相对弯曲旋转。

法向验证

[ct=contact.find("ball-ball",1)] 查找接触

pb模型在纯压状态永远不会破坏。

model.clean命令是一个全局命令，可以在模型创建过程中的任何时候给出，以创建模型中的所有接触，初始化所有件属性，确保数据结构存在以进行空间搜索，并根据几何形状更新接触活动状态。

pb_stage 在右侧plot窗口中选择查看破坏状态。

new 

;纯压，8s应力值为8Mpa
;纯拉，2s应力值为1Mpa
domain extent -100 100 

ball create position 0 -10 radius 10
ball create position 0 10 radius 10

ball fix vel spin range id 1
ball fix vel spin range id 2
ball attribute vel 0 1 range id 2 ；可以多次调试竖向速度为1或者-1，横向速度为1或者-1

ball attribute density 2.3e3 damp 0.7

;力fbond=pb_kn*disp*A---------bond（对于bp胶结模型）其中A表示面积
;力flinear=kn*disp-------linear（对于线性模型）

[pb_kn=1e7/20.0]
；pb_kn=kn/A=1e7/20.0
cmat default model linearpbond property ...
                    kn 0 pb_kn @pb_kn pb_ks @pb_kn ... 
                    pb_ten 1e6 pb_fa 45 pb_coh 1e6   
;pb_kn和pb_ks控制变形，pb_ten 1e6 pb_fa 45 pb_coh 1e6 控制强度
clean
contact method bond gap 1
set timestep fix 1e-3
[ct=contact.find("ball-ball",1)]
def jiance   ；进行监测
    whilestepping  ；进行回调
    faxiang_force=contact.force.normal(ct)  ；法向力
    faxiang_yingli=faxiang_force/20.0     ；法向应力
    
    qiexiang_force=contact.force.shear(ct)
    qiexiang_yingli=qiexiang_force/20.0
    
    qiexiang_fengzhi=1e6+math.tan(math.pi*0.25)*faxiang_yingli
；   计算切向力的峰值，来给出包络线c+σ*tanφ
;    lin_stress=comp.x(contact.prop(ct,"lin_force"))/20.0
;    pb_stress=contact.prop(ct,"pb_sigma")
end
history id 1 @faxiang_yingli
history id 2 @qiexiang_yingli
history id 3 @qiexiang_fengzhi
solve time 8
save faxiang1

全量计算与增量计算

全量：累积之前的记录

增量：不考虑之前的记录

线性部分可以有全量和增量，而胶结部分只有增量。

def jiance   ；进行监测
    whilestepping  ；进行回调
    faxiang_force=contact.force.normal(ct)  ；法向力
    faxiang_yingli=faxiang_force/20.0     ；法向应力
    
    qiexiang_force=contact.force.shear(ct)
    qiexiang_yingli=qiexiang_force/20.0
    
    qiexiang_fengzhi=1e6+math.tan(math.pi*0.25)*faxiang_yingli
    lin_stress=comp.x(contact.prop(ct,"lin_force"))/20.0
    pb_stress=contact.prop(ct,"pb_sigma")
    ；comp.x返回矢量的x方向的数值，应力等于力除以面积
    ；线性力不是从0开始计算的
    ；胶结力是从0开始计算的
    ；pb_sigma指的是Normal stress at bond periphery [stress]
end

切向验证

计算切向力和切向应力，并绘制强度包络线。

def jiance   ；进行监测
    whilestepping  ；进行回调
    faxiang_force=contact.force.normal(ct)  ；法向力
    faxiang_yingli=faxiang_force/20.0     ；法向应力
    
    qiexiang_force=contact.force.shear(ct)
    qiexiang_yingli=qiexiang_force/20.0
    
    qiexiang_fengzhi=1e6+math.tan(math.pi*0.25)*faxiang_yingli
    ；计算切向力的峰值，来给出包络线c+σ*tanφ

end

加胶结方式

首先是无胶结生成颗粒，然后再给颗粒添加胶结。

无胶结生成颗粒

new 
domain extent -100 100 

ball create position 0 -9 radius 10
ball create position 0 10 radius 10
；生成颗粒有一定的重叠量

wall generate box -10 10 -19 20    ；给一个墙体约束

ball attribute density 2.3e3 damp 0.7
cmat default model linear property kn 1e7
；生成线性模型

cycle 1
solve
save sample
；松散体内部达到平衡，此时模型内部存在力

加胶结

边界和颗粒之间无胶结，但是颗粒和颗粒之间有胶结。随着颗粒运动，线性弹簧由压迫状态变为松弛；而胶结弹簧一开始就是松弛的，由于可以运动，弹簧被拉伸，产生法向力，当大于一定程度时，将会产生破坏。

负号表示拉伸，正号表示压缩。

把时步调快能够增快计算速度。

cmat对新生成的接触赋值，如果已经有旧的cmat，则需要使用cmat apply。

restore sample
[pb_kn=1e9/20.0]
cmat default type ball-facet model linear property ...
                    kn 1e8 
                    ；设置墙体刚度
cmat default type ball-ball model linearpbond property ...
                    kn 1e9 pb_kn @pb_kn pb_ks @pb_kn ... 
                    pb_ten 1e6 pb_fa 45 pb_coh 1e6 
                    ；kn和ks控制变形
                    ；ten，fa，coh控制强度
cmat apply 
；更新模型后，将会把原先的力进行重置为0，而实际计算时会根据kn*A来重新计算
clean
set timestep fix 1e-5   ；减小时步来查看问题
[ct=contact.find("ball-ball",1)] ；去看contact的标号，查找后为1
def jiance
    whilestepping    
    lin_stress=comp.x(contact.prop(ct,"lin_force"))/20.0 
    ；线性应力，应该有原始积累，全量计算
    pb_stress=contact.prop(ct,"pb_sigma")                
    ；胶结应力，增量计算
end
history id 1 @lin_stress
history id 2 @pb_stress
cycle 1
solve

contact method bond gap 1

cycle 1
solve
save problem

避免胶结破坏的方法

① 将接触力清零

calm
contact property lin_force 0.0 0.0 
ball attribute contactforce multiply 0.0 contactmoment multiply 0.0
;弊端是线性力和胶结力增加的模式不太符合物理实际

② 在改完刚度（cmat apply）之后，再加一个胶结——刚度变化带来的不平衡力计算平衡之后，再添加胶结。

history id 1 @lin_stress
history id 2 @pb_stress
cycle 1
solve

contact method bond gap 1
；之前改变刚度带来的不平衡力已经抵消
；此时再添加胶结，形成bond

cycle 1
solve
save problem