饒登宇，白 冰

（北京交通大學 土木建筑工程學院，北京 100044）

1 研究背景

流體在多孔介質(zhì)中流動的現(xiàn)象廣泛存在于工業(yè)制造、能源開發(fā)、農(nóng)業(yè)生產(chǎn)、環(huán)境治理等各個方面［1］。認識多孔介質(zhì)滲流的溶質(zhì)遷移擴散規(guī)律，將為水土污染治理、垃圾填埋場污染評估、核廢料處置庫的安全性評估等環(huán)境巖土工程問題，提供新的理論依據(jù)和解決辦法。如果能夠在孔隙尺度下從物理本質(zhì)出發(fā)模擬多孔介質(zhì)中溶質(zhì)的運移彌散過程，有助于理解彌散現(xiàn)象產(chǎn)生的物理機制，厘清多孔介質(zhì)彌散度的影響因素。

溶質(zhì)在土體中的遷移主要包括對流和水動力彌散兩個過程［1］。溶質(zhì)隨平均流速運移的過程稱為對流運移；水動力彌散包括分子熱運動引起的分子擴散，以及由于多孔介質(zhì)中流體非均勻運動而產(chǎn)生的機械彌散，機械彌散效果與彌散度和流速相關(guān)。機械彌散現(xiàn)象的出現(xiàn)從微觀尺度上可歸結(jié)為：（1）同一孔隙內(nèi)，地下水質(zhì)點實際流速不等于平均流速u；（2）不同孔隙中孔徑不同，地下水質(zhì)點的實際流速不同，因而不同孔道內(nèi)的流體質(zhì)點出現(xiàn)速度差；（3）在直線距離相同時，多孔介質(zhì)中有的流體質(zhì)點會比其他質(zhì)點沿著更長的路徑運動。若不考慮彌散度的尺寸效應，則影響機械彌散的兩個主要因素可概括為流體質(zhì)點的速度差以及迂曲路徑差。本文將就這兩種因素對彌散的影響展開研究。

受限于實驗尺度和孔隙結(jié)構(gòu)的復雜性，對于孔隙介質(zhì)的溶質(zhì)遷移問題，傳統(tǒng)的研究方法一般在表征體元（REV）尺度下采用平均化的狀態(tài)方程加以概括，認為流體在整個多孔介質(zhì)的空間內(nèi)滲流，并不關(guān)注孔隙內(nèi)的實際流動情況，對溶質(zhì)彌散現(xiàn)象的描述停留在擴散系數(shù)、平均流速等若干宏觀參數(shù)［2］。近年來，孔隙尺度下溶質(zhì)運移的相關(guān)研究有了新的思路：一方面，一些無損測量技術(shù)（如計算機斷層掃描（CT）技術(shù)，X 射線微觀成像等）可直觀得到多孔介質(zhì)中的電導率、滲透率以及其他微觀結(jié)構(gòu)特征，已經(jīng)逐步被用于流體在多孔介質(zhì)內(nèi)的分布、流動、擴散等研究中［3-4］；另一方面各種新型計算流體力學方法（如SPH（光滑粒子流體動力學）、LBM（格子玻爾茲曼方法）等）不斷發(fā)展成熟，使得采用數(shù)值方法從細觀尺度下研究滲流問題成為了可能。數(shù)值方法與實驗法相比，具有周期短、耗費小且不易受外界環(huán)境干擾等優(yōu)勢；與解析法相比，它既可以求解復雜微分方程，又可進行機理研究，并對理論模型進行精確驗證。

本文采用的SPH 方法是一種拉格朗日形式的無網(wǎng)格粒子算法，最初被提出是用于解決天體物理學問題，現(xiàn)已被廣泛地應用于流體力學。Monaghan［5］將其首先延伸至水動力學的模擬，后續(xù)學者進一步采用SPH 方法計算潰壩、水下爆炸、流體與建筑物的相互作用以及孔隙流等問題［6-7］。Morris［8］模擬了低雷諾數(shù)下的不可壓縮流，并將SPH 方法的計算結(jié)果與有限元法進行了對比； Zhu［9］建立孔隙尺度下周期多孔介質(zhì)滲流計算的SPH 模型，用于求解泰勒彌散問題。Tartakovsky［10］提出多相流模型，拓展了SPH 法在多孔介質(zhì)和斷裂介質(zhì)中模擬多相流和反應性遷移的應用。Ryan［11］基于SPH 方法建立界面反應模型，模擬溶質(zhì)的競爭吸附問題，并將孔隙尺度下的模擬結(jié)果與達西尺度模型進行對比驗證。另一種可用于介觀尺度流體力學研究的LBM 方法也受到較多關(guān)注。Yang［12］基于LBM 方法模擬討論了多孔結(jié)構(gòu)隨機性、粒徑分布和顆粒形狀對滲透率的影響，但并未探究多孔結(jié)構(gòu)特性對溶質(zhì)運移及彌散度的影響。需要指出的是，LBM 方法仍然基于網(wǎng)格，在多尺度問題中，如果計算區(qū)域復雜（如大畸變、存在物質(zhì)交界面或自由表面等）易產(chǎn)生網(wǎng)格扭曲等問題［2］，從而導致數(shù)值誤差過大或計算不收斂。在處理動量驅(qū)動流及多孔介質(zhì)等復雜幾何圖形時，SPH 比基于網(wǎng)格的方法有天然優(yōu)勢。首先，線性動量守恒方程的SPH 離散化格式滿足伽利略不變性［11］；其次，SPH 的拉格朗日框架能夠便捷生成可移動和變形的邊界，而不需要復雜的跟蹤算法；并且允許在界面上實現(xiàn)復雜的化學或物理作用，比如擴散和吸附作用。

本文基于SPH 方法，從描述流體質(zhì)點運動和溶質(zhì)擴散的物理本質(zhì)出發(fā)，在孔隙尺度模擬多孔介質(zhì)中溶質(zhì)運移的彌散現(xiàn)象，分析流體質(zhì)點的速度差和迂曲路徑差對機械彌散效應的影響，并討論彌散度與多孔介質(zhì)結(jié)構(gòu)參數(shù)的關(guān)系。

2 SPH 方法基本原理

在SPH 方法中，問題域的狀態(tài)用一系列粒子來描述，粒子之間不需要任何連接，粒子兼具近似點和材料成分的功能。利用SPH 求解，首先需用積分表示法來近似場函數(shù)（場函數(shù)核近似法），然后將其離散為級數(shù)求和的形式（粒子近似法），則f（x）的積分表達式可寫成［10，13］：

式中：W (x-x′，h)為光滑函數(shù)，或稱核函數(shù)，形態(tài)如圖1中所示，積分表達式中的核函數(shù)項可理解為點x′處的函數(shù)值在x 處的權(quán)重，h 為核函數(shù)的光滑長度；Ω為積分區(qū)域。圖1中κh 為核函數(shù)的影響域半徑。

將式（1）轉(zhuǎn)化為支持域內(nèi)所有粒子疊加求和的離散化形式，則第i 粒子點對應的粒子近似式為［13-14］：

圖1 粒子近似圖

式中：xi和xj分別為粒子點i 和j 對應的坐標，粒子間距xij=xi-xj；相應的mj和ρj為粒子點的質(zhì)量和密度；以下將W (xij，h)簡寫為Wij。

本文的核函數(shù)選用一種常見的三次B-樣條函數(shù)［6］：

3 SPH 水動力擴散模型

3.1 質(zhì)量方程對密度直接使用SPH 近似法，粒子i 密度可通過其支持域內(nèi)所有粒子的密度加權(quán)平均近似得到［6，8］：

式中：n 為粒子i 的支持域內(nèi)粒子總數(shù)；mj為粒子的質(zhì)量；Wij為粒子j 對粒子i 產(chǎn)生影響的光滑函數(shù)。

Randles［15］對式（4）進行了標準化改進：

式（5）的改進可以提高自由邊界處和相同材料粒子密度不連續(xù)交界面處的精度。

3.2 動量方程對流體微元，動量方程（不考慮重力）可寫為：

式中：u 為流速；等式右邊兩項分別表示單位質(zhì)量上的各向同性壓力差產(chǎn)生的力和黏性應力差產(chǎn)生的加速度，其中：

式中：上標a 和b 表示坐標方向；p 為壓強；δ為狄拉克函數(shù)；τ為黏性剪應力；μ為黏性系數(shù)；ε為剪應變率。

綜上，SPH 離散化形式可寫為［6］：

3.3 溶質(zhì)擴散方程在孔隙尺度下，流體質(zhì)點間的溶質(zhì)交換僅表現(xiàn)為分子擴散，采用經(jīng)典的Fick 擴散方程描述，其控制方程為［1］：

式中De為分子擴散系數(shù)。

將彌散系數(shù)考慮為常數(shù)，參考熱傳導方程［16-17］的離散形式，離散化式（10），得到在粒子點xi處的溶質(zhì)遷移二維近似表達式：

3.4 人工壓縮性對于不可壓縮流體，采用實際的狀態(tài)方程會導致計算時間步長過小，因此理論上不可壓的流體在數(shù)值計算時一般考慮為微可壓縮流體。這里所采用的人工壓縮性狀態(tài)方程為［8］：

式中：c 為人工速度，一般取最大流速的10 倍。

為保證粒子位置的相對齊整，采用Monaghan［18］推導的“XSPH”方法平均化粒子的相對速度：

式中γ 為常數(shù)，常取0.3。XSPH 方法可有效保持粒子的運動速度與相鄰粒子的平均速度相近，從而避免相近粒子的相互穿透。

3.5 固液邊界處理SPH 計算中，常用的邊界處理方法有：邊界力法、鏡像粒子法、耦合邊界法等。為處理多動介質(zhì)中復雜的固液邊界，本文采用的邊界力模型為描述粒子間相互作用的Lennard-Jones 力模型［2，5］：

式中：r0為截止距離，流體質(zhì)點只有運動到這一距離以內(nèi)，才會受到邊界力的作用；n1、n2為常數(shù)，通常取12 和4；F 為經(jīng)驗參數(shù)，根據(jù)具體問題確定，一般取為與速度最大值的平方相等的量級［2］。

3.6 時間步為保證數(shù)值計算的穩(wěn)定性，時間步長需滿足CFL 條件（Courant-Friedrichs-Levy condition）。對SPH 的水動力模型，需滿足下列規(guī)則［8，10］：

式中：h 為光滑長度；c 為人工速度；fa為粒子加速度，ν=μ/ρ為運動黏度。

而對于擴散問題，時間步長則需滿足［16-17］：

式中：β為計算參數(shù)，相關(guān)研究表明當β＜0.15 時計算是穩(wěn)定的，一般可取β=0.1。計算中，密度、速度和位置時間積分采用具有二階精度的跳蛙法推進［10］，濃度的時間積分則采用Euler 公式［17］。

4 仿真實驗方法

首先介紹無量綱數(shù)Peclet，Peclet 數(shù)表示對流與擴散的相對比例［19］，常用P 或Pe 表示：

其中：u 為滲流速度，m/s；R 為特征長度，m；De為分子擴散系數(shù)，m2/s。Pe 數(shù)增大，輸運中分子擴散輸運比例下降，對流輸運增大。當Pe 小于1 時，分子擴散占優(yōu)；而當Pe 超過100 時，可基本忽略分子擴散的影響［20］。

4.1 模型處理對計算模型做如下規(guī)定：（1）流場內(nèi)飽和；（2）多孔介質(zhì)骨架不可壓縮且不發(fā)生移動；（3）假定流體在多孔介質(zhì)中流速穩(wěn)定，不考慮重力對流動的直接影響。

為使邊界條件實現(xiàn)速度控制和溶質(zhì)濃度控制，將進入流場前后（圖2中L1 至L2，L3 至L4）范圍內(nèi)視為邊界區(qū)域（邊界區(qū)域尺寸的選取需保證流場內(nèi)流體粒子的影響域完整），進入該區(qū)域內(nèi)的流體粒子在水平方向的速度值設定為注入速度、濃度值設定為注入濃度，其他性質(zhì)與正常流體粒子相同。為實現(xiàn)下游流出粒子的重復利用和溶液的持續(xù)注入，設計了粒子傳輸模塊：統(tǒng)計越過L3 界面的粒子個數(shù)記為Nout，以及粒子攜帶的溶質(zhì)濃度Cout，并將流出L4 界面的粒子初始化后遷移至注入端，以此實現(xiàn)模擬流體持續(xù)流動。流場的上下界面均為周期邊界［9］，越過下邊界的粒子會進入上邊界，即可視作流場沿周期邊界的正交方向重復延伸。

圖2 二維多孔介質(zhì)流場的邊界處理示例

定義孔隙水體積數(shù)pv（Pore water volume），它是無量綱時間變量，反映流出水體積與總孔隙水體積的比值［21］，在宏觀尺度下，pv與時間t 有以下關(guān)系：

式中：u 為注入速度，m/s；A 為橫截面積，m；Vp為孔隙體積。在多孔介質(zhì)中，Vp=φV，φ為孔隙率，V 為總體積，m3；而在規(guī)則流場中有，pv=ut/L，L 為流距，m。

在這里，pv即為流出介質(zhì)區(qū)域的粒子總數(shù)Nout與流場內(nèi)總流體粒子數(shù)Ntotal的比值：

4.2 計算過程首先需將計算域離散為粒子點，并賦初值，包括坐標、質(zhì)量、密度、光滑長度等基本參數(shù)，以及設定粒子類別（如圖2中的流場內(nèi)流體粒子、兩側(cè)邊界區(qū)域內(nèi)流體粒子以及固相粒子）。在每一時間步中，先采用N-S 方程的SPH 離散形式（式（5）、式（8）、式（11）、式（12））計算得到所有粒子的速度、加速度、坐標、密度、壓力等，對進入固相粒子影響域范圍內(nèi)的流體粒子還需計算罰力（式（13）），再由溶質(zhì)擴散方程（式（10）），計算得到粒子點的濃度值。最后根據(jù)粒子的坐標，判定是否進入兩側(cè)的邊界區(qū)域內(nèi)，更新粒子類別；若粒子位置超出最右側(cè)邊界（L4），則將其參數(shù)初始化（密度、壓強恢復初始值）后遷移到左側(cè)注入端（L1）。對于兩側(cè)邊界區(qū)域內(nèi)的粒子，需規(guī)定水平速度等于注入速度，左側(cè)注入?yún)^(qū)（圖2中L1 至L2）內(nèi)粒子濃度等于注入濃度，再進入下一時間步。時間步的推進方式如3.6 節(jié)所述。

4.3 解析解檢驗若將圖2中的流場區(qū)域全部填充為流體粒子，則該算例可簡化為一維定解問題?；谄骄魉俚囊痪S多孔介質(zhì)對流-彌散方程為［1，21］：

式中：D 是水動力彌散系數(shù)，m2/s； u 為平均流速，m/s。其中D=De+Dm，De和Dm分別為分子擴散系數(shù)和機械彌散系數(shù)。一般認為［21］，Dm=αuλ，α為機械彌散度，m，簡稱彌散度，參數(shù)λ一般取1［22］。

假定多孔介質(zhì)初始濃度為0，注入濃度為C0的溶液，邊界條件為：C（0，t）=C0，C（∞，t）=0，則該一維注入問題的解析解為：

表1 計算參數(shù)

式中：erfc 函數(shù)為互補誤差函數(shù)，在Matlab、Excel 等計算軟件中均有內(nèi)置函數(shù)可直接調(diào)用。

相關(guān)計算參數(shù)見表1。

此時Pe 數(shù)較小，可不考慮機械彌散系數(shù)，解析解的水動力彌散系數(shù)D 按照分子擴散系數(shù)De計算。經(jīng)過仿真實驗，得到了不同分子擴散系數(shù)時流場末端x=L 處的穿透曲線，并與解析解進行對照，結(jié)果見圖3。對比發(fā)現(xiàn)，SPH 模擬計算結(jié)果與解析解吻合較好，說明該SPH 水動力擴散模型可以用于模擬溶質(zhì)的對流擴散現(xiàn)象。

圖3 流場末端穿透曲線與解析解（x=L）

5 高Peclet 數(shù)流場模擬

5.1 多孔介質(zhì)孔隙流的模擬計算中采用的多孔介質(zhì)為堆積圓形砂礫土，假定為理想的交錯堆積方式。選取某一理想正方形截面作為研究對象，將其置入單向注入的流場中。

模型處理方法如4.1 節(jié)所述，左右兩側(cè)為速度控制的邊界區(qū)域，上下均為周期邊界條件（如圖2），流出右側(cè)界面的粒子將被初始化后遷移至注入端，該算例可視為恒定流速黏性流體穿透二維多孔介質(zhì)薄層的仿真實驗。

圖2中選用的填充顆粒直徑在1.0～2.5 mm 之間，介質(zhì)域為0.012 m×0.012 m 的方形區(qū)域，注入流速恒定u0=0.04 m/s，孔隙率φ=0.58，粒子初始間距為0.0002 m，時間步間隔Δt=0.00002 s。計算中固相粒子密度按2 g/cm3計算，位置固定，其加速度、速度以及壓強值始終為0，計算結(jié)果見圖4至圖7。由圖4可知，在x 方向上，流速大的區(qū)域集中在基質(zhì)顆粒之間的孔道內(nèi)，孔道內(nèi)流速分布具有拋物線形泊肅葉流特征，不同粗細的孔隙通道內(nèi)質(zhì)點流速有明顯差異。

圖5分 別 為 注 入0.02s、0.05s、0.1s 以 及0.2s 后攜帶溶質(zhì)的流體質(zhì)點在多孔介質(zhì)中的彌散現(xiàn)象。若追蹤注入界面上分布的若干顆粒子的運動軌跡，可得到圖6所示的流體粒子跡線圖（圖中將越過周期邊界的粒子畫在了邊界外側(cè)）。根據(jù)該跡線圖，可計算得到該多孔介質(zhì)的迂曲度［9］：

圖4 速度矢量與等值點圖

圖6 流體粒子跡線圖

若溶質(zhì)分子擴散系數(shù)De= 1.35×10-9m2/s（取為NaCl 在20 ℃水中的擴散系數(shù)），在多孔介質(zhì)中計算Pe 數(shù)的特征長度一般取為顆粒粒徑［9］，可得Pe 數(shù)3×104～7.4×104?100，故忽略分子擴散的影響。則水動力彌散系數(shù)D=Dm=αu，彌散度α的大小可采用式（20）對穿透曲線進行曲線擬合得到（如圖7）。彌散度與迂曲度、孔隙率的關(guān)系將在5.3 節(jié)統(tǒng)一討論。

圖7 介質(zhì)末端的穿透曲線及彌散度的擬合曲線（R2=0.98）

5.2 速度差對彌散度的影響為進一步研究溶質(zhì)滲透遷移過程中速度差對機械彌散效應的影響，建立了三段理想化的孔隙通道模型，如圖8所示。三段曲折管的管徑、容積相同，以保證三種曲折方式中流體質(zhì)點的迂曲路徑基本一致。仿真實驗中，從左側(cè)注入含某溶質(zhì)的等濃度溶液，邊界區(qū)域粒子速度按泊肅葉流速方程設定［8］：

本算例中，注入和流出端管徑l=1 mm，動力粘度ν=10-6m2/s，單位質(zhì)量的傳動力F=0.8 m/s2，由式（22）可求得相應的中心水平流速u0=0.1 m/s，雷諾數(shù)Re=100，若分子擴散De= 1.35×10-9m2/s，得Pe =7.4×105?100，故只考慮機械彌散。以上三種曲折管在注入恒定流速的黏性流體時，產(chǎn)生了明顯速度分布差異，其穿透曲線如圖9所示。

圖8 溶液注入的機械彌散過程

圖9 三種曲折管流場的穿透曲線圖

由圖9中溶質(zhì)的穿透曲線，可擬合得到表2中彌散度α。根據(jù)表2中彌散度與粒子流速變異系數(shù)的變化趨勢，可知在迂曲路徑相同的情況下，流體質(zhì)點的速度差越大，機械彌散效應越顯著，表現(xiàn)出更大的彌散度。

5.3 彌散度與多孔介質(zhì)結(jié)構(gòu)參數(shù)的關(guān)系彌散度被認為是反映多孔介質(zhì)骨架結(jié)構(gòu)特征的物理量，孔隙的疏密程度以及顆粒的排布方式等均可影響多孔介質(zhì)的彌散效果。常用于描述多孔介質(zhì)的結(jié)構(gòu)參數(shù)包括孔隙率、顆粒不均勻系數(shù)、迂曲度等，但這些參數(shù)與彌散度的聯(lián)系未見有具體的研究。因此，本節(jié)將討論多孔介質(zhì)的孔隙率φ，顆粒不均勻系數(shù)Cu和迂曲度τ與彌散度α的相關(guān)性。這里顆粒的不均勻系數(shù)Cu=d60/d10，反映了大小不同粒組的分布情況，d60和d10通過對顆粒粒徑的插值近似取值。

類似圖2生成一系列具有隨機粒徑分布和不同孔隙率的二維多孔介質(zhì)計算單元，進行了多組多孔介質(zhì)溶質(zhì)穿透仿真實驗，得到的多孔介質(zhì)結(jié)構(gòu)參數(shù)與彌散度的關(guān)系如圖10所示。

速度差的存在，使得攜帶溶質(zhì)的流體質(zhì)點在單位時間內(nèi)的遷移距離存在差異，由此導致了溶質(zhì)彌散現(xiàn)象的產(chǎn)生。隨著孔隙水流速分布的離散性越大，這種彌散作用更加強烈，因而彌散度與速度的變異系數(shù)大致成正相關(guān)（圖10（a））。

表2 三種曲折管彌散度比較

圖10 多孔介質(zhì)結(jié)構(gòu)參數(shù)與彌散度的關(guān)系

在孔隙率、彌散度和迂曲度的關(guān)系中，目前，已有大量的實驗結(jié)果和理論依據(jù)來證明迂曲度隨孔隙率的增大而減小［23-24］（圖10（b））。由圖10（c）（d），隨著孔隙率的增大，彌散度有減小的趨勢；此外隨迂曲度增大，彌散度大致線性增加。這是由于孔隙率降低，將使得基質(zhì)顆粒間的孔隙通道減小，導致流體質(zhì)點的滲透遷移面臨更多阻礙，從而迂曲度增大。而迂曲度越大，孔隙流體迂回曲折的路徑越復雜，彌散效果則越顯著，彌散度越大，因而彌散度的變化趨勢與孔隙率呈負相關(guān)；圖10（e）中，隨著迂曲度的標準差增大，流體質(zhì)點穿透單位長度的多孔介質(zhì)時，所走過的迂曲路徑差異越大，因而彌散度越大。圖10（f）中，隨顆粒不均勻系數(shù)增大，彌散度有增大趨勢。

綜上，根據(jù)多組隨機粒徑的二維多孔介質(zhì)溶質(zhì)穿透仿真實驗結(jié)果，迂曲度、孔隙率、迂曲度標準差以及顆粒不均勻系數(shù)與彌散度均存在明顯關(guān)聯(lián)。迂曲度作為反映多孔介質(zhì)內(nèi)部孔隙通道復雜程度的重要參數(shù)，可通過實驗手段測量［23］；且迂曲度與多孔介質(zhì)彌散度相關(guān)性較好，今后可進一步研究采用迂曲度定量估算彌散度的可能性。

6 結(jié)論

本文基于SPH 方法，從溶質(zhì)遷移擴散的物理本質(zhì)出發(fā)，實現(xiàn)了孔隙尺度下多孔介質(zhì)溶質(zhì)彌散的仿真實驗。仿真實驗設計中，充分利用了SPH 方法的拉格朗日粒子特性，可便捷得到仿真實驗中的質(zhì)點速度分布、迂曲路徑以及溶質(zhì)穿透曲線。此外，采用周期邊界和上下游粒子傳輸模塊可實現(xiàn)介質(zhì)域的重復拓展和溶液的持續(xù)注入功能。研究中，通過建立三段迂曲路徑一致的孔隙通道模型分析速度差對機械彌散效應的影響。然后在多組溶質(zhì)穿透多孔介質(zhì)的仿真實驗基礎上，討論多孔介質(zhì)結(jié)構(gòu)參數(shù)對彌散度的影響。結(jié)果表明：（1）在控制流體質(zhì)點迂曲路徑相同的情況下，流動速度差對溶質(zhì)機械彌散效應仍有顯著影響，流速的變異系數(shù)越大，則介質(zhì)的彌散度越大；（2）彌散度與迂曲度、迂曲路徑差以及不均勻系數(shù)大致呈正相關(guān)，與孔隙率呈負相關(guān)；（3）鑒于迂曲度與多孔介質(zhì)彌散度相關(guān)性較好，且迂曲度作為反映多孔介質(zhì)內(nèi)部孔隙通道復雜程度的重要參數(shù)可通過實驗手段測量，在今后的研究中可考慮采用迂曲度定量計算彌散度。

由于多孔介質(zhì)本身的復雜性，介質(zhì)彌散度與結(jié)構(gòu)特征參數(shù)的關(guān)系仍需大量研究支持，而在實驗中涵蓋各特征參數(shù)往往十分困難，本文僅提供了一種采用SPH 進行孔隙尺度數(shù)值實驗的思路。此外，目前對彌散度影響因素的認識較不明確，本文研究結(jié)論也可為理論模型的參數(shù)選擇提供參考。