馬玖辰，崔阿鳳，呂林海，楊 杰

(1.天津城建大學(xué) 能源與安全工程學(xué)院，天津 300384；2.天津大學(xué) 中低溫?zé)崮芨咝Ю媒逃恐攸c(diǎn)實(shí)驗(yàn)室，天津 300072；3.天津城建大學(xué) 地?zé)岣咝Ю眉夹g(shù)研究中心，天津 300384)

含水層儲能作為節(jié)能減排、能源轉(zhuǎn)型的關(guān)鍵技術(shù)，在工業(yè)余熱回收利用、電力消納、清潔建筑供熱等領(lǐng)域得到了廣泛的應(yīng)用，為促進(jìn)“碳達(dá)峰、碳中和”目標(biāo)做出重要貢獻(xiàn)。然而，黏土礦物、砂巖顆粒等懸浮顆粒會在水動力或應(yīng)變力的作用下，侵入流道或孔隙空間，堵塞多孔介質(zhì)的孔喉，對地下含水層的滲透性造成不可逆的損害[1]。因此，研究懸浮顆粒在多孔介質(zhì)中的遷移和沉積對含水層儲能技術(shù)[2]、地下水回灌[3]和地下污染物擴(kuò)散[4]等領(lǐng)域具有重要的理論研究與實(shí)踐應(yīng)用價值。目前，國內(nèi)外學(xué)者主要采用理論模型、柱狀層析柱試驗(yàn)、顯微可視化技術(shù)及數(shù)值模擬計(jì)算等方法開展含水層內(nèi)顆粒運(yùn)移研究。

由于具有相同滲透率的多孔介質(zhì)通常存在不同的孔隙結(jié)構(gòu)，因此，國內(nèi)外學(xué)者通常開展層析柱試驗(yàn)確定多孔介質(zhì)的孔隙度、粒徑及孔隙結(jié)構(gòu)，作為在宏觀尺度下建立含水層顆粒運(yùn)移理論模型的基礎(chǔ)。薛傳成等[5]通過對懸浮顆粒的遷移過程進(jìn)行研究，認(rèn)為滲流溶液溫度與pH值、顆粒種類是影響多孔介質(zhì)中懸浮顆粒遷移規(guī)律的重要因素。Chen等[6]分析了地下水滲流過程對于顆粒釋放的影響。張鵬遠(yuǎn)等[7]通過砂柱試驗(yàn)，得到顆粒滯留程度隨著粒徑增大而提高。Bai[8-9]等通過理論與試驗(yàn)研究相結(jié)合，得到顆粒遷移過程隨著溶液離子濃度與含水介質(zhì)孔隙結(jié)構(gòu)的變化而改變。 馬玖辰等[10]結(jié)合電鏡掃描與顆粒表面ζ電位測試，研究了球形與非球形微納米顆粒形貌變化與其釋放、運(yùn)移、沉積過程的內(nèi)在關(guān)聯(lián)。蔣思晨等[11]通過開展不同形狀顆粒在多種溶液離子強(qiáng)度下的滲透試驗(yàn)，將穿透曲線與DLVO理論相結(jié)合，確定顆粒形貌對其遷移特性具有重要影響。雖然，粒子圖像測速（РIV）[12]、粒子陰影圖測速（РSV）[13]、CT掃描[14]等顯微可視化技術(shù)得以在試驗(yàn)中應(yīng)用，但是仍難以揭示流體與顆粒的相互作用機(jī)理。

當(dāng)前，針對懸浮顆粒在多孔介質(zhì)中運(yùn)移、重組過程的研究，國內(nèi)外學(xué)者提出格子Boltzmann方法（lattice boltzmann method，LBM）、計(jì)算流體力學(xué)（CFD）與離散單元法（DEM）耦合計(jì)算方法。格子Boltzmann方法從孔隙尺度研究多孔介質(zhì)中的顆粒滯留和堵塞，多用于納米流體動力學(xué)[15]與非均質(zhì)多相流[16]。盡管LBM方法能很好地表示彎曲顆粒的無滑移邊界條件局部速度分布，但是通常忽略了粒子運(yùn)動對流場的影響[17]。CFD-DEM方法可以在動力場、顆粒運(yùn)動軌跡、應(yīng)力變化等方面對邊坡降雨[18]、泥沙運(yùn)移[19]、顆粒污染[20]等領(lǐng)域開展有效的研究。在DEM計(jì)算方法中，分別針對每個懸浮顆粒的受力形式與運(yùn)移過程建立控制方程，與LBM耦合計(jì)算可以精準(zhǔn)描述孔隙中固相顆粒與液相粒子的運(yùn)移行為[21]。因此，LBM-DEM耦合計(jì)算方法在簡化控制方程、并行計(jì)算性能方面具有顯著優(yōu)勢[22]。國內(nèi)外學(xué)者采用LBM-DEM方法在路基土壤流態(tài)化[23]、固體流變[24]和隧道突泥[25]等方面開展深入研究。Zhou等[26]研究了滲流溶液濃度、流速和pH對多孔介質(zhì)堵塞機(jī)理的耦合影響。Zhou等[27]通過建立不同滲透率的雙通道模型，揭示了顆粒直徑、流速、顆粒體積分?jǐn)?shù)和注入量對非均質(zhì)儲層顆粒截留率和滲透率損傷的影響。

綜上，由于多孔介質(zhì)孔隙結(jié)構(gòu)具有多樣性，宏觀尺度的運(yùn)移模型與層析柱試驗(yàn)難以有效揭示含水層孔隙結(jié)構(gòu)變化的誘導(dǎo)機(jī)制；而數(shù)值模擬計(jì)算存在并行計(jì)算能力弱、收斂性差等問題。同時，采用LBMDEM方法從介觀尺度研究顆粒形貌對其在含水層中的釋放—運(yùn)移—沉積過程影響的研究尚少見報(bào)道。本文基于多孔介質(zhì)滲流溶液與懸浮顆粒運(yùn)移、受力控制方程，采用LBM-DEM耦合方法，將含水層流場劃分成規(guī)則的網(wǎng)格，通過引入流體虛擬粒子，將滲流速度離散到多個方向。利用浸沒移動邊界法（IMB）處理固相顆粒與液相粒子的相互作用，從介觀尺度對懸浮顆粒在多孔介質(zhì)中的運(yùn)動過程開展數(shù)值計(jì)算。搭建滲流層析柱試驗(yàn)系統(tǒng)，預(yù)測含水層水動力場演化過程，同時驗(yàn)證理論模型與求解方法的正確性。本文將介觀尺度數(shù)值計(jì)算與試驗(yàn)研究相結(jié)合，探究顆粒形貌、地下水動力和多孔介質(zhì)尺寸效應(yīng)等因素對含水層滲透性能的影響及其內(nèi)在誘導(dǎo)機(jī)制。

1 理論基礎(chǔ)與控制方程

1.1 含水層滲流控制方程

本文以連續(xù)介質(zhì)運(yùn)移理論為基礎(chǔ)，運(yùn)用連續(xù)性方程（式（1））和Navier-Stokes方程（式（2））[28]來描述滲流溶液運(yùn)動的流場，在Navier-Stokes方程中，引入外部體積力Ff（式（3））體現(xiàn)外部動力和流體相互作用：

式（1）～（4）中：ρ為流體密度，kg/m3；ε為多孔介質(zhì)孔隙率；t為時間，s；v為流體速度矢量，m/s；P為流體壓力，Рa；υ為流體運(yùn)動黏性系數(shù)，m2/s；K為多孔介質(zhì)滲透率；為多孔介質(zhì)顆粒平均直徑，m；G為外力項(xiàng)，N。

1.2 基于LBM的流體控制方程

通過引入虛擬粒子對流體速度進(jìn)行離散，采用LBM模擬含水層中流體粒子的運(yùn)動過程，對格子Boltzmann演化方程（式（5））進(jìn)行求解，確定各網(wǎng)格節(jié)點(diǎn)處粒子的遷移和碰撞過程[29]。

式中，x為位置矢量。

在每個時間步長Δt內(nèi)，x處流體粒子以離散速度矢量ci運(yùn)動到相鄰的節(jié)點(diǎn)上，繼而在該節(jié)點(diǎn)處與其他粒子發(fā)生碰撞，碰撞后更新粒子速度分布函數(shù)fi(x,t)。流體粒子發(fā)生一系列遷移碰撞，逐漸趨向平衡態(tài)分布，平衡態(tài)分布函數(shù)(x,t)（式（6））可由離散的Maxwell分布函數(shù)表示：

式中，wi為i方向的權(quán)重因子，cs為格子聲速。

基于格子單松弛模型（lattice bhatnagar-grosskrook，LBGK），將粒子分布函數(shù)向平衡態(tài)分布函數(shù)的時間松弛簡化為單個BGK算子，用于代替復(fù)雜碰撞過程；考慮到定壓物理邊界條件和流體相互作用力，在演化方程中引入離散力源項(xiàng)Fi：

式中，τ為無量綱松弛時間。

利用演化方程求解虛擬粒子的碰撞和遷移過程后，由速度分布函數(shù)fi(x,t)的第0和第1速度矩求得流體的宏觀參數(shù)，包括局部密度ρ(x,t)、局部流速v(x,t)：

運(yùn)動黏性系數(shù)υ與格子間距Δx、時間步長Δt有關(guān)。通過Chapman-Enskog展開[30]，由演化方程與Navier-Stokes方程可得到運(yùn)動黏性系數(shù)υ與無量綱松弛時間τ的關(guān)系式為：

式中，c為格子速度，即格子間距與時間步長之比。

1.3 基于DEM的顆粒運(yùn)移控制方程

采用基于軟球模型的離散單元法（discrete element method，DEM）模擬固體懸浮顆粒間的相互作用。所有的顆粒運(yùn)動都遵循牛頓定律，允許彼此有一定重合度，將顆粒相互作用簡化為彈簧-阻尼器[31]過程。因此，顆粒運(yùn)動的動力學(xué)方程分為平移和旋轉(zhuǎn)，考慮到固-液相互作用、固-固接觸以及重力等因素，顆粒的平移由法向接觸力Fcn、切向接觸力Fct、拖拽力Fd、流體與顆粒相互作用力Ff,p、重力Fg、浮力Fb等力共同作用；顆粒的旋轉(zhuǎn)由顆粒碰撞力矩Tp和流體流動施加在顆粒上的力矩TH共同作用，如式（11）、（12）所示：

式（11）～（12）中，u為顆粒平移速度，m為顆粒質(zhì)量，IР為轉(zhuǎn)動慣量，ωp為角速度。

圖1為顆粒在孔隙中運(yùn)動的受力示意圖。

圖1 懸浮顆粒受力示意圖Fig.1 Schematic diagrams of the forces acting on the suspended particles

圖2 試驗(yàn)裝置示意圖Fig.2 Schematic diagram of experimental device

根據(jù)Hooke定律可知，顆粒間的法向接觸力與切向接觸力可以通過顆粒間的法向彈簧與切向彈簧來模擬，如式（13）、（14）所示：

當(dāng)區(qū)域Re數(shù)較低且顆粒不靠近通道壁時，流體受到Stokes阻力，此時顆粒受到拖拽力Fd：

在滲流溶液中，顆粒與流體粒子相互碰撞，從而產(chǎn)生的布朗運(yùn)動會影響懸浮顆粒的遷移，導(dǎo)致顆粒受到相互作用力Ff,p：

顆粒在地下含水層中受到的重力Fg和浮力Fb的關(guān)系可由式（17）表示：

式（13）～（17）中：kn和kt為法向和切向接觸剛度；δn和δt為法向重疊量和切向位移；ηn和ηt為法向和切向阻尼系數(shù)；Δun和Δut為法向和切向相對速度；d為當(dāng)量直徑；kB為Boltzmann常數(shù)；T為溫度；α、σ分別為高斯分布函數(shù)的均值和標(biāo)準(zhǔn)差，設(shè)置為0和1；V為顆粒等效體積；ρp為懸浮顆粒密度；g為重力加速度。

顆粒碰撞力矩與轉(zhuǎn)動慣量可用式（18）和（19）計(jì)算：

式（18）～（19）中，R為從顆粒中心指向接觸點(diǎn)的位置向量。

2 試驗(yàn)系統(tǒng)與研究內(nèi)容

2.1 試驗(yàn)裝置

試驗(yàn)采用長300 mm、內(nèi)徑25 mm、壁厚3.5 mm的有機(jī)玻璃層析柱進(jìn)行強(qiáng)制滲流試驗(yàn)。層析柱水平放置，在出、入口處均放置0.01 mm厚、孔徑為0.18 mm的濾網(wǎng)。選用RSC型雙頭蠕動泵將水容器里高純度去離子水以恒定轉(zhuǎn)速注入層析柱；在LSР01-2A型注射泵的驅(qū)動下，將懸浮液以相同的速度注入層析柱；使用精確度為 ±2%的TL2350濁度儀測量滲流液濁度。

將高純度去離子水不斷注入特定濃度懸浮液進(jìn)行稀釋，測量相應(yīng)懸浮液濁度，得到顆粒濃度與濁度的關(guān)系曲線如圖3所示。將試驗(yàn)結(jié)果的標(biāo)定濃度與濁度值進(jìn)行線性回歸分析，其判定系數(shù)R2＞0.99，說明兩者具有高度相關(guān)性。由圖3可知：顆粒的濃度隨濁度的增大而增大；隨著濁度的增加，球形硅微粉與非球形硅微粉的濃度差隨之增加。由于非球形硅微粉以片狀、多面體結(jié)構(gòu)為主，比表面積大于球形硅微粉，對光線的散射作用增強(qiáng)。因此，在相同濃度下，非球形硅微粉懸浮液的濁度遠(yuǎn)高于球形硅微粉懸浮液。

圖3 懸浮顆粒濁度與濃度的關(guān)系曲線Fig.3 Relationship curves between turbidity and concentration of suspended particles

2.2 試驗(yàn)材料

選用人工制備砂樣作為層析柱中填充的多孔介質(zhì)，為了確保填充材料顆粒級配連續(xù)性與礦物構(gòu)成與地下含水層原始砂樣相近，根據(jù)原砂的機(jī)械組成[10]，將等效粒徑75 μm和120 μm的天然石英砂，按照1∶4比例均勻混合。試驗(yàn)前，將選取的石英砂進(jìn)行酸洗，使用去離子水反復(fù)沖洗，烘干備用。采用密度為2.26 g/cm3、中位粒徑為13 μm的球形和非球形硅微粉作為微納米顆粒物質(zhì)[10]。利用場發(fā)射掃描電子顯微鏡（SEM），觀察微納米顆粒形貌特征與粒徑尺寸分布發(fā)現(xiàn)：球形硅微粉形狀統(tǒng)一為圓球狀，比表面積小，具有良好的分散性；非球形硅微粉形貌、幾何尺寸呈不規(guī)律分布，多為多面體狀顆粒，如圖4所示。

2.3 試驗(yàn)方案與結(jié)果分析

采用分段濕填法填充4根層析柱，確保每根層析柱中含水介質(zhì)重度達(dá)到（1.65±0.02） kg/L，反演計(jì)算得到孔隙率ε近似為0.37。分別配置800 mL濃度為0.25 mg/mL的球形與非球形硅微粉懸浮液，在滲流速度為2.5×10-5m/s與5.0×10-5m/s這兩種工況下，開展4組強(qiáng)制滲流試驗(yàn)。首先，關(guān)閉閥門C，打開閥門A、B，采用蠕動泵以固定流速向?qū)游鲋羞B續(xù)通入高純度去離子水，當(dāng)流出液的滲流速度穩(wěn)定，并且無顆粒出現(xiàn)時，認(rèn)為層析柱中含水介質(zhì)的飽水、排氣過程完成。隨即關(guān)閉閥門A、B，打開閥門C，使用注射泵分別以2.5×10-5m/s與5.0×10-5m/s這兩種流速將制備好的球形與非球形硅微粉懸浮液連續(xù)注入層析柱。待懸浮液注入完畢后，關(guān)閉閥門C，打開閥門A、B，繼續(xù)以相同流速無間斷通入高純度去離子水。4組試驗(yàn)運(yùn)行時長為180 min，每組試驗(yàn)均在室溫（22～24 ℃）下進(jìn)行。試驗(yàn)過程中每4 min采用秒表量筒法測試滲出液的流速；同時，收集滲出液，使用濁度儀測量溶液濁度，根據(jù)擬合關(guān)系曲線（圖3）得到滲出液中的顆粒濃度，通過抽濾稱重法測定4組試驗(yàn)在不同時段內(nèi)的滲出液中懸浮顆粒的累計(jì)質(zhì)量。

釋放顆粒濃度與累計(jì)質(zhì)量動態(tài)變化曲線如圖5所示。

由圖5（a）可知：以滲流速度為5.0×10-5m/s向?qū)游鲋⑸淝蛐喂栉⒎蹜腋∫簳r，滲出液顆粒濃度變化曲線呈雙峰狀分布狀態(tài)，當(dāng)試驗(yàn)進(jìn)行到88 min時，釋放顆粒濃度達(dá)到最大值0.379 mg/mL。當(dāng)滲流速度降低至2.5×10-5m/s時，滲出液顆粒濃度無明顯峰值出現(xiàn)，在試驗(yàn)進(jìn)行的52～132 min，釋放顆粒濃度曲線在0.087～0.167 mg/mL范圍內(nèi)呈不規(guī)則波動。向?qū)游鲋鶅?nèi)注入非球形硅微粉懸浮液時，滲出液的顆粒濃度變化規(guī)律與注入球形硅微粉懸浮液時一致。由圖5（b）可知：滲流速度為5.0×10-5m/s時，釋放顆粒濃度最大值降低到0.11 mg/mL，同時峰值時間出現(xiàn)延遲；滲流速度為2.5×10-5m/s時，釋放顆粒濃度曲線在0.019～0.035 mg/mL范圍內(nèi)波動，波動時間增長20 min。

表1為4組層析柱滲流試驗(yàn)結(jié)果。由表1可知，隨著滲流速度的提高，顆粒釋放率增大。以滲流速度為5.0×10-5m/s為例，注射球形與非球形硅微粉懸浮液相比較，微納米顆粒釋放率提高了47.9%。試驗(yàn)結(jié)果表明，顆粒形貌與水動力對顆粒的遷移特性有重要影響。與球形硅微粉相比，非球形硅微粉比表面積更大，沉積在多孔介質(zhì)表面或被孔喉捕獲的機(jī)率更高。由于滲流剪切應(yīng)力與顆粒法向截面面積成正比[10]，因此，在相同流速下，沉積于多孔介質(zhì)表面的非球形硅微粉難以脫離釋放。

表1 不同形貌顆粒釋放率與滲出液顆粒濃度Tab.1 Release rates of particles and particle concentration of the seepage solution with different morphology

3 數(shù)值計(jì)算與分析

3.1 求解方法

滲流溶液流體粒子、流體粒子與懸浮顆粒及懸浮顆粒之間的受力模式都是引起多孔介質(zhì)內(nèi)懸浮顆粒遷移狀態(tài)改變的誘導(dǎo)機(jī)制。因此，采用歐拉坐標(biāo)下的LBM模擬計(jì)算流體粒子流動過程[16]；采用拉格朗日坐標(biāo)下的DEM軟球模型模擬計(jì)算顆粒遷移過程[32]；同時引入浸沒移動邊界方法（immersed moving boundary method，IMB）處理流體粒子-固體顆粒邊界，從而克服了LMB與DEM模型的動量不連續(xù)問題，實(shí)現(xiàn)對懸浮顆粒運(yùn)移過程的準(zhǔn)確數(shù)值計(jì)算。

基于LBGK模型提出包含外力項(xiàng)的流體粒子速度場演化方程（式（5）），其中，權(quán)重因子wi和離散速度矢量ci均與網(wǎng)格離散化相關(guān)[31]，取決于使用的格子模型。因此，根據(jù)層析柱實(shí)際滲流過程，選擇2維平面具有9個離散速度的D2Q9模型對演化方程進(jìn)行離散求解。每個方向的離散速度和權(quán)重因子由式（20）、（21）確定：

由于IMB方法可以在離散格子中同時包含流體粒子和固體顆粒，因此，在LBM-DEM耦合求解中，通過無滑移邊界條件來實(shí)現(xiàn)固-液相互作用。在IMB中引入格子固含率εs表征格子節(jié)點(diǎn)被顆粒覆蓋的體積分?jǐn)?shù)，其中，流體內(nèi)部節(jié)點(diǎn)εs為0，顆粒內(nèi)部節(jié)點(diǎn)εs為1；當(dāng)流體與顆粒邊界碰撞時，可以確定每個節(jié)點(diǎn)的εs取值（0～1）。根據(jù)IMB計(jì)算原理，將流體粒子速度場演化方程（式（5））進(jìn)行優(yōu)化修正（式（22）），引入碰撞項(xiàng)Ωis（式（23））用以實(shí)現(xiàn)固體顆粒對流體粒子流動的影響：

式（22）～（23）中：±i表示正（反）方向；B為Ωis的加權(quán)函數(shù)，由式（24）確定：

基于優(yōu)化的速度場演化方程（式（22）），結(jié)合流體粒子相互作用的離散力、流體粒子對于懸浮顆粒的拖拽力、顆粒間相互作用力，確定流體粒子施加于固體顆粒上的力FH，即對附加碰撞項(xiàng)在流體邊界、顆粒邊界及顆粒內(nèi)部節(jié)點(diǎn)的所有方向所引起的動量變化求和，流體流動施加在顆粒上的力FH和力矩TH可表達(dá)為：

式中（25）～（26），n為該懸浮顆粒所標(biāo)記節(jié)點(diǎn)的數(shù)量，(xk-xp)為從顆粒中心指向第k個節(jié)點(diǎn)的位置向量。

3.2 LBM-DEM計(jì)算平臺與驗(yàn)證

根據(jù)強(qiáng)制滲流試驗(yàn)系統(tǒng)中層析柱的幾何結(jié)構(gòu)，將計(jì)算模型設(shè)置為300 mm×25 mm的長方形，如圖6所示。根據(jù)所確定的格子距離，將計(jì)算區(qū)域離散為1.2×1011個節(jié)點(diǎn)（1 200 000×100 000）。將計(jì)算區(qū)域的上下邊緣均設(shè)置為封閉的實(shí)體邊界，采用反彈邊界格式處理；滲流溶液與懸浮顆粒從左側(cè)進(jìn)入、右側(cè)流出，因此左、右側(cè)為速度邊界，采用周期邊界格式處理。根據(jù)試驗(yàn)結(jié)果，通過隨機(jī)四參數(shù)法在區(qū)域內(nèi)生成多孔介質(zhì)，確定含水介質(zhì)的孔隙率；同時分別設(shè)置7.68×107個中位粒徑為13 μm的球形和非球形顆粒作為注入含水介質(zhì)中的懸浮顆粒。模擬計(jì)算的基本參數(shù)設(shè)置見表2。

表2 數(shù)值模擬基本參數(shù)Tab.2 Basic parameters of the numerical simulation

圖6 2維數(shù)值計(jì)算模型Fig.6 Two-dimensional numerical model

根據(jù)所建立的演化方程與求解方法，在數(shù)值計(jì)算平臺MATLAB上，開發(fā)了LBM-IMB-DEM耦合求解程序。設(shè)定流體粒子與固體顆粒的初始條件，在每個時間步長上，利用LBM更新流場變化，利用DEM循環(huán)計(jì)算顆粒的碰撞。通過計(jì)算流體粒子作用在固體顆粒上的力和力矩，更新顆粒位置信息，完成流體與顆粒的耦合，耦合計(jì)算流程如圖7所示。由于在LBM計(jì)算中，網(wǎng)格劃分形式與數(shù)量對計(jì)算結(jié)果影響很大，因此通過反復(fù)的預(yù)處理計(jì)算確定格子距離在5×10-7～1×10-6m內(nèi)，與孔隙率及粒子作用力無關(guān)。由于LBM與DEM均選用顯式循環(huán)算法，并且具有各自的時間步長，因此在求解過程中取Δt/ΔtD=100，將DEM作為LBM的子循環(huán)從而確保耦合計(jì)算時間的統(tǒng)一。

圖7 LBM-IMB-DEM耦合流程圖Fig.7 Flowchart of the LBM-IMB-DEM coupling scheme

利用LBM-IMB-DEM耦合求解程序，計(jì)算得到4組試驗(yàn)下層析柱出口處滲流速度vout逐時變化過程。選擇滲流速度vout的絕對誤差與均方根誤差作為計(jì)算值與試驗(yàn)結(jié)果的相似度判定指標(biāo)。如圖8所示，4組試驗(yàn)下模擬結(jié)果在計(jì)算周期中的動態(tài)變化與試驗(yàn)數(shù)據(jù)基本一致，均可全程跟蹤試驗(yàn)過程，vout的均方根誤差均小于10%，最大標(biāo)準(zhǔn)誤差為0.26。為進(jìn)一步驗(yàn)證模型的正確性，利用本文所寫程序?qū)hou[26]和Feng[33]等中的DKT算例進(jìn)行模擬，圖9為不同時刻顆粒流速分布對比，結(jié)果顯示三者的模擬結(jié)果吻合度較高。因此，認(rèn)為所建立的演化方程與耦合求解方法可以準(zhǔn)確地反映懸浮顆粒在含水介質(zhì)中運(yùn)移、沉積及再釋放的重組過程。

圖8 層析柱出口滲流速度動態(tài)變化曲線Fig.8 Curves of the seepage velocity determined at chromatographic column outlet

圖9 DKT算例不同時刻的顆粒流速分布Fig.9 Velocity of the particles at different for DKT case

3.3 計(jì)算結(jié)果與分析

圖10～11表示在不同時刻下，層析柱含水介質(zhì)中滲流速度的演化過程。由于在LBM的計(jì)算過程中需要對流場初始化，將流體粒子節(jié)點(diǎn)均設(shè)置為初始滲流速度，因此在初始狀態(tài)（0 min）下含水介質(zhì)流體粒子的滲流速度與入口邊界處賦值相同。隨著懸浮顆粒溶液的注入，LBM-IMB-DEM耦合求解程序啟動，對比20 min與40 min這兩個時刻，在4組試驗(yàn)中含水介質(zhì)中下游區(qū)域的滲流速度均出現(xiàn)了顯著降低。當(dāng)再次注入去離子水時，滲流速度出現(xiàn)回升，截至計(jì)算周期結(jié)束（180 min），含水介質(zhì)流場的整體分布趨于一致。

圖10 注射球形硅微粉工況含水介質(zhì)流速分布Fig.10 Distribution of liquid velocity in the aquifer medium with injection of spherical silica powders

由圖10可知，向含水介質(zhì)注入球形硅微粉懸浮液時，隨著進(jìn)水流速的提高，滲流速度恢復(fù)程度增加。懸浮顆粒會與多孔介質(zhì)反復(fù)碰撞，即顆粒之間以及流體粒子與顆粒發(fā)生碰撞，造成動能的損失，若在此過程中不能到達(dá)分離點(diǎn)，顆粒最終將沉積在多孔介質(zhì)表面，導(dǎo)致滲流速度下降。提高滲流速度，將使顆粒所受法向接觸力Fcn（式（13））、切向接觸力Fct（式（14））、拖拽力Fd（式（15））增大，使顆?？焖俚竭_(dá)分離點(diǎn)，增強(qiáng)溶液對顆粒的攜帶效果，使更多顆粒流出含水介質(zhì)，流體粒子運(yùn)動空間增加，滲流速度得以恢復(fù)。

由圖11可知，向含水介質(zhì)注入非球形硅微粉懸浮液時，滲流速度降幅增大，僅在高流速工況下的進(jìn)口處出現(xiàn)微弱回升。相同等效體積V下，圓球狀的球形硅微粉比多面體狀的非球形硅微粉的當(dāng)量直徑d小，導(dǎo)致球形硅微粉的轉(zhuǎn)動慣量Ip（式（19））減小，同時顆粒中心到含水介質(zhì)接觸點(diǎn)的距離增加，使顆粒碰撞力矩Tp（式（18））增大，因此球形硅微粉易發(fā)生滾動、脫離。注入非球形硅微粉懸浮液時，由于相互作用力Ff,p（式（16））增大，引起流體粒子與顆粒更劇烈的碰撞，動能損失增大，滲流速度降低。非球形硅微粉比表面積大，易絮凝成團(tuán)，在狹窄孔隙喉道及凹陷處被捕獲的機(jī)率更大，有效孔隙空間的減少增強(qiáng)了對顆粒的阻礙作用，導(dǎo)致顆粒的自由遷移路徑更少，從而引起含水介質(zhì)滲透性下降。

圖11 注射非球形硅微粉工況含水介質(zhì)流速分布Fig.11 Distribution of liquid velocity in the aquifer medium with injection of non spherical silica powders

選擇含水介質(zhì)距進(jìn)水端50、150、250 mm處作為研究對象，計(jì)算確定各斷面全部流體粒子的平均滲流速度vˉ，從而擬合生成動態(tài)變化曲線如圖12所示。

圖12 vin=2.5×10-5 m/s時含水介質(zhì)不同斷面平均滲流速度動態(tài)變化曲線Fig.12 Curves of the average seepage velocity in different sections of the aquifer medium with vin=2.5×10-5 m/s

向含水介質(zhì)注入球形硅微粉懸浮液時（圖12（a）），速度變化曲線呈高斯分布，各斷面平均流速最大降幅分別達(dá)到12.6%、29.2%、38.3%，計(jì)算周期結(jié)束時穩(wěn)定在2.43×10-5、2.36×10-5、2.34×10-5m/s。向含水介質(zhì)注入非球形硅微粉懸浮液時（圖12（b）），速度變化呈玻爾茲曼分布，各斷面平均流速依次下降至2.0×10-5、1.7×10-5、1.3×10-5m/s，并且持續(xù)保持穩(wěn)定。

圖13為進(jìn)水速度提高到5.0×10-5m/s時，含水介質(zhì)不同斷面平均滲流速度動態(tài)變化曲線。如圖13所示，當(dāng)進(jìn)水流速提高到5.0×10-5m/s 時，含水介質(zhì)滲流速度變化均呈現(xiàn)下降、回升、穩(wěn)定3個連續(xù)階段。以距進(jìn)水端150 mm處為例：當(dāng)向含水介質(zhì)注入球形硅微粉懸浮液時，平均流速下降至3.1×10-5m/s，在117 min便回升至4.8×10-5m/s；當(dāng)向含水介質(zhì)注入非球形硅微粉懸浮液時，平均流速最低下降至2.7×10-5m/s，計(jì)算周期結(jié)束時僅回升至3.4×10-5m/s。

圖13 vin=5.0×10-5 m/s時含水介質(zhì)不同斷面平均滲流速度動態(tài)變化曲線Fig.13 Curves of the average seepage velocity in different sections of the aquifer medium with vin=5.0×10-5 m/s

以進(jìn)口流速為5.0×10-5m/s，向含水介質(zhì)注入非球形硅微粉懸浮液為例，3組斷面平均流速最大降幅依次提高，分別為27.1%、43.7%、56.3%；計(jì)算周期結(jié)束時平均流速恢復(fù)率則逐漸降低，依次達(dá)到62.5%、30.4%、28.1%。通過分析可知，當(dāng)懸浮顆粒大量注入含水介質(zhì)時，流體與顆粒的運(yùn)動空間在一定時間內(nèi)迅速減少，隨著顆粒在含水介質(zhì)中的運(yùn)動路程增大，顆粒之間以及顆粒與多孔介質(zhì)的碰撞率提高，動能損失增大，從而引起含水介質(zhì)中、后部平均流速降幅增大，導(dǎo)致懸浮顆粒沉積量增加。

4 結(jié) 論

1）通過引入離散力源項(xiàng)，得到具有外力項(xiàng)的LBGK速度演化方程；根據(jù)DEM的軟球模型，對含水介質(zhì)中懸浮顆粒的平移和旋轉(zhuǎn)過程進(jìn)行受力分析，建立LBMDEM耦合計(jì)算模型。利用IMB引入固含率εs、碰撞項(xiàng)Ωis完成固-液相互作用的動量交換過程，開發(fā)LBMIMB-DEM求解程序。通過層析柱滲流試驗(yàn)，將模擬結(jié)果與試驗(yàn)數(shù)據(jù)相比較，發(fā)現(xiàn)4組試驗(yàn)的vout均方根誤差均小于10%，而且DKT算例與相關(guān)文獻(xiàn)結(jié)果擬合度良好，認(rèn)為耦合計(jì)算模型與求解方法可以準(zhǔn)確地反映懸浮顆粒在含水介質(zhì)中運(yùn)移過程。

2）當(dāng)進(jìn)口流速為5.0×10-5m/s時，將注入含水介質(zhì)的懸浮顆粒溶液由球形硅微粉替換為非球形硅微粉后，所選3組斷面的平均流速最大降幅分別提升15.0%、11.4%、2.6%。研究表明：球形硅微粉轉(zhuǎn)動慣量Ip較小，易發(fā)生滾動、脫離；非球形硅微粉比表面積大，易絮凝成團(tuán)，沉積在多孔介質(zhì)表面后難以脫離釋放。隨著進(jìn)水流速的降低，顆粒形貌的變化對于含水層水動力場演化過程的影響程度增強(qiáng)。

3）在進(jìn)口流速為5.0×10-5m/s，向含水介質(zhì)注入非球形硅微粉懸浮液條件下，當(dāng)截面距進(jìn)口距離由50 mm增加至250 mm時，平均流速最大降幅由27.1%上升至56.3%，恢復(fù)率由62.5%下降至28.1%。隨著含水介質(zhì)計(jì)算區(qū)域的增加，顆粒的運(yùn)動路程增大，顆粒之間以及顆粒與多孔介質(zhì)的碰撞幾率增加，動能損失增大，導(dǎo)致滲流速度降幅增大，恢復(fù)率降低。