目录
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