离散元模拟C30混凝土构件
目前离散元模拟构件（混凝土）有两个方法，第一个是和标定岩石参数一样去标定混凝土参数；第二个是用规则排列的颗粒去模拟混凝土。
第二个方法可以直接将宏观参数和微观参数联系在一起，并且提取内力很方便，所以也被广泛采用。
本文主要是对用规则排列的颗粒进行讲解，一方面讲其宏观参数与微观参数的联系，另一方面讲解如何提取内力。
首先我们采用规则排列的颗粒，加上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_ballbp1=ball.near(vector(-width*0.5,-height*0.5))bp2=ball.near(vector(width*0.5,-height*0.5))ball.fix(bp1,2)=trueball.fix(bp2,2)=trueend@fix_ballwall attribute yvelocity [-height*1e-2] range id 1
给出加载下的位移场：
三、监测变量
这里监测的量有：
1）dip------实测的加载板位移值
2) force------实测的加载板上的力
3）disptheroy------根据加载板的力以及材料力学知识计算的梁中心的扰度。
4） fengzhiP------根据抗拉强度2e6以及材料力学知识计算的最大的力。
我们需要对比的是dip和disptheroy、fengzhiP和实际梁破坏时候的力
def jiancedip=wall.disp.y(wp)force=wall.force.contact.y(wp)I=1*height^3/12.0EI=30e9*Idisptheroy=-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_01def get_jianli_wanjutb_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=0wanju=0loop 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)endlooptable(tb_jianli,-width*0.5 rd*2*(n-1))=posjianlitable(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)=jj1jj2=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)=jj2endloopend@get_jianli_wanju
算0.1s的加载力变化为，可以看出来此时的内力为3.9KN。
根据材料力学或者结构力学知识，剪力应该为其一半，即1.95KN。我们运行上述剪力弯矩代码可以看到，其剪力值中心最大，为1.96KN，两边最小，为1.87KN，不仅吻合了材料力学理论，还体现出了有限元计算结果的规律。
弯矩为剪力的积分，理论最大值应该是1*p*l/4=1.58e3，实测最大值为1.5e3，有一点点出入，这个是由于平截面假定的局限性引起的，这里读者可以自己思考一下。
综上：我们论述了将宏观参数赋值给规则排列颗粒的可行性，也给出了内力的计算方法，所以在离散元中使用颗粒去模拟构件也是可以操作的事情。
