基于PFC3D与Flac3D耦合的SHPB压杆模拟

2年前浏览6029

0 引言

最近在熟悉PFC6.0的内容，对于6.0来说，其最大的突破就是实现了软件内与FLAC3D的耦合。这其实解决了PFC很大的一个问题，就是颗粒数太多时计算力的限制，我们可以用网格单元作为边界来弱化边界效应，当然也可以和本文一样，使用网格单元模拟金属等连续性材料。

本文试着去模拟一个结构方向同学经常遇到的一个工况，霍普金森压杆（SHPB）实验，对于测量混凝土材料在高应变率下的力学特性有很大的参考价值。因为我这里也没有看过这方面的实验，仅从有限的资料大概理解其边界设置，利用有限差分和离散元耦合来实现SHPB压杆模拟，希望能给各位有一定的参考价值。

1 成样、预压、加胶结

首先是我们的离散元部分，这里进行常规的方式来生成一个圆柱形式样。这里不再赘述了。

成样代码：

model new

fish define Par    sample_width=0.05      radmin=0.8e-3    radmax=radmin*1.8    poro=0.28      sample_rad=sample_width*0.5      poleLength=sample_width*40end@Parmodel domain extent [-sample_width]  [sample_width]

 wall generate cylinder base [-sample_width*0.5*1.4] 0 0 axis 1 0 0  ...                    height [sample_width*1.4] radius [sample_rad] resolution 0.3 cap false false

wall generate plane position [sample_width*0.5] 0 0  dip 90 dip-dir 90wall generate plane position [-sample_width*0.5] 0 0  dip 90 dip-dir 90

ball distribute group "shiyang" radius [radmin] [radmax] porosity @poro ...    range cylinder end-1 [sample_width*0.5-radmax] 0 0  ...    end-2 [-sample_width*0.5 radmax] 0 0  radius [sample_rad-radmax]

cmat default type ball-facet model linear method deform emod 100e6 kratio 1.5 cmat default type ball-ball model linear method deform emod 100e6 kratio 1.5  ball attribute density 2.7e3 damp 0.2model cycle 2000 calm 50

model solve

model save "sample"

预压代码：

model restore "sample"

ball property "fric" 0.5

fish def wp_wall    wp_up=wall.find(2)    wp_down=wall.find(3)    wp_rr=wall.find(1)        loop foreach vt wall.vertexlist(wp_rr)        vert_in_ce=vt             endloopend@wp_wall

[txx=-1e6][trr=-1e6][sevro_fac=0.5]

[do_xservo=true][do_rservo=true]



[sevro_freq=100][timestepNow=global.step-1]fish def servo_walls    computer_wallStress    if timestepNow<global.step then        get_gain(sevro_fac)        timestepNow =sevro_freq    endif    if do_xservo=true then        x_vel=gx*(wsxx-txx)        wall.vel.x(wp_up)=-x_vel        wall.vel.x(wp_down)=x_vel    endif     if do_rservo=true then        r_vel_mag=(-1)*gr*(wsrr-trr)        loop foreach vt wall.vertexlist(wp_rr)            mag=math.sqrt(wall.vertex.pos.y(vt)^2 wall.vertex.pos.z(vt)^2)            fang_normal_y=wall.vertex.pos.y(vt)/mag            fang_normal_z=wall.vertex.pos.z(vt)/mag                        r_vel=vector(0,fang_normal_y,fang_normal_z)*r_vel_mag                        wall.vertex.vel(vt)=r_vel         endloop    endif     end



fish def computer_chicun    y_pos=wall.vertex.pos.y(vert_in_ce)    z_pos=wall.vertex.pos.z(vert_in_ce)    wlr=math.sqrt(y_pos^2 z_pos^2)    wlx=wall.pos.x(wp_up)-wall.pos.x(wp_down)end

fish def computer_wallStress    computer_chicun    ding_yuanmianji=math.pi*wlr^2    wsxx=(wall.force.contact.x(wp_down)-wall.force.contact.x(wp_up))*0.5/ding_yuanmianji    ce_mianji=2*math.pi*wlr*wlx    wsrr=0    loop foreach ft wall.facetlist(wp_rr)        ft_fangxiang=wall.facet.normal(ft)        loop foreach ct wall.facet.contactmap(ft)            force_in_facet=contact.force.global(ct)            wsrr =(math.dot(force_in_facet,ft_fangxiang))/ce_mianji               endloop    endloop end

def get_gain(fac)    gx=0    gr=0    zonggangX=1e7    zonggangR=1e7    loop foreach ct wall.contactmap(wp_up)        zonggangX =contact.prop(ct,"kn")    endloop     loop foreach ct wall.contactmap(wp_down)        zonggangX =contact.prop(ct,"kn")    endloop     loop foreach ct wall.contactmap(wp_rr)        zonggangR =contact.prop(ct,"kn")    endloop     gX=fac*ding_yuanmianji/(zonggangX*global.timestep)    gr=fac*ce_mianji/(zonggangR*global.timestep)    end





fish callback add  @servo_walls -1.0

history id 1 @wsXXhistory id 2 @wsrr

model cycle 1model solvemodel save "yuya"

加胶结代码

model restore "yuya"

[pb_coh=100e6][ten_coh=2.7]contact cmat default type ball-facet model linear method deformability emod 10e8 kratio 1.5



contact cmat default type ball-ball model linearpbond method deformability ...            emod 12e8 kratio 1.5 pb_deformability emod 54e8 kratio 1.5 ...    property "pb_coh" [pb_coh] "pb_ten" [pb_coh*ten_coh] "pb_fa" 50 "fric" 0contact cmat apply



model cycle 1model solve

contact method bond gap [radmin*0.5]

model cycle 1 model solve

model save "jiajiaojie"

生成的式样如图

2、卸载

这里模拟混凝土从模具中取出，或者岩石从地下取出的过程，所以预压的值是有讲究的。

model restore "jiajiaojie"

[txx=-1e3][trr=-1e3]

model cycle 1 model solve

model save "xiezai"

3、生成压杆

这部分需要生产入射杆和透射杆的单元，目前我知道的方式是通过四个1/4的cylinder拼起来，不知道会不会有更便捷的方法，可以后台交流。

之后就是耦合部分的定义，PFC和FLAC提供了两种耦合方式，wall-zone耦合与ball-zone耦合，就使用来看，本文应当用wall-zone耦合来实现zone与ball之间的力学联系。

需要注意的一点是，由于我们的杆件是拼起来的，所以zone生成wall时，会出现重复点的情况，这里可以直接指定skip-errors跳过即可。

model restore "xiezai"

model domain extent [-poleLength*1.5] [poleLength*1.5] ...                    [-sample_width]  [sample_width] wall delete fish callback remove  @servo_walls -1.0[materialMod=200e9]

model largestrain on

[bullet_Length=sample_width*8][midu=5]

zone create cylinder point 0 ([-poleLength-wlx*0.5],0,0) ...                    point 1 ([-poleLength-wlx*0.5],[-sample_rad],0) ...                    point 2 ([-wlx*0.5],0,0) ...                    point 3 ([-poleLength-wlx*0.5],0,[sample_rad]) ...                    size   @midu [math.floor(midu*poleLength/sample_rad)]  @midu group "ImportPole"zone create cylinder point 0 ([-poleLength-wlx*0.5],0,0) ...                    point 1 ([-poleLength-wlx*0.5],0,[-sample_rad]) ...                    point 2 ([-wlx*0.5],0,0) ...                    point 3 ([-poleLength-wlx*0.5],[-sample_rad],0) ...                    size    @midu [math.floor(midu*poleLength/sample_rad)]  @midu group "ImportPole"zone create cylinder point 0 ([-poleLength-wlx*0.5],0,0) ...                    point 1 ([-poleLength-wlx*0.5],[sample_rad],0) ...                    point 2 ([-wlx*0.5],0,0) ...                    point 3 ([-poleLength-wlx*0.5],0,[-sample_rad]) ...                    size   @midu [math.floor(midu*poleLength/sample_rad)]  @midu group "ImportPole"zone create cylinder point 0 ([-poleLength-wlx*0.5],0,0) ...                    point 1 ([-poleLength-wlx*0.5],0,[sample_rad]) ...                    point 2 ([-wlx*0.5],0,0) ...                    point 3 ([-poleLength-wlx*0.5],[sample_rad],0) ...                    size   @midu [math.floor(midu*poleLength/sample_rad)]  @midu group "ImportPole"





; toushezone create cylinder point 0 ([wlx*0.5],0,0)   ...                    point 1 ([wlx*0.5],[-sample_rad],0) ...                    point 2 ([poleLength wlx*0.5],0,0) ...                    point 3 ([wlx*0.5],0,[sample_rad]) ...                    size   @midu [math.floor(midu*poleLength/sample_rad)]  @midu group "ExportPole"zone create cylinder point 0 ([wlx*0.5],0,0) ...                    point 1 ([wlx*0.5],0,[-sample_rad]) ...                    point 2 ([poleLength wlx*0.5],0,0) ...                    point 3 ([wlx*0.5],[-sample_rad],0) ...                    size   @midu [math.floor(midu*poleLength/sample_rad)]  @midu group "ExportPole"zone create cylinder point 0 ([wlx*0.5],0,0) ...                    point 1 ([wlx*0.5],[sample_rad],0) ...                    point 2 ([poleLength wlx*0.5],0,0) ...                    point 3 ([wlx*0.5],0,[-sample_rad]) ...                    size   @midu [math.floor(midu*poleLength/sample_rad)]  @midu group "ExportPole"zone create cylinder point 0 ([wlx*0.5],0,0) ...                    point 1 ([wlx*0.5],0,[sample_rad]) ...                    point 2 ([poleLength wlx*0.5],0,0) ...                    point 3 ([wlx*0.5],[sample_rad],0) ...                    size   @midu [math.floor(midu*poleLength/sample_rad)]  @midu group "ExportPole"

zone cmodel assign elasticzone property young [materialMod] poisson 0.25



wall-zone create skip-errors range position-x [-wlx*0.5-0.001] [-wlx*0.5 0.001]wall-zone create skip-errors range position-x [wlx*0.5-0.001] [wlx*0.5 0.001]



model save "pole"

生产的压杆如图：

细节放大：

可以发现此时颗粒是有一定速度的，但是由于我们研究的是高应变率所以这里可以不用管。

4、加载

这里直接施加一个冲击荷载，在入射杆的端部指定一个200MPa的应力，持续200us。之后接触边界条件，运算1500us。

model restore "pole"

model mechanical time-total 0model configure dynamic



ball attribute displacement multiply 0

zone face  group "rushe" range position-x [-poleLength-wlx*0.5-0.001] ...    [-poleLength-wlx*0.5 0.001]

zone face apply stress-xx [(-1)*(200e6)] range group "rushe"



zone initialize density 7.8e3model cycle 1

measure create position 0 0 0 radius [wlr*0.5]def cedian    zoneImportPole=zone.near(vector(-wlx*0.5-poleLength*0.5,[-wlr],0))    zoneExportPole=zone.near(vector(wlx*0.5 poleLength*0.5,[-wlr],0))    mp1=measure.find(1)       zone.group(zoneImportPole,"cedianGroup")="cedian"    zone.group(zoneExportPole,"cedianGroup")="cedian"    stressShiyang0=measure.stress.xx(mp1)end@cedian

def jiance    time=mech.time.total    stressImport=zone.stress.xx(zoneImportPole)    stressExport=zone.stress.xx(zoneExportPole)    stressShiyang=measure.stress.xx(mp1)-stressShiyang0endfish callback add  @jiance -1.0history deletehistory id 1 @stressImporthistory id 2 @stressExporthistory id 3 @stressShiyang

history id 4 @time[baocunpinlv=500e-6][time_record=-10][count=0]def savefile        if time-time_record >= baocunpinlv then        filename=string.build("jieguo%1",count)        command            model save @filename        endcommand        time_record=time        count  =1    endif    end                        fish callback add @savefile -1.0 program call "fracture.p3fis"@track_initmodel solve time 200e-6zone face apply-remove stress-xx range group "rushe"model solve time 1500e-6model save "jieguo"