混凝土构件三点弯曲二维模拟
首先祝大家新年快乐，2020年也算过去了，对于我而言每年都有新的挑战，也有新的惊喜。今年最大的挑战就是防疫和毕业了，最大的惊喜就是自己闲暇时候分享的所学知识会有很多人的赞同和认可。
这个从九月份建立已有40多篇原创文章，分享了离散元建模的案例与技巧，几乎开源了所有命令流和代码。也收到很多反馈，也很开心能够帮助很多人学习离散元。迄今为止，几乎包括了岩土工程中的大部分问题，基坑、边坡、隧道、土工试验、裂纹扩展、振动等，也对离散元建模中的一些细致的问题进行了介绍，比如成样方式、阶段保存文件等等。希望大家能够更好的面对2021，也会迎来新的挑战，转角也会有新的惊喜。
今天2021第一天，作为一个岩土学生，也稍微做一点新的东西。因为本科是地质的，所以对混凝土这一块不是特别了解，只是之前在做桩的时候，学了一点混凝土设计原理，了解了一点点。这篇文章简单讲解一下混凝土的三点弯曲试验，主要是解决建模可能出现的问题，对于后面的分析没有涉及，这个可能需要很多比较专业的知识支持。
首先分析了一下这个案例需要解决的技术难点：
1、使用clump模拟骨料
2、解决clump和ball在胶结破坏时候的裂纹问题
一步步来，我们整个程序包括七个步骤，分别是：
生成球颗粒，替换clump，预压、加胶结、卸载、生成加载棒，加载
一、生成球颗粒
之前也介绍过直接生成土石混合体的方法，这里使用另一个方法，先生成球颗粒，之后使用clump替换。命令流为：
newset random 10001def parwidth=2.0height=width*0.2poro=0.03radius1min=0.01radius1max=0.016radius2min=0.002radius2max=0.004end@pardomain extent [-width] [width] [-height] [height]wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5cmat default model linear method deformability emod 100e6 kratio 1.5ball distribute ...porosity @poro ...numbin 2 ...bin 1 ...radius @radius1min @radius1max ...volumefraction 0.3 ...group ball1 ...bin 2 ...radius @radius2min @radius2max ...volumefraction 0.7 ...group ball2 ...range x [-width*0.5 0.001] [width*0.5-0.001] ...y [-height*0.5 0.001] [height*0.5-0.001]ball attribute density 7800 damp 0.7cycle 1000 calm 10solveball delete range x [-width] [-width*0.5]ball delete range x [width*0.5] [width]ball delete range y [-height] [-height*0.5]ball delete range y [height*0.5] [height]save ball_sample
这里使用两个级配表示骨料和砂浆，粗的粒径代表骨料颗粒，生成完后的试样为：
可以看到绿色代表的砂浆颗粒与蓝色代表的骨料颗粒随机分布在试样中。
需要注意的一点是这里没有太关注尺寸和粒径范围，大家可以根据自己需要去取。
二、替换clump
这步算是解决第一个难点了，首先需要导入骨料的形状文件，这里需要读者自己画了，之后生成clump的模板(template)，之后写一个函数用于替换粗颗粒为clump，主要是找出其位置和大小。注意到这里将wall的范围扩大然后在压缩到原来的位置，这个是为了解决clump的自锁现象，当clump有三四个pebble在wall的另一边的时候，clump就会被wall锁住。下面为其命令流：
restore ball_samplegeometry import guliao.dxfclump template create name guliao ...geometry guliao ...bubblepack ratio 0.2 distance 160 ...surfcalculatedef create_clumploop foreach bp ball.groupmap("ball1")x_pos = ball.pos(bp,1)y_pos = ball.pos(bp,2)bvol = math.pi*ball.radius(bp)^2ball.delete(bp)angle= math.random.uniform*180commandclump replicate ...name guliao ...x @x_pos y @y_pos ...volume @bvol angle @angleendcommandendloopend@create_clumpclump attribute density 7800 damp 0.7wall attribute yposition [-height*0.5*1.05] range id 1wall attribute yposition [height*0.5*1.05] range id 3wall attribute xposition [-width*0.5*1.05] range id 4wall attribute xposition [width*0.5*1.05] range id 2wall attribute yvel [height*0.5*0.05*100] range id 1wall attribute yvel [-height*0.5*0.05*100] range id 3wall attribute xvel [width*0.5*0.05*100] range id 4wall attribute xvel [-width*0.5*0.05*100] range id 2solve time 0.01wall attribute vel 0 0cycle 2000 calm 10solvesave clump_sample
替换后的试样如图：
红色的不规则形状便是骨料颗粒。
三、预压
胶结材料都是在一定的压力下形成的，混凝土也不例外。与岩石的高应力下的形成条件不同，混凝土都是在地面上用容器装好凝固的，所以这里的预压值不用设置的很大。
restore clump_sampleball property fric 0.5[txx=-10e3][tyy=-10e3][sevro_factor=0.2][do_xSevro=true][do_ySevro=true][sevro_freq=100][timestepNow=global.step-1]def sevro_wallscompute_stressif timestepNow<global.step thenget_g(sevro_factor)timestepNow =sevro_freqendifif do_xSevro=true thenXvel=gx*(wxss-txx)wall.vel.x(wpRight)=-Xvel; suduwall.vel.x(wpLeft)=Xvelendifif do_ySevro=true thenYvel=gy*(wyss-tyy)wall.vel.y(wpUp)=-Yvelwall.vel.y(wpDown)=Yvelendifenddef wp_iniwpDown=wall.find(1)wpRight=wall.find(2)wpUp=wall.find(3)wpLeft=wall.find(4)end@wp_inidef computer_chiCunwlx=wall.pos.x(wpRight)-wall.pos.x(wpLeft)wly=wall.pos.y(wpUp)-wall.pos.y(wpDown)enddef compute_stresscomputer_chiCunwxss=-(wall.force.contact.x(wpRight)-wall.force.contact.x(wpLeft))*0.5/wlywyss=-(wall.force.contact.y(wpUp)-wall.force.contact.y(wpDown))*0.5/wlxend@compute_stressdef get_g(fac)computer_chiCungx=0gy=0zongKNX=100e6*10zongKNY=100e6*10loop foreach ct wall.contactmap(wpLeft)zongKNX =contact.prop(ct,"kn")endlooploop foreach ct wall.contactmap(wpRight)zongKNX =contact.prop(ct,"kn")endlooploop foreach ct wall.contactmap(wpUp)zongKNY =contact.prop(ct,"kn")endlooploop foreach ct wall.contactmap(wpDown)zongKNY =contact.prop(ct,"kn")endloopgx=fac*wly/(zongKNX*global.timestep)gy=fac*wlx/(zongKNY*global.timestep)end@compute_stressset fish callback -1.0 @sevro_wallshistory id 1 @wxsshistory id 2 @wysscycle 1solvesave yuya
预压后的颗粒以及没有很大的重叠量了。
四、加胶结
模拟的是混凝土的凝固过程，这里和岩石差不多，但是不要忘记给pebble参数：
restore yuya[ball_emod=72e8][pb_emod=ball_emod*0.5][pb_coh=20e6][ten_coh=2.7]cmat default type ball-facet model linear method deformability emod [ball_emod*2] kratio 1.5cmat default type ball-ball model linearpbond method deformability ...emod @ball_emod kratio 1.5 pb_deformability emod @pb_emod kratio 1.5 ...property pb_coh [pb_coh] pb_ten [pb_coh*ten_coh] pb_fa 50 fric 0.5cmat default type ball-pebble model linearpbond method deformability ...emod @ball_emod kratio 1.5 pb_deformability emod @pb_emod kratio 1.5 ...property pb_coh [pb_coh] pb_ten [pb_coh*ten_coh] pb_fa 50 fric 0.5cmat default type pebble-pebble model linearpbond method deformability ...emod @ball_emod kratio 1.5 pb_deformability emod @pb_emod kratio 1.5 ...property pb_coh [pb_coh] pb_ten [pb_coh*ten_coh] pb_fa 50 fric 0.5cmat applycleancycle 1solvecontact method bond gap [radius2min]cycle 1solvesave jiajiaojie
每次加完胶结一定要检查pb_state，确保胶结都加上去了。
五、卸载
这里模拟的是混凝土从容器中取出的过程，给一个比较小的围压即可。
restore jiajiaojie[txx=-1e2][tyy=-1e2]cycle 1solvesave xiezai
六、加加载棒
这里根据试样的尺寸计算加载棒的位置即可
restore xiezaiwall deleteball attribute displacement multiply 0set fish callback -1.0 remove @sevro_walls[wall_radius=wlx*0.04]wall generate id 1 circle position [-wlx*0.5 wall_radius*3] [-wly*0.5-wall_radius] radius [wall_radius]wall generate id 2 circle position 0 [wly*0.5 wall_radius] radius [wall_radius]wall generate id 3 circle position [wlx*0.5-wall_radius*3] [-wly*0.5-wall_radius] radius [wall_radius]cycle 1solvesave addjiazai
加完如图：
七、加载
下面这一步比较简单，直接给中间的加载棒一个竖向的速度即可，注意最好给加载棒一个摩擦，这样可以使得给试样的力是向下的，不然会由于颗粒的离散性导致横向出现移动。
restore addjiazaiwall property fric 0.5[strainrate=1e-2]set mech age 0ball attribute displacement multiply 0wall attribute yvel [-strainrate*wly] range id 2[wp_jiazai=wall.find(2)]def jiancewhilesteppingwall_force=wall.force.contact.y(wp_jiazai)weiyi=wall.disp.y(wp_jiazai)end@jiance[baocunpinlv=0.01][time_record=mech.age-100][count=0]def savefileif mech.age-time_record >= baocunpinlv thenfilename=string.build("jieguo%1",count)commandsave @filenameendcommandtime_record=mech.agecount =1endifendset fish callback -1.0 @savefilehistory deletehistory id 1 @wall_forcehistory id 2 @weiyicall fracture.p2fis@track_initsolve time 2save result
这里最重要的是裂纹文件的改进，原理很简单，对接触两边的颗粒进行判定即可，注意要修改裂纹长度和位置两个地方。
define add_crack(entries)local contact = entries(1)local mode = entries(2)local frac_pos = contact.pos(contact)local norm = contact.normal(contact)local dfn_label = 'crack'local frac_sizelocal bp1 = contact.end1(contact)local bp2 = contact.end2(contact)local ret=0if(type.pointer.name(bp1)=="Ball")if(type.pointer.name(bp2)=="Ball")ret = math.min(ball.radius(bp1),ball.radius(bp2));contact.method(contact,'pb_radius')else if (type.pointer.name(bp2)=="Pebble")ret = math.min(ball.radius(bp1), clump.pebble.radius(bp2))endifelse if(type.pointer.name(bp1)=="Pebble")if(type.pointer.name(bp2)=="Ball")ret = math.min(clump.pebble.radius(bp1),ball.radius(bp2));contact.method(contact,'pb_radius')else if (type.pointer.name(bp2)=="Pebble")ret = math.min(clump.pebble.radius(bp1), clump.pebble.radius(bp2))endifendiffrac_size = retlocal inDir = vector(-comp.y(norm),comp.x(norm))local vert1 = frac_pos inDir * frac_sizelocal vert2 = frac_pos - inDir * frac_sizelocal arg = array.create(4)arg(1) = 'vertices'arg(2) = 2arg(3) = vert1arg(4) = vert2crack_num = crack_num 1jiaojiepohuaineng = entries(4)pb_cohesion=contact.prop(contact,"pb_coh")if mode = 1 then; failed in tensioncrack_tension_num =1dfn_label = dfn_label '_tension'else if mode = 2 then; failed in shearcrack_shear_num =1failforce= entries(3)if failforce<pb_cohesion thendfn_label = dfn_label '_shear_t'crack_shearT_num =1elsedfn_label = dfn_label '_shear_c'crack_shearC_num =1endifendifglobal dfn = dfn.find(dfn_label)if dfn = null thendfn = dfn.add(0,dfn_label)endiflocal fnew = dfn.addfracture(dfn,arg)if mode =1 thencrack_tension_length = dfn.fracture.len(fnew)endifif mode =2 thencrack_shear_length = dfn.fracture.len(fnew)endifcrack_length =dfn.fracture.len(fnew)dfn.fracture.prop(fnew,'age') = mech.agedfn.fracture.extra(fnew,1) = bp1dfn.fracture.extra(fnew,2) = bp2crack_accum = 1if crack_accum > 25if frag_time < mech.agefrag_time = mech.agecrack_accum = 0commandfragment computeendcommandgetfragInfo; go through and update the fracture positionsloop for (local i = 0, i < 3, i = i 1)local name = 'crack_tension'if i = 1name = 'crack_shear_t'endifif i = 2name = 'crack_shear_c'endifdfn = dfn.find(name)if dfn # nullloop foreach local frac dfn.fracturelist(dfn)local ball1 = dfn.fracture.extra(frac,1)local ball2 = dfn.fracture.extra(frac,2)if ball1 # nullif ball2 # nulllocal len = dfn.fracture.len(frac)/2.0local pos=0if(type.pointer.name(ball1)=="Ball")if(type.pointer.name(ball2)=="Ball")pos = (ball.pos(ball1) ball.pos(ball2))/2.0;contact.method(contact,'pb_radius')else if (type.pointer.name(ball2)=="Pebble")pos = (ball.pos(ball1) clump.pebble.pos(ball2))/2.0endifelse if(type.pointer.name(ball1)=="Pebble")if(type.pointer.name(ball2)=="Ball")pos = ( clump.pebble.pos(ball1) ball.pos(ball2))/2.0;contact.method(contact,'pb_radius')else if (type.pointer.name(ball2)=="Pebble")pos = (clump.pebble.pos(ball1) clump.pebble.pos(ball2))/2.0endifendifif comp.x(pos)-len > xminif comp.x(pos) len < xmaxif comp.y(pos)-len > yminif comp.y(pos) len < ymaxdfn.fracture.pos(frac) = posendifendifendifendifendifendifendloopendifendloopendifendifenddefine track_initcommanddfn deleteball result clearfragment clearfragment register ball-ballfragment register ball-pebblefragment register pebble-pebblefragment computeendcommand; activate fishcallscommandset fish callback bond_break remove @add_crackset fish callback bond_break @add_crackendcommand; reset global variablesglobal crack_accum = 0global crack_num = 0global crack_tension_num=0global crack_shear_num=0global crack_shearT_num=0global crack_shearC_num=0global crack_tension_length=0global crack_shear_length=0global crack_length=0global fragNum=0global max_frage_vol=1global jiaojiepohuaineng=0global track_time0 = mech.ageglobal frag_time = mech.ageglobal xmin = domain.min.x()global ymin = domain.min.y()global xmax = domain.max.x()global ymax = domain.max.y()commandhistory id 50 @crack_numhistory id 51 @crack_tension_numhistory id 52 @crack_shear_numhistory id 53 @crack_lengthhistory id 54 @crack_tension_lengthhistory id 55 @crack_shear_lengthhistory id 56 @fragNumhistory id 57 @max_frage_volhistory id 58 @jiaojiepohuainenghistory id 59 @crack_shearT_numhistory id 60 @crack_shearC_numendcommandenddef getfragInfolocal frag1 = fragment.find(1)local catalognow=fragment.catalog.num(mech.age)fragNum=fragment.num(catalognow)max_frage_vol = fragment.vol(frag1,catalognow) / fragment.vol(frag1,1)end
结果分析：
1、位移场
下图为加载0.7s颗粒的位移场，可以发现中间位移为负值，两边位移为正值，也是符合我们的常见的简支梁的位移场的。
2、裂纹发展
这个是使用离散元做混凝土最关心的问题，这种工况毫无疑问裂纹是从底部开始发展的。混凝土作为一种脆性材料，其破坏肯定是突然的，我们先看一下加载到0.6s试样的状况，这时候还没有破坏，从力的图也可以看出来。
加载到0.68s时，试样依然没有破坏
加载到0.69s时，试样在底部出现了裂纹，力也出现了衰减：
放大点看，可以看到混凝土中裂纹是绕着骨料发育的，并且和岩石这样的胶结材料一样以拉剪裂纹为主。
加载到0.7s，裂纹有了进一步的发育
加载到0.8s，裂纹试样几乎被完全折断，并且拉裂纹也开始发育：
3、破坏能
混凝土破坏研究这个应该是必不可少的，这里使用胶结破坏能代替，这个单个工况只能发现胶结材料的突变破坏，横向对比还是很有意义的。
4、拉剪接触分布
打开力链图，便可以看到拉剪接触的分布情况，下面是加载0.6s的状态，和常规的材料力学所解释的一样，下方为受拉侧，只是失去了平截面假设后，这种分布变得不均匀。但是规律性很好。
5、应力方向
通过tensor绘制应力十字架图，可以发现主应力的方向，这个和常识也是没有偏差的。
6、接触破坏情况
这是算了1s的状态图，通过接触看试样破坏是比较直观的，毕竟颗粒是离散的，但是接触是连续的。很明显的发现试样在三点弯曲加载下的破坏情况。
好了，新的一年开始了，也很感谢各位的关注，出去哈皮去了hhhhhh
登录后免费查看全文
著作权归作者所有，欢迎分享，未经许可，不得转载
首次发布时间：2022-07-20
最近编辑：1年前
硕士
|
无
擅长颗粒流PFC
作者推荐
相关推荐
热门文章
