砂雨模型制備PFC3D 的數值模擬

王浩然，王永志，王 海，湯兆光，王體強

(1.中國地震局工程力學研究所 地震工程與工程振動重點試驗室，黑龍江 哈爾濱 150080；2.濟南大學 水利與環境學院，山東 濟南 250022)

砂雨法，也稱落雨法，因能較真實地模擬自然風力、水力等作用下天然粗粒土的沉積過程，從而獲得與天然粗粒土地基較為近似的力學特性；同時，其與插搗法、振動法、擊實法等相比，可有效避免顆粒破碎和獲得較寬密實度范圍，被廣泛應用于室內物理試驗模型制備，尤其是傳統單元試驗和土工離心模型試驗[1-6]。通常，砂雨法利用落距、移速、出砂口尺寸、環境介質、篩網數量等因素控制，實現預期粗粒土模型孔隙率和均勻性的控制目標。Lagioia 等[4]利用模型制備試驗，探討了空氣、水、真空等環境因素對砂雨法制備模型孔隙率的影響，結果表明水及真空環境分別可獲得最大、最小孔隙率。Cresswell 等[5]通過流量、移速和落距等因素影響規律的對比，提出了砂雨法制備最小孔隙率的基本原則與要求。Miura 等[6-7]為認識砂雨法控制因素影響特征，開展了一系列粗粒土制備試驗，指出流量為控制孔隙率的首要因素，而孔隙率的控制精度和均勻性隨孔隙率增大而降低。馬險峰等[8-9]介紹了我國8 所院校開展的砂雨法制備平行試驗，給出了出砂頭類型及孔徑、流量、移速、落距等因素影響基本規律。王海等[10]對比了干砂與飽和砂模型砂雨法制備控制因素的影響差異，指出了兩類模型制備的主控因素。方薇等[11]基于PFC3D離散元軟件，探討了砂雨法成型過程結拱現象與孔隙率形成機制，給出了不同控制因素與孔隙率的關聯特征。總體而言，現有砂雨法制備技術研究主要集中于控制因素與影響規律的物理試驗認識，數值模擬方法的研究尚十分鮮見，且缺少環境介質因素對砂雨法控制影響規律的認識和方法可行性的論證。

在前人研究基礎上，以PFC3D離散元數值軟件為工具，建立砂雨法制備過程數值模擬方法，通過探討落距、篩孔粒徑比、環境介質等因素影響特征與規律，結合已有物理試驗認識，對砂雨法制備數值模擬分析方法的可行性進行論證，給出環境因素對孔隙率影響基本規律，為進一步研究砂雨法控制因素耦合機制與綜合參數選取提供指導借鑒。

1 砂雨法原理

砂雨法制備粗粒土模型的基本原理：利用土體的自由下落與重力勢能-動能轉化，模擬自然界風力、水力等作用下土體顆粒的碰撞與沉積過程，通過落距、移速、流量、環境介質等因素調整流量強度與動能，從而獲得一定孔隙率的試樣。砂雨法自身源于室內單元試驗重塑試樣，如：三軸、共振柱、直剪等試驗，應用于土工離心試驗模型制備時，由于模型試樣尺寸遠遠大于單元試驗，通常需采用移動式漏斗和交叉“S”型行走路徑實現模型試樣的制備，但因增加了移動速度、反復堆積等制模因素，加大了孔隙率精度和均勻性的控制難度。

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

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

2 PFC3D 建模方法

2.1 PFC 簡介

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

作為一種基于非連續力學理論開發的顆粒流分析軟件，PFC 建模的實體組成一般包含基本形態顆粒、墻體、多形態顆粒等；介質材料分為散體材料和黏結材料，包含柔性、剛性顆粒/群體/幾何體等。PFC 數值運算特點為通過各類宏-細觀力平衡關系和功、能量守恒關系，實現流體、固體、流-固耦合等動力與靜力問題的復雜計算分析，適用于顆粒集合體流動破壞、裂紋發展、堆積擴散、開挖填埋等變形過程及力學性質研究。

2.2 建模方法

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

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

圖2 不同篩網Fig.2 Different sieves

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

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

3 控制要素與影響規律

基于已有砂雨法物理制備試驗的研究成果與認識，對落距、篩孔粒徑比、環境介質等3 個因素進行數值模擬研究。

3.1 落 距

落距是各類砂雨法制備模型的主要控制因素之一，指模型表面至出砂口的垂直高度。落距40、60、90 和120 cm 條件下，3 組砂雨法制備模擬試驗孔隙率的變化如圖4 所示。由圖4 可見，孔隙率與落距呈負相關性，當落距由40 cm 增大到120 cm，3 組孔隙率分別由37.75%、37.86%、38.48%減小到36.24%、36.65%、37.14%；此外，還可發現孔隙率的下降趨勢隨落距增大而逐漸變緩，原因為受空氣阻力作用，顆粒落速隨下降高度增加而趨近臨界值，同時，顆粒撞擊模型表面時的反彈次數增加，也進一步限制了孔隙率的減小趨勢。

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

3.2 篩孔粒徑比

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

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

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

3.3 環境介質

理論上而言，環境介質對砂雨法制備模型的孔隙率具有顯著影響，已有國際學者采用真空、水兩種環境介質獲得了土體最小和最大孔隙率，并建議以此方法替代現行國際土工試驗方法標準[4]。但實際中，關于環境介質對砂雨法制備模型孔隙率的影響規律研究和認識尚且不足。以前期開展的水沉式飽和方法和砂雨法制備的松散飽和砂模型試驗為樣本[10]，對比和檢驗空氣、水兩種介質對模型孔隙率的影響特征與規律。

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

圖6 考慮環境介質時孔隙率變化曲線Fig.6 Porosity variation in different environmental media of the simulation

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

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

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

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

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

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

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

4 結語

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

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

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