CPT贯入试验

2年前浏览4193

这里粗糙的做了一下CPT贯入的模拟，感觉还有些地方可以斟酌一下，但是整个框架应该是没问题的。

CPT试验是我们经常做的原位测试试验，不用取样进入实验室就可以得到砂土或者黏土的力学特性，通常是分析锥尖阻力来测试土体的力学性质。

模拟的顺序为：成样 - 预压 - 自重 - 染色 - 加锥头 - 加载

一、成样

这部分和之前一样，生成一个矩形试样就可以。


new def par    width=2.5    height=width*0.5        rdmax=0.009    rdmin=0.006        poro=0.08end@pardomain extent [-height*5] [height*5]set random 10001wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5ball distribute porosity @poro radius @rdmin @rdmax box [-width*0.5] [width*0.5] ...                                                        [-height*0.5] [height*0.5]cmat default model linear method deformability emod 100e6 kratio 1.5ball attribute density 2.7e3 damp 0.7cycle 2000 calm 50solve save sample

二、预压

这部分不一定是必要的，但我觉得还是预压一下好，因为一开始试样中的内应力太大了，释放一下不至于自重的时候飞走，而且预压一下试样也分布的更加均匀点。注意这里的预压值为10KPa。



restore sample[txx=-10e3][tyy=-10e3][sevro_factor=0.5][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        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    endifenddef wp_ini    wpDown=wall.find(1)    wpRight=wall.find(2)    wpUp=wall.find(3)    wpLeft=wall.find(4)end@wp_inidef computer_chiCun    wlx=wall.pos.x(wpRight)-wall.pos.x(wpLeft)    wly=wall.pos.y(wpUp)-wall.pos.y(wpDown)enddef 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/wlxenddef get_g(fac)    gx=0    gy=0    zongKNX=100e6*2*10    zongKNY=100e6*2*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)    endset fish callback -1.0 @sevro_wallshistory id 1 @wxsshistory id 2 @wysshistory id 3 @gxhistory id 4 @gycycle 1solve
save yuya

三、加自重

这里直接加上重力就好，如果想做更大尺度的，参考离心机原理即可。


restore yuyaset fish callback -1.0 remove @sevro_wallswall attribute velocity multiply 0 range id 1wall attribute velocity multiply 0 range id 2wall attribute velocity multiply 0 range id 3wall attribute velocity multiply 0 range id 4wall delete walls range id 3set gravity 9.8cycle 1solve aratio 1e-5save zizhong

到这一步我们原始的地基就生成了。

四、染色

这个也不是必要的，但是为了直观点，染个色会好很多。


restore zizhongdef sousuoFanwei    local xyMinMax=array.create(2,2)        local x_min=1e100    local y_min=1e100    local x_max=-1e100    local y_max=-1e100    loop foreach bp ball.list        if x_min>ball.pos.x(bp)-ball.radius(bp) then            x_min=ball.pos.x(bp)-ball.radius(bp)         endif        if x_max<ball.pos.x(bp) ball.radius(bp) then            x_max=ball.pos.x(bp) ball.radius(bp)         endif                if y_min>ball.pos.y(bp)-ball.radius(bp) then            y_min=ball.pos.y(bp)-ball.radius(bp)         endif        if y_max<ball.pos.y(bp) ball.radius(bp) then            y_max=ball.pos.y(bp) ball.radius(bp)         endif    endloop    xyMinMax(1,1)=x_min    xyMinMax(1,2)=x_max    xyMinMax(2,1)=y_min    xyMinMax(2,2)=y_max    sousuoFanwei=xyMinMaxend
def ranse    local xy_min_max=sousuoFanwei    local NWidth=20.0 ;横向方块数目    local NHeight=15.0;纵向数目       local stepwidth=(xy_min_max(1,2)-xy_min_max(1,1))/NWidth    local stepheight=(xy_min_max(2,2)-xy_min_max(2,1))/NHeight        loop local n (1,NWidth)        local minwidth=xy_min_max(1,1) stepwidth*(n-1)                loop local m (1,NHeight)            local minheight=xy_min_max(2,1) stepheight*(m-1)                        if (n m)/2-(n m)/2.0=0 then                command                    ball group gg1 slot 5 range x [minwidth] [minwidth stepwidth] ...                                        y [minheight] [minheight stepheight]                endcommand            else              command                    ball group gg2 slot 5 range x [minwidth] [minwidth stepwidth] ...                                        y [minheight] [minheight stepheight]                endcommand            endif        endloop    endloopend@ransesave ranse

五、生成锥头

这里主要是修改锥杆的宽度和高度，这里默认是45度了，这个需要读者自己琢磨修改一下了。这里考虑半结构。

restore ranse[zuoweizhi=wall.pos.x(wpLeft)][CPTkuandu=width*0.05][CPTHeight=height*0.5*1.02]wall create vertices [zuoweizhi] [CPTHeight] [zuoweizhi CPTkuandu] [CPTHeight CPTkuandu] id 10wall create vertices   [zuoweizhi CPTkuandu] [CPTHeight CPTkuandu]  [zuoweizhi CPTkuandu] [CPTHeight CPTkuandu*5] id 11[wpl=wall.find(10)][wpu=wall.find(11)]cmat add 1 model linear method deformability emod 100e6 kratio 1.5 range contact type ball-facetcmat add 2 model linear method deformability emod 100e6 kratio 1.5 property fric 0.5 range contact type ball-ballcmat applycycle 1solvesave jiaCPT

六、加载

这里直接给速度就可以了。一共运行30s，每1s保存一个sav用于后处理。


restore jiaCPTset mech age 0ball attribute displacement multiply 0wall attribute yvelocity [-height*1e-2] range id 10 11[time_record=mech.age][baocunpinlv=1]def savefile        if mech.age-time_record > baocunpinlv then        filename=string.build("jieguo%1",count)        command            save @filename            plot bitmap filename @filename        endcommand        time_record=mech.age        count  =1    endif    endset fish callback -1.0 @savefiledef jiance    faxiang=math.sqrt(wall.force.contact.x(wpl)*wall.force.contact.x(wpl) wall.force.contact.y(wpl)*wall.force.contact.y(wpl))    cemozuli=wall.force.contact.y(wpu)endset fish callback -1.1 @jiance
history deletehistory id 1 @faxianghistory id 2 @cemozulisolve time 30save result