砂雨模型制備PFC3D 的數(shù)值模擬

王浩然，王永志，王 海，湯兆光，王體強(qiáng)

(1.中國(guó)地震局工程力學(xué)研究所 地震工程與工程振動(dòng)重點(diǎn)試驗(yàn)室，黑龍江 哈爾濱 150080；2.濟(jì)南大學(xué) 水利與環(huán)境學(xué)院，山東 濟(jì)南 250022)

砂雨法，也稱(chēng)落雨法，因能較真實(shí)地模擬自然風(fēng)力、水力等作用下天然粗粒土的沉積過(guò)程，從而獲得與天然粗粒土地基較為近似的力學(xué)特性；同時(shí)，其與插搗法、振動(dòng)法、擊實(shí)法等相比，可有效避免顆粒破碎和獲得較寬密實(shí)度范圍，被廣泛應(yīng)用于室內(nèi)物理試驗(yàn)?zāi)Ｐ椭苽洌绕涫莻鹘y(tǒng)單元試驗(yàn)和土工離心模型試驗(yàn)[1-6]。通常，砂雨法利用落距、移速、出砂口尺寸、環(huán)境介質(zhì)、篩網(wǎng)數(shù)量等因素控制，實(shí)現(xiàn)預(yù)期粗粒土模型孔隙率和均勻性的控制目標(biāo)。Lagioia 等[4]利用模型制備試驗(yàn)，探討了空氣、水、真空等環(huán)境因素對(duì)砂雨法制備模型孔隙率的影響，結(jié)果表明水及真空環(huán)境分別可獲得最大、最小孔隙率。Cresswell 等[5]通過(guò)流量、移速和落距等因素影響規(guī)律的對(duì)比，提出了砂雨法制備最小孔隙率的基本原則與要求。Miura 等[6-7]為認(rèn)識(shí)砂雨法控制因素影響特征，開(kāi)展了一系列粗粒土制備試驗(yàn)，指出流量為控制孔隙率的首要因素，而孔隙率的控制精度和均勻性隨孔隙率增大而降低。馬險(xiǎn)峰等[8-9]介紹了我國(guó)8 所院校開(kāi)展的砂雨法制備平行試驗(yàn)，給出了出砂頭類(lèi)型及孔徑、流量、移速、落距等因素影響基本規(guī)律。王海等[10]對(duì)比了干砂與飽和砂模型砂雨法制備控制因素的影響差異，指出了兩類(lèi)模型制備的主控因素。方薇等[11]基于PFC3D離散元軟件，探討了砂雨法成型過(guò)程結(jié)拱現(xiàn)象與孔隙率形成機(jī)制，給出了不同控制因素與孔隙率的關(guān)聯(lián)特征?？傮w而言，現(xiàn)有砂雨法制備技術(shù)研究主要集中于控制因素與影響規(guī)律的物理試驗(yàn)認(rèn)識(shí)，數(shù)值模擬方法的研究尚十分鮮見(jiàn)，且缺少環(huán)境介質(zhì)因素對(duì)砂雨法控制影響規(guī)律的認(rèn)識(shí)和方法可行性的論證。

在前人研究基礎(chǔ)上，以PFC3D離散元數(shù)值軟件為工具，建立砂雨法制備過(guò)程數(shù)值模擬方法，通過(guò)探討落距、篩孔粒徑比、環(huán)境介質(zhì)等因素影響特征與規(guī)律，結(jié)合已有物理試驗(yàn)認(rèn)識(shí)，對(duì)砂雨法制備數(shù)值模擬分析方法的可行性進(jìn)行論證，給出環(huán)境因素對(duì)孔隙率影響基本規(guī)律，為進(jìn)一步研究砂雨法控制因素耦合機(jī)制與綜合參數(shù)選取提供指導(dǎo)借鑒。

1 砂雨法原理

砂雨法制備粗粒土模型的基本原理：利用土體的自由下落與重力勢(shì)能-動(dòng)能轉(zhuǎn)化，模擬自然界風(fēng)力、水力等作用下土體顆粒的碰撞與沉積過(guò)程，通過(guò)落距、移速、流量、環(huán)境介質(zhì)等因素調(diào)整流量強(qiáng)度與動(dòng)能，從而獲得一定孔隙率的試樣。砂雨法自身源于室內(nèi)單元試驗(yàn)重塑試樣，如：三軸、共振柱、直剪等試驗(yàn)，應(yīng)用于土工離心試驗(yàn)?zāi)Ｐ椭苽鋾r(shí)，由于模型試樣尺寸遠(yuǎn)遠(yuǎn)大于單元試驗(yàn)，通常需采用移動(dòng)式漏斗和交叉“S”型行走路徑實(shí)現(xiàn)模型試樣的制備，但因增加了移動(dòng)速度、反復(fù)堆積等制模因素，加大了孔隙率精度和均勻性的控制難度。

目前國(guó)內(nèi)外各單位砂雨裝置多為自主研發(fā)，裝置種類(lèi)及出砂頭型式各異，出現(xiàn)了單孔式、篩網(wǎng)式、鴨嘴式、短線(xiàn)式等出砂頭，但其流量相對(duì)于模型尺寸而言，尚屬于“點(diǎn)”式砂雨法，且孔隙率和均勻性的控制因素及影響規(guī)律尚未形成統(tǒng)一認(rèn)識(shí)。一般而言，砂雨裝置包含移動(dòng)裝置、盛砂容器、軟管和出砂頭等（圖1），可根據(jù)需要，配置一些其他量測(cè)裝置。砂雨法控制因素一般包含落距、水平移速、移動(dòng)軌跡、出砂頭形狀及尺寸、篩網(wǎng)數(shù)量等。本文將以篩網(wǎng)式出砂頭為例，探討PFC3D數(shù)值模擬方法研究砂雨法控制因素與影響規(guī)律的可行性。

圖1 砂雨法工作原理示意[10]Fig.1 Schematic drawing of pluviation[10]

2 PFC3D 建模方法

2.1 PFC 簡(jiǎn)介

PFC（Particle Flow Code）是一款基于顆粒流模擬的離散元數(shù)值分析軟件，分為PFC3D和PFC2D兩個(gè)模塊，由美國(guó)Itasca 公司開(kāi)發(fā)，主要用于模擬各種顆粒相關(guān)的工程與科學(xué)仿真問(wèn)題。目前，PFC 廣泛應(yīng)用于材料工程、礦業(yè)工程、土木工程、食品工程、制藥工程及農(nóng)業(yè)等領(lǐng)域，取得了較豐富的研究成果，為科學(xué)研究提供了一種有效的數(shù)值方法支撐。

作為一種基于非連續(xù)力學(xué)理論開(kāi)發(fā)的顆粒流分析軟件，PFC 建模的實(shí)體組成一般包含基本形態(tài)顆粒、墻體、多形態(tài)顆粒等；介質(zhì)材料分為散體材料和黏結(jié)材料，包含柔性、剛性顆粒/群體/幾何體等。PFC 數(shù)值運(yùn)算特點(diǎn)為通過(guò)各類(lèi)宏-細(xì)觀(guān)力平衡關(guān)系和功、能量守恒關(guān)系，實(shí)現(xiàn)流體、固體、流-固耦合等動(dòng)力與靜力問(wèn)題的復(fù)雜計(jì)算分析，適用于顆粒集合體流動(dòng)破壞、裂紋發(fā)展、堆積擴(kuò)散、開(kāi)挖填埋等變形過(guò)程及力學(xué)性質(zhì)研究。

2.2 建模方法

采用PFC3D軟件進(jìn)行建模，模擬砂雨法在不同情況下制備粗粒土試樣的成型過(guò)程，分析試樣孔隙率的變化規(guī)律。由于軟件自身一些參數(shù)和屬性設(shè)置限制，考慮本文主要為影響規(guī)律研究，對(duì)土體物理力學(xué)參數(shù)和試驗(yàn)工況作了一定程度簡(jiǎn)化和合理調(diào)整。

采用PFC3D軟件進(jìn)行砂雨法成型過(guò)程建模和孔隙率計(jì)算，模型分為兩組，編號(hào)X1、X2，主要方法和設(shè)置參數(shù)如下：顆粒與漏斗及模型箱面之間的接觸采用線(xiàn)性接觸模型，其中法向剛度為100 MPa，切向剛度為100 MPa，摩擦系數(shù)為0.1。顆粒與顆粒之間的接觸采用赫茲接觸模型，顆粒剪切模量為20 MPa，泊松比為0.2，摩擦系數(shù)為0.35。顆粒密度為2 600 kg/m3，初始生成顆?？紫堵蕿?.45。X1 組模擬粒徑范圍為1.8～2.2 cm，主要研究高度與篩孔粒徑比對(duì)孔隙率的影響，其對(duì)應(yīng)漏斗長(zhǎng)、寬、高尺寸為200 mm×200 mm×900 mm，模型箱長(zhǎng)、寬、高尺寸為600 mm×200 mm×1 300 mm，落距為40、60、90 和120 cm。X2 組模擬粒徑范圍為1.8～2.2 mm，主要研究環(huán)境介質(zhì)對(duì)孔隙率的影響，其對(duì)應(yīng)漏斗長(zhǎng)、寬、高尺寸為24 mm×24 mm×60 mm，模型箱長(zhǎng)、寬、高尺寸為60 mm×24 mm×1 000 mm，總落距為30、50、70 和90 cm，水層厚度為0 和20 cm。模型箱長(zhǎng)度大于漏斗長(zhǎng)度，使顆粒下落有一定的自由擴(kuò)散空間，模型箱寬度與漏斗相同，模型箱高度只需高于漏斗下端高度即可。其他參數(shù)：重力加速度為9.81 m/s2，空氣密度為1.29 kg/m3，黏滯系數(shù)取1.79×10?5Pa·s，水的密度為1 000 kg/m3，黏滯系數(shù)取0.001 Pa·s。模擬土顆粒由空氣落入水中工況時(shí)，水面高度恒定設(shè)為20 cm，不考慮水面隨顆粒入水而上升的情況。模擬的篩網(wǎng)網(wǎng)格為3×3、4×4、5×5 三種方形網(wǎng)格，網(wǎng)格尺寸L與平均顆粒半徑r之比為6.7、5.0 和4.0，分別代表稀疏、中密和密集3 種篩網(wǎng)，如圖2 所示。

圖2 不同篩網(wǎng)Fig.2 Different sieves

試樣孔隙率由設(shè)定的多個(gè)不同勘測(cè)球，求解孔隙與勘測(cè)球體積之比的平均值所得。需要說(shuō)明的是，由于模擬顆粒均為圓形球體，在模型箱邊緣處分布較為整齊，對(duì)孔隙率略有影響，因此，勘測(cè)點(diǎn)設(shè)置的測(cè)量球半徑不應(yīng)過(guò)大，宜與模型箱邊界保持一定距離。同時(shí)，選取的測(cè)量球半徑也不應(yīng)過(guò)小，否則可能剛好測(cè)到一些顆粒過(guò)于密集或稀疏的區(qū)域，使數(shù)據(jù)偏差較大，出現(xiàn)較大離散。模擬過(guò)程中孔隙率勘測(cè)點(diǎn)與測(cè)量球的設(shè)計(jì)如圖3 所示，測(cè)量球距離模型邊界有適當(dāng)距離。最終制備的試樣均勻性良好，四周孔隙率略大于中間孔隙率。

圖3 孔隙率測(cè)量球分布Fig.3 Distribution of porosity measuring balls

3 控制要素與影響規(guī)律

基于已有砂雨法物理制備試驗(yàn)的研究成果與認(rèn)識(shí)，對(duì)落距、篩孔粒徑比、環(huán)境介質(zhì)等3 個(gè)因素進(jìn)行數(shù)值模擬研究。

3.1 落 距

落距是各類(lèi)砂雨法制備模型的主要控制因素之一，指模型表面至出砂口的垂直高度。落距40、60、90 和120 cm 條件下，3 組砂雨法制備模擬試驗(yàn)孔隙率的變化如圖4 所示。由圖4 可見(jiàn)，孔隙率與落距呈負(fù)相關(guān)性，當(dāng)落距由40 cm 增大到120 cm，3 組孔隙率分別由37.75%、37.86%、38.48%減小到36.24%、36.65%、37.14%；此外，還可發(fā)現(xiàn)孔隙率的下降趨勢(shì)隨落距增大而逐漸變緩，原因?yàn)槭芸諝庾枇ψ饔?，顆粒落速隨下降高度增加而趨近臨界值，同時(shí)，顆粒撞擊模型表面時(shí)的反彈次數(shù)增加，也進(jìn)一步限制了孔隙率的減小趨勢(shì)。

圖4 不同高度下孔隙率變化曲線(xiàn)Fig.4 Porosity variation at different heights of the simulation

3.2 篩孔粒徑比

篩孔尺寸是篩網(wǎng)式砂雨法的關(guān)鍵控制因素之一，已有經(jīng)驗(yàn)認(rèn)為單孔式砂雨法孔徑與粒徑之比以3～8 為宜。為認(rèn)識(shí)篩孔尺寸與平均粒徑相關(guān)性對(duì)孔隙率的影響規(guī)律，將以篩孔尺寸與平均粒徑之比作為研究對(duì)象，并簡(jiǎn)稱(chēng)為篩孔粒徑比。以L(fǎng)/r≈6.7、5.0 和4.0 分別模擬稀疏、中密和密集篩網(wǎng)3 種工況，PFC3D的模擬分析結(jié)果如圖5 所示。由圖5 可以看出，孔隙率與篩孔粒徑比呈正相關(guān)性，隨篩孔粒徑比增大，孔隙率在不同落距下均表現(xiàn)為一定的增大現(xiàn)象。原因?yàn)楹Y孔粒徑比增大，將導(dǎo)致土顆粒流截面的平均空氣阻力增大，使平均落速和顆粒動(dòng)能減??；同時(shí)，撞擊模型表面時(shí)，顆粒流截面擴(kuò)大有助于加強(qiáng)顆粒之間相互碰撞，從而縮小顆粒跳躍半徑不利于卷吸流體的釋放，形成相對(duì)較大孔隙率。需要注意，篩網(wǎng)式砂雨法制備模型時(shí)，篩孔尺寸選擇不宜過(guò)小，否則容易出現(xiàn)阻塞現(xiàn)象。

與已有研究成果[8-9]相比，圖4 和圖5 中PFC3D數(shù)值模擬方法獲得了與現(xiàn)有物理試驗(yàn)認(rèn)識(shí)十分吻合的分析結(jié)果，并且獲得的曲線(xiàn)更加光滑、規(guī)律性更強(qiáng)，一定程度證明了PFC3D數(shù)值仿真方法模擬砂雨法成型過(guò)程和探討其控制因素影響規(guī)律的可行性。

圖5 不同篩孔粒徑比下孔隙率曲線(xiàn)Fig.5 Porosity variation of different mesh-grain size ratios of the simulation

3.3 環(huán)境介質(zhì)

理論上而言，環(huán)境介質(zhì)對(duì)砂雨法制備模型的孔隙率具有顯著影響，已有國(guó)際學(xué)者采用真空、水兩種環(huán)境介質(zhì)獲得了土體最小和最大孔隙率，并建議以此方法替代現(xiàn)行國(guó)際土工試驗(yàn)方法標(biāo)準(zhǔn)[4]。但實(shí)際中，關(guān)于環(huán)境介質(zhì)對(duì)砂雨法制備模型孔隙率的影響規(guī)律研究和認(rèn)識(shí)尚且不足。以前期開(kāi)展的水沉式飽和方法和砂雨法制備的松散飽和砂模型試驗(yàn)為樣本[10]，對(duì)比和檢驗(yàn)空氣、水兩種介質(zhì)對(duì)模型孔隙率的影響特征與規(guī)律。

總落距30、50、70 和90 cm，模型表面水層厚度0 和20 cm 條件下，PFC3D的模擬分析結(jié)果如圖6 所示。從圖6 容易發(fā)現(xiàn)，模型表面水層存在前后，孔隙率發(fā)生了較大變化；增加水層時(shí)，孔隙率整體上升約2%，對(duì)應(yīng)相對(duì)密實(shí)度約10%。同時(shí)，在模型表面存在水層時(shí)，隨著空中落距增加，孔隙率未發(fā)生較大變化，表明模型表面水層對(duì)孔隙率發(fā)揮了主導(dǎo)控制作用。這些現(xiàn)象可以解釋為水的密度遠(yuǎn)大于空氣密度，相應(yīng)產(chǎn)生的流體阻力、浮力也遠(yuǎn)高于空氣，而顆粒達(dá)到的臨界速度將遠(yuǎn)低于空氣，從而表現(xiàn)出水層“緩沖降速”機(jī)制。此外，圖6 中PFC3D模擬分析結(jié)果與已有物理制備試驗(yàn)規(guī)律相符，說(shuō)明其能較好地模擬環(huán)境介質(zhì)對(duì)砂雨法制備模型孔隙率的影響規(guī)律。

圖6 考慮環(huán)境介質(zhì)時(shí)孔隙率變化曲線(xiàn)Fig.6 Porosity variation in different environmental media of the simulation

為驗(yàn)證模型表面水層的“緩沖減速”機(jī)制，圖7 給出了總落距相等條件下模型表面水層厚度0 和20 cm 時(shí)，PFC3D模擬砂雨法成型過(guò)程中相同時(shí)刻的土顆粒流下落狀態(tài)，以顆粒的疏密程度判斷顆粒下落速度，可知模型表面無(wú)水時(shí)，顆粒速度隨落距增加而增大，并且達(dá)到一定落距時(shí)，稀疏趨勢(shì)不再明顯；而在模型表面有水層時(shí)，顆粒速度入水之前隨落距增加而增大，但在入水之后，顆粒速度快速減小，顆粒流變得密集，且密集程度很快達(dá)到穩(wěn)定和獲得臨界速度。

圖7 同時(shí)刻顆粒在空氣和水中下落對(duì)比Fig.7 Contrast of particle falling in the air and water at the same time

為進(jìn)一步揭示土體顆粒入水之后的速度變化過(guò)程，借鑒相關(guān)文獻(xiàn)[7,12-14]，采用單個(gè)球體和根據(jù)式（1）力平衡方程，對(duì)顆粒由空氣入水后的速度發(fā)展進(jìn)行了計(jì)算，結(jié)果如圖8 所示。等式（1）中左側(cè)為顆粒下落的擬動(dòng)力，右側(cè)由左起3 項(xiàng)分別為顆粒下落的自然重力、流體浮力和流體阻力。

圖8 顆粒落水后速度變化曲線(xiàn)Fig.8 Diagram of velocity change of particle after falling into water

式中：ad為顆粒加速度；V為顆粒體積；ρs為 顆粒密度；ρf為環(huán)境介質(zhì)密度；A為顆粒迎流面面積；v為顆粒速度。Cd為阻力系數(shù)，Cd=24/Re（Re<1）；Cd=18.5/Re0.6（1≤Re≤1 000）；Cd=0.44（1 000

由圖8 計(jì)算結(jié)果可知，相同總落距90 cm 條件下粒徑1.4、1.8、2.2 和2.6 mm 的顆粒，雖然入水的落速略有不同，但是入水之后，一致出現(xiàn)顆粒速度顯著減小趨勢(shì)，并在水深約5 cm 時(shí)趨于臨界值，使得空中落距對(duì)顆粒速度和孔隙率的影響基本可忽略，而圖6 中PFC3D模擬工況中模型表面水層厚度為20 cm，已足夠讓所有顆粒達(dá)到臨界速度。

需特別說(shuō)明，基于PFC3D軟件進(jìn)行土顆粒流由空氣入水模擬時(shí)，為真實(shí)模擬土體顆粒速度在兩種流體介質(zhì)變換條件下速度的變化規(guī)律，采用自編嵌入式模塊對(duì)兩種流體介質(zhì)變化和顆粒入水過(guò)程進(jìn)行了計(jì)算分析，具體為：以式（1）中兩種流體介質(zhì)條件的力平衡方程為基礎(chǔ)，使用ball.list 命令給所有顆粒增加浮力和流體阻力，并使用callback 命令將其加入顆粒的受力系統(tǒng)，使軟件每個(gè)計(jì)算步都會(huì)對(duì)所有顆粒進(jìn)行重力、浮力、阻力的計(jì)算，以保證計(jì)算的合理性。浮力分為3 種情況：（1）顆粒位于水面以上；（2）顆粒部分落入水中；（3）顆粒完全落入水中。其中顆粒部分落入水中時(shí)，根據(jù)PFC3D軟件內(nèi)置命令ball.radius(ball)、ball.pos.z(ball)命令可以得到每個(gè)顆粒半徑、顆粒的球心高度和水面高度，而建立相應(yīng)力平衡關(guān)系。

4 結(jié)語(yǔ)

（1）基于PFC3D離散元分析軟件，建立了砂雨法成型過(guò)程數(shù)值模擬與孔隙率影響因素分析方法，包含接觸面模型的選取、對(duì)象屬性參數(shù)的設(shè)計(jì)、孔隙率的求解和力平衡方程的建立。通過(guò)探討落距、篩孔粒徑比、環(huán)境介質(zhì)等因素對(duì)孔隙率影響規(guī)律的模擬，獲得與以往物理試驗(yàn)十分吻合的結(jié)果，證明了建立數(shù)值分析方法的可行性。

（2）試樣孔隙率與落距呈負(fù)相關(guān)性，孔隙率隨落距增大而下降，且下降速率隨落距增大而趨勢(shì)變緩。孔隙率與篩孔粒徑比呈正相關(guān)性，孔隙率隨篩孔粒徑比增大而增長(zhǎng)，且增長(zhǎng)速率隨二者比值增大而趨勢(shì)變緩。兩個(gè)因素的影響趨勢(shì)可通過(guò)自由落速、模型表面撞擊次數(shù)及跳躍半徑等變化合理解釋。

（3）以水沉法和砂雨法制備飽和砂土試樣為例，對(duì)比了空氣、水兩種環(huán)境介質(zhì)對(duì)孔隙率的影響特征與規(guī)律，分別從砂雨法成型過(guò)程和單顆粒力平衡方程，解釋了模型表面水層的“緩沖減速”機(jī)制，指出了孔隙率與環(huán)境介質(zhì)密度呈負(fù)相關(guān)性，模型表面水層對(duì)孔隙率起主導(dǎo)控制作用，當(dāng)水層厚度不小于5 cm 時(shí)，試樣孔隙率受落距影響可忽略。