基于PFC3D與Flac3D耦合的SHPB壓桿模擬

lobby

2022年3月7日 13:59

瀏覽：5345 評論：3 收藏：28

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"

生成的式樣如圖

基于PFC3D與Flac3D耦合的SHPB壓桿模擬的圖1

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"

生產的壓桿如圖：

基于PFC3D與Flac3D耦合的SHPB壓桿模擬的圖2

細節放大：

基于PFC3D與Flac3D耦合的SHPB壓桿模擬的圖3

可以發現此時顆粒是有一定速度的，但是由于我們研究的是高應變率所以這里可以不用管。

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"