混凝土构件三点弯曲二维模拟

2年前浏览3327

首先祝大家新年快乐，2020年也算过去了，对于我而言每年都有新的挑战，也有新的惊喜。今年最大的挑战就是防疫和毕业了，最大的惊喜就是自己闲暇时候分享的所学知识会有很多人的赞同和认可。

这个从九月份建立已有40多篇原创文章，分享了离散元建模的案例与技巧，几乎开源了所有命令流和代码。也收到很多反馈，也很开心能够帮助很多人学习离散元。迄今为止，几乎包括了岩土工程中的大部分问题，基坑、边坡、隧道、土工试验、裂纹扩展、振动等，也对离散元建模中的一些细致的问题进行了介绍，比如成样方式、阶段保存文件等等。希望大家能够更好的面对2021，也会迎来新的挑战，转角也会有新的惊喜。

今天2021第一天，作为一个岩土学生，也稍微做一点新的东西。因为本科是地质的，所以对混凝土这一块不是特别了解，只是之前在做桩的时候，学了一点混凝土设计原理，了解了一点点。这篇文章简单讲解一下混凝土的三点弯曲试验，主要是解决建模可能出现的问题，对于后面的分析没有涉及，这个可能需要很多比较专业的知识支持。

首先分析了一下这个案例需要解决的技术难点：

1、使用clump模拟骨料

2、解决clump和ball在胶结破坏时候的裂纹问题

一步步来，我们整个程序包括七个步骤，分别是：

生成球颗粒，替换clump，预压、加胶结、卸载、生成加载棒，加载

一、生成球颗粒

之前也介绍过直接生成土石混合体的方法，这里使用另一个方法，先生成球颗粒，之后使用clump替换。命令流为：

new

set random 10001def par    width=2.0    height=width*0.2    poro=0.03     radius1min=0.01    radius1max=0.016     radius2min=0.002    radius2max=0.004  end@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.5 ball 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 10solve ball 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 ...                        surfcalculate def create_clump    loop foreach bp ball.groupmap("ball1")                       x_pos   = ball.pos(bp,1)            y_pos   = ball.pos(bp,2)            bvol = math.pi*ball.radius(bp)^2               ball.delete(bp)               angle= math.random.uniform*180             command                                clump replicate    ...                            name guliao ...                            x @x_pos y @y_pos   ...                           volume @bvol  angle @angle            endcommand                 endloopend@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 2

wall attribute yvel [height*0.5*0.05*100] range id 1wall attribute yvel [-height*0.5*0.05*100] range id 3

wall 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 0 cycle 2000 calm 10solvesave clump_sample

替换后的试样如图：

红色的不规则形状便是骨料颗粒。

三、预压

胶结材料都是在一定的压力下形成的，混凝土也不例外。与岩石的高应力下的形成条件不同，混凝土都是在地面上用容器装好凝固的，所以这里的预压值不用设置的很大。

restore clump_sample

ball 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_walls    compute_stress    if timestepNow<global.step then        get_g(sevro_factor)        timestepNow =sevro_freq    endif    if do_xSevro=true then        Xvel=gx*(wxss-txx)        wall.vel.x(wpRight)=-Xvel;   sudu        wall.vel.x(wpLeft)=Xvel    endif     if do_ySevro=true then        Yvel=gy*(wyss-tyy)        wall.vel.y(wpUp)=-Yvel        wall.vel.y(wpDown)=Yvel    endifend

def wp_ini    wpDown=wall.find(1)    wpRight=wall.find(2)    wpUp=wall.find(3)    wpLeft=wall.find(4)end@wp_ini

def computer_chiCun    wlx=wall.pos.x(wpRight)-wall.pos.x(wpLeft)    wly=wall.pos.y(wpUp)-wall.pos.y(wpDown)end

def compute_stress    computer_chiCun    wxss=-(wall.force.contact.x(wpRight)-wall.force.contact.x(wpLeft))*0.5/wly    wyss=-(wall.force.contact.y(wpUp)-wall.force.contact.y(wpDown))*0.5/wlxend@compute_stress

def get_g(fac)    computer_chiCun    gx=0    gy=0    zongKNX=100e6*10    zongKNY=100e6*10    loop foreach ct wall.contactmap(wpLeft)        zongKNX =contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpRight)        zongKNX =contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpUp)        zongKNY =contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpDown)        zongKNY =contact.prop(ct,"kn")    endloop    gx=fac*wly/(zongKNX*global.timestep)    gy=fac*wlx/(zongKNY*global.timestep) end

@compute_stress

set fish callback -1.0 @sevro_walls

history id 1 @wxsshistory id 2 @wysscycle 1solve save yuya

预压后的颗粒以及没有很大的重叠量了。

微信截图_20220720153009.png

四、加胶结

模拟的是混凝土的凝固过程，这里和岩石差不多，但是不要忘记给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.5 cmat 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 applyclean

cycle 1solvecontact method bond gap [radius2min]

cycle 1 solvesave jiajiaojie

每次加完胶结一定要检查pb_state，确保胶结都加上去了。

五、卸载

这里模拟的是混凝土从容器中取出的过程，给一个比较小的围压即可。

restore jiajiaojie

[txx=-1e2][tyy=-1e2]

cycle 1 solvesave xiezai

六、加加载棒

这里根据试样的尺寸计算加载棒的位置即可

restore xiezai

wall 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 1 solvesave addjiazai

加完如图：

七、加载

下面这一步比较简单，直接给中间的加载棒一个竖向的速度即可，注意最好给加载棒一个摩擦，这样可以使得给试样的力是向下的，不然会由于颗粒的离散性导致横向出现移动。



restore addjiazaiwall property fric 0.5[strainrate=1e-2]set mech age 0ball attribute displacement multiply 0

wall attribute yvel [-strainrate*wly] range id 2

[wp_jiazai=wall.find(2)]

def jiance    whilestepping    wall_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 savefile        if mech.age-time_record >= baocunpinlv then        filename=string.build("jieguo%1",count)        command            save @filename        endcommand        time_record=mech.age        count  =1    endif    endset fish callback -1.0 @savefile

history deletehistory id 1 @wall_forcehistory id 2 @weiyicall fracture.p2fis@track_init

solve time 2

save 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_size    local bp1 = contact.end1(contact)    local bp2 = contact.end2(contact)    local ret=0     if(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))        endif    else 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))        endif    endif    frac_size = ret    local inDir = vector(-comp.y(norm),comp.x(norm))    local vert1 = frac_pos   inDir * frac_size    local vert2 = frac_pos - inDir * frac_size    local arg = array.create(4)    arg(1) = 'vertices'    arg(2) = 2    arg(3) = vert1    arg(4) = vert2    crack_num = crack_num   1    jiaojiepohuaineng = entries(4)      pb_cohesion=contact.prop(contact,"pb_coh")   if mode = 1 then        ; failed in tension

         crack_tension_num =1        dfn_label = dfn_label   '_tension'

    else if mode = 2 then        ; failed in shear        crack_shear_num =1        failforce= entries(3)        if failforce<pb_cohesion then             dfn_label = dfn_label   '_shear_t'            crack_shearT_num =1        else             dfn_label = dfn_label   '_shear_c'            crack_shearC_num =1        endif    endif      global dfn = dfn.find(dfn_label)    if dfn = null then        dfn = dfn.add(0,dfn_label)    endif    local fnew = dfn.addfracture(dfn,arg)     if mode =1 then        crack_tension_length = dfn.fracture.len(fnew)    endif     if mode =2 then        crack_shear_length = dfn.fracture.len(fnew)    endif    crack_length =dfn.fracture.len(fnew)    dfn.fracture.prop(fnew,'age')  = mech.age    dfn.fracture.extra(fnew,1) = bp1    dfn.fracture.extra(fnew,2) = bp2    crack_accum  = 1    if crack_accum > 25        if frag_time < mech.age            frag_time = mech.age            crack_accum = 0            command                 fragment compute            endcommand            getfragInfo            ; go through and update the fracture positions            loop for (local i = 0, i < 3, i = i   1)                local name = 'crack_tension'                if i = 1                    name = 'crack_shear_t'                endif                if i = 2                    name = 'crack_shear_c'                endif                 dfn = dfn.find(name)                if dfn # null                    loop foreach local frac dfn.fracturelist(dfn)                        local ball1 = dfn.fracture.extra(frac,1)                        local ball2 = dfn.fracture.extra(frac,2)                        if ball1 # null                            if ball2 # null                                local len = dfn.fracture.len(frac)/2.0                                local pos=0                                 if(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.0                                    endif                                else 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.0                                    endif                                endif                                if comp.x(pos)-len > xmin                                    if comp.x(pos) len < xmax                                        if comp.y(pos)-len > ymin                                            if comp.y(pos) len < ymax                                                dfn.fracture.pos(frac) = pos                                            endif                                        endif                                    endif                                endif                            endif                        endif                    endloop                endif            endloop        endif    endifend

define track_init    command        dfn delete        ball result clear        fragment clear        fragment register ball-ball        fragment register ball-pebble        fragment register pebble-pebble           fragment compute    endcommand  ; activate fishcalls    command        set fish callback bond_break remove @add_crack        set fish callback bond_break @add_crack    endcommand    ; reset global variables    global crack_accum = 0    global crack_num = 0    global crack_tension_num=0    global crack_shear_num=0    global crack_shearT_num=0    global crack_shearC_num=0    global crack_tension_length=0    global crack_shear_length=0    global crack_length=0     global fragNum=0    global max_frage_vol=1    global jiaojiepohuaineng=0     global track_time0 = mech.age    global frag_time = mech.age    global xmin = domain.min.x()    global ymin = domain.min.y()    global xmax = domain.max.x()    global ymax = domain.max.y()     command        history id 50 @crack_num        history id 51 @crack_tension_num        history id 52 @crack_shear_num        history id 53 @crack_length        history id 54 @crack_tension_length        history id 55 @crack_shear_length

        history id 56 @fragNum        history id 57 @max_frage_vol                history id 58 @jiaojiepohuaineng          history id 59 @crack_shearT_num        history id 60 @crack_shearC_num    endcommand

end

def getfragInfo    local 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