首页/文章/ 详情

离散元模拟C30混凝土构件

lobby

2年前浏览4418

目前离散元模拟构件（混凝土）有两个方法，第一个是和标定岩石参数一样去标定混凝土参数；第二个是用规则排列的颗粒去模拟混凝土。

第二个方法可以直接将宏观参数和微观参数联系在一起，并且提取内力很方便，所以也被广泛采用。

本文主要是对用规则排列的颗粒进行讲解，一方面讲其宏观参数与微观参数的联系，另一方面讲解如何提取内力。

首先我们采用规则排列的颗粒，加上pb模型去生成混凝土桩，这里的混凝土参数取C30混凝土的参数。弹性模量30e9，抗拉强度2MPa。因为桩都是受弯发生拉破坏，所以这里不考虑压坏。

下面一步步讲解建模思路。

一、生成混凝土

第一步我们需要混凝土试样，试样首先需要规则排列的颗粒，这里定义了颗粒的半径，竖向颗粒个数，并计算出竖向的混凝土的高，之后又定义了长是高的15倍。

我们给出了混凝土的密度，并且设置了默认接触参数，这个参数会在后面加载板与混凝土间起作用。

new[rd=0.006][shuxianggeshu=9][height=rd*2*shuxianggeshu][width=height*15]domain extent [-width*2] [width*2] [-height*2] [height*2]ball generate box [-width*0.5 rd] [width*0.5-rd] [-height*0.5 rd] [height*0.5-rd] radius [rd] cubicball attribute density 3e3 damp 0.7cmat default model linear method deform emod 10e9 kratio 1.5

这部分代码运算完会生成1215个颗粒

第二步我们需要给定胶结，规则排列的颗粒是可以直接将宏观模量给模型的，这里不给定颗粒的刚度(emod)，而是将所有的模量直接给胶结刚度(pb_emod)，这样就避免了受压和受拉模量不一样的问题。泊松比这个对规则排列好像是没有用的，这里也是随意赋值，这里的pb_kra也能看出其与泊松比的理论关系。


[ehongguan=30e9][pb_em=ehongguan][vhongguan=0.4][pb_kra=1/vhongguan]cmat add 1 model linearpbond method deformability emod 0 kratio 1.5 pb_deformability emod [pb_em] kratio [pb_kra] ...        property pb_coh 20e6 pb_ten 2e6 pb_fa 50 fric 0.5 range contact type ball-ballcleancontact method bond gap [rd*0.2]

运行完后胶结加上，pb_state=3，说明胶结都加上了。到这里完整的混凝土试样完成了。

二、简支梁加载条件

为了去验证模型参数，我们对其做一个简支梁的模型验证。简支梁为两端铰接，中间加载。

这里生成一个矩形的box作为加载板，对左下角和右下角的颗粒施加竖向约束模拟铰接情况，如果约束横向的话，整个一列的颗粒都无法转动，这样就变成固定了，这里可以自己尝试理解一下。

最后给加载板一个向下的速度


[box_height=height*0.2]wall generate box [-rd] [rd] [height*0.5] [height box_height] cleandef fix_ball    bp1=ball.near(vector(-width*0.5,-height*0.5))    bp2=ball.near(vector(width*0.5,-height*0.5))         ball.fix(bp1,2)=true          ball.fix(bp2,2)=true end @fix_ballwall attribute yvelocity [-height*1e-2] range id 1

给出加载下的位移场：

三、监测变量

这里监测的量有：

1）dip------实测的加载板位移值

2) force------实测的加载板上的力

3）disptheroy------根据加载板的力以及材料力学知识计算的梁中心的扰度。

4） fengzhiP------根据抗拉强度2e6以及材料力学知识计算的最大的力。

我们需要对比的是dip和disptheroy、fengzhiP和实际梁破坏时候的力


def jiance     dip=wall.disp.y(wp)    force=wall.force.contact.y(wp)        I=1*height^3/12.0    EI=30e9*I    disptheroy=-force*width^3/(EI*48.0)        fengzhiP=(2e6*I/(height*0.5))*4/widthendset fish callback -1.0 @jiancehistory id 1 @diphistory id 2 @forcehistory id 3 @disptheroysolve time 3e-1save result

四、扰度对比

绿色线为实测值，蓝色线为理论（因为力在变化，所以理论扰度也在变），可以发现蓝色线在破坏前一直围绕绿色线变化，这个是因为加载速度过快，体现出了离散元计算的动力特性。总的来说理论扰度和实测扰度还是比较吻合的。

这个说明直接将宏观弹性模量给pb_emod是可行的

五、强度对比

根据force的变化，可以看出其峰值力为9.7KN，理论值为9.6KN，有一点点出入，但是误差只有1%左右，还是可以接受的。

这个说明直接将宏观抗拉强度给pb_ten是可行的

六、剪力弯矩对比

因为破坏后的内力分布没办法和理论做对比，这里只算0.1s，然后分析内力。

这里给出剪力和弯矩的计算方法，大概的逻辑是剪力就是横向两个颗粒间的Y向力，而弯矩是x向力对中性轴的积分，这里遵循平截面假定，认为中性轴在中间。

这里使用UDVector来绘制剪力和弯矩，注意一次只能绘制其中一个，不然会重叠在一起。


restore result_01
def get_jianli_wanju      tb_jianli=table.create("jianli")    tb_wanju=table.create("wanju")    hengxiang=math.ceiling( width*0.5/rd)    zongxiang=math.ceiling(height*0.5/rd)        array jianliarr(200)    array wanjuarr(200)    loop  n(1,hengxiang-1)        posjianli=0        wanju=0        loop  m(1,zongxiang)                        ct=contact.find("ball-ball",hengxiang*(m-1) n,hengxiang*(m-1) n 1)            posjianli =contact.force.local.y(ct)            wanju =contact.pos.y(ct)*contact.force.local.x(ct)        endloop        table(tb_jianli,-width*0.5 rd*2*(n-1))=posjianli        table(tb_wanju,-width*0.5 rd*2*(n-1))=wanju                       ; jj1=user.vector.create(vector(-width*0.5 rd*2*n,0));        user.vector.value(jj1)=vector(0, posjianli);        user.vector.group(jj1) = "jianli";        jianliarr(n)=jj1                jj2=user.vector.create(vector(-width*0.5 rd*2*n,0))        user.vector.value(jj2)=vector(0, wanju)        user.vector.group(jj2) = "wanju"        wanjuarr(n)=jj2    endloopend@get_jianli_wanju