李 維 仲， 趙 月 帥， 宋 永 臣

（大連理工大學(xué) 海洋能源利用與節(jié)能教育部重點實驗室，遼寧 大連 116024）

0 引 言

在自然界和很多工業(yè)應(yīng)用中，都包含了多孔介質(zhì)內(nèi)流體流動問題，如石油開采、地下水凈化和流化床等領(lǐng)域.因此，深入了解在不同工況下，流體在多孔介質(zhì)中的流動特性有著重大的實際意義.在實際工程中，流體在多孔介質(zhì)內(nèi)流動問題可以用達(dá)西定律進(jìn)行描述，然而達(dá)西定律是一個基于宏觀觀測的實驗定律，流體和多孔介質(zhì)之間復(fù)雜的相互作用只能體現(xiàn)在一個宏觀的統(tǒng)計參數(shù)k之中.因而，要詳細(xì)研究這個參數(shù)，了解多孔介質(zhì)具體構(gòu)造對k的影響，發(fā)展一種新的能夠在孔隙尺度下直接高效地模擬流體在復(fù)雜多孔介質(zhì)中流動細(xì)節(jié)的數(shù)值模擬方法是非常必要的.

然而，采用傳統(tǒng)的數(shù)值模擬方法，如有限元法、有限體積法等，模擬多孔介質(zhì)內(nèi)流體流動問題比較困難，在處理復(fù)雜的孔隙幾何結(jié)構(gòu)、多尺度、多相流動等問題時，都需要額外的處理技巧，同時容易出現(xiàn)精確度不高和數(shù)值穩(wěn)定性較差等問題.值得一提的是，近年來發(fā)展非?？斓母褡覤oltzmann方法（LBM）也常被應(yīng)用于孔隙尺度下多孔介質(zhì)內(nèi)流動問題的研究.LBM是基于玻爾茲曼方程，介于微觀流體分子動力學(xué)模型和宏觀連續(xù)介質(zhì)模型之間的介觀方法，它既能捕捉到流體的微觀特性，又能再現(xiàn)流體的宏觀流動特性.在處理復(fù)雜邊界問題時，LBM具有極大的優(yōu)勢，因此它是一種可以用來模擬具有任意幾何邊界系統(tǒng)（例如多孔介質(zhì)）的非常有效的模擬方法.與有限差分方法等相比，LBM算法簡單、易于程序化處理、計算效率高和具有天然的并行性，可以在計算機(jī)或者集群上運行，減少了計算代價和計算時間.但是，LBM仍然是一個基于網(wǎng)格的方法，在多尺度問題中，如果計算區(qū)域復(fù)雜，尤其是在解決具有大畸變、運動物質(zhì)交界面、變形邊界以及自由表面問題時，容易產(chǎn)生網(wǎng)格扭曲與網(wǎng)格重構(gòu)問題.

一個可以替代LBM的方法是光滑粒子動力學(xué)（SPH）方法［1］.SPH 方法是一種無網(wǎng)格粒子方法，最初主要用于解決天文學(xué)問題［2－3］.經(jīng)過數(shù)十年的發(fā)展，它已經(jīng)廣泛應(yīng)用于不同領(lǐng)域內(nèi)的計算流體力學(xué)問題.Monaghan［4］對SPH方法應(yīng)用于流體力學(xué)的計算做了一系列的研究.Morris等［5］采用SPH方法研究了不可壓縮流動問題和多孔介質(zhì)中的流體流動問題，Zhu等［6－8］應(yīng)用 SPH 方法模擬了孔隙尺度下流體的流動和擴(kuò)散問題.

與基于網(wǎng)格的數(shù)值模擬方法相比，SPH方法主要具有以下兩方面的優(yōu)勢［9］：一是它的物理意義簡單明確.通過添加成對的粒子間相互作用，只需更改少量的代碼，即可將復(fù)雜的物理和化學(xué)變化包含在SPH模型中.因為總粒子數(shù)在模擬中保持恒定，同時相互作用是對稱的，質(zhì)量、動量和能量方程守恒性好，而且使用者可以簡單地在一個指定區(qū)域內(nèi)添加更多粒子從而得到更加精確的結(jié)果.二是可以相對簡單地處理復(fù)雜幾何邊界問題.該方法僅僅處理一些離散的點，在粒子的運動過程中，可以通過粒子和骨架結(jié)構(gòu)的相互作用反映固體邊界的影響，無需進(jìn)行額外的處理，因此它非常適用于研究復(fù)雜幾何邊界問題.

本文采用SPH方法，在孔隙尺度下模擬多孔介質(zhì)內(nèi)不可壓縮流體流動，研究流體流動的細(xì)節(jié)及其特性.

1 光滑流體動力學(xué)基礎(chǔ)

SPH方法是一種基于插值理論的無網(wǎng)格拉格朗日方法.在SPH方法中，系統(tǒng)狀態(tài)和運動由一系列隨機(jī)分布的離散粒子來表示.每一個粒子都具有質(zhì)量m、密度ρ、速度v，以及其他的流體性質(zhì).通過使用一系列任意分布的節(jié)點（或者粒子）來求解積分方程組或者偏微分方程組，從而得到精確穩(wěn)定的數(shù)值解.

在SPH方法中，函數(shù)f（r）的積分近似式為

式中：r為粒子坐標(biāo)；h為光滑長度；W（r－r′，h）為核函數(shù).

相應(yīng)的函數(shù)在某一插值點i的粒子近似式為

式中：求和操作應(yīng)用于核函數(shù)在插值點i支持域內(nèi)的所有N個粒子；Wij＝W（ri－rj，h），為核函數(shù).核函數(shù)選取應(yīng)用最廣泛且精度較高的五次光滑核函數(shù)：

SPH方法中的控制方程組是不可壓縮流體N－S方程組的特殊離散形式，其中，連續(xù)方程為

動量方程：

式中：p為壓力；μ為流體動力黏度；fb為單位質(zhì)量體積力；∑f為粒子與多孔介質(zhì)的相互作用力為一個很小的人工力，用于抑制張力不穩(wěn)定性［10］，其中

式中：1＝0.2，2＝0.05，n＝3；ΔS一般取初始粒子間距.

通過引用人工壓縮方法來處理不可壓縮流體，采用的狀態(tài)方程為

粒子按照以下形式運動：

為了實現(xiàn)無滑移邊界條件，采用一種反射對稱布置的虛粒子來模擬壁面.這些虛粒子具有與相對應(yīng)的實粒子相同的密度和壓力，但速度方向相反.這些虛粒子在每個計算步中由對應(yīng)的實粒子對稱產(chǎn)生.為阻止粒子穿越邊界，每個邊界上的虛粒子都添加了一個類似于計算分子力時所使用的Lennard－Jones排斥力，如下所示：

式中：n1＝12，n2＝4；D為由具體問題所定的參數(shù)，一般取與速度的最大值的平方相等的量級；r0為截止半徑.

本文采用周期性邊界條件模擬無限長通道中的流體流動問題.這意味著在粒子運動過程中，當(dāng)粒子即將離開指定區(qū)域穿越邊界時，有一個具有相同流體性質(zhì)（速度和密度）的新粒子在另一個邊界產(chǎn)生并進(jìn)入計算區(qū)域，進(jìn)而保證計算區(qū)域的粒子總數(shù)不變.

2 SPH方法的數(shù)值驗證

用經(jīng)典算例Poiseuille流和Couette流對計算方法進(jìn)行驗證，并將模擬結(jié)果和解析解進(jìn)行對比.本文中，計算域為0.001m×0.001m，流體為水，密度ρ＝1×103kg／m3，運動黏度ν＝1×10－6m2／s，為了減少計算誤差，使流體流動的特征速度接近mm／s的數(shù)量級，Poiseuille流中的水平方向單位質(zhì)量體積力fb＝0.02N／kg，Couette流中的頂邊運動速度v0＝2.5×10－3m／s，粒子數(shù)量為40×40，光滑長度取定值，為粒子初始間距的1.1倍，時間步長取為1×10－4s.

圖1 Poiseuille流速度分布Fig.1 Velocity profiles for the Poiseuille flow

圖2 Couette流速度分布Fig.2 Velocity profiles for the Couette flow

在經(jīng)歷6 000步計算后，模擬達(dá)到穩(wěn)定狀態(tài).圖1和2分別給出了Poiseuille流與Couette流，由SPH方法和級數(shù)求解法（解析法）得到的在t＝0.01、0.1s和最終狀態(tài)t＝∞時的速度分布圖，可見SPH方法所得到的結(jié)果和解析法所得到的結(jié)果相當(dāng)吻合，相對誤差在1%以內(nèi).

3 SPH方法在多孔介質(zhì)中流體流動的應(yīng)用

3.1 模擬規(guī)則多孔介質(zhì)中流體流動

計算域為0.001m×0.001m，初始粒子分布為64×64，x和y方向均采用周期性邊界條件，流體為水，運動黏度和密度分別為ν＝1×10－6m2／s和ρ＝1×103kg／m3，為了減少計算誤差，使流體流動的特征速度接近mm／s的數(shù)量級，水平方向單位質(zhì)量體積力取為fb＝0.2N／kg.本文中，分別采用方柱和圓柱兩種骨架結(jié)構(gòu)來構(gòu)造規(guī)則多孔介質(zhì)，方柱邊長和圓柱直徑分別為0.5mm和0.564mm，兩種多孔介質(zhì)的孔隙率均為0.75.

同時采用SPH方法和FEM對穩(wěn)態(tài)不可壓縮流體流動進(jìn)行模擬，得到的水平速度分布如圖3所示，曲線a和曲線b分別表示在x＝l／2和出口處的水平速度分布.兩種模擬方法所得到的結(jié)果吻合良好.

圖3 多孔介質(zhì)中SPH方法和FEM得到的水平速度分布比較Fig.3 Comparison of horizontal velocity obtained by SPH method and FEM in porous media

3.2 模擬貝雷巖心模型與人工隨機(jī)多孔介質(zhì)中的流體流動

計算域為0.001m×0.001m，初始粒子分布為64×64，x和y方向均采用周期性邊界條件，流體為水，運動黏度和密度分別為ν＝1×10－6m2／s和ρ＝1×103kg／m3.根據(jù)規(guī)則多孔介質(zhì)的模擬結(jié)果，水平方向單位質(zhì)量體積力分別取為fb＝0.2、0.4、0.6、0.8N／kg.在本文中，采用兩種不同的多孔介質(zhì)模型，第1種模型是基于貝雷巖心的核磁成像圖像數(shù)據(jù)構(gòu)建的，圖4為由核磁共振成像儀得到的自旋密度圖及將成像數(shù)據(jù)二值化處理后得到的圖像.考慮計算效率，只選取其中的局部區(qū)域作為研究對象，考察其孔隙結(jié)構(gòu)內(nèi)部的流體流動情況，所選取的巖心模型如圖5所示.

圖4 貝雷巖心自旋密度圖像與二值化圖像Fig.4 Spin density image and binary image of Berea sandstone

圖5 貝雷巖心模型（孔隙率為0.577）Fig.5 Berea sandstone model（porosity is 0.577）

第2種模型為采用隨機(jī)方法，以同尺度的圓柱為基元生成的多孔介質(zhì)模型，其孔隙率與貝雷巖心模型相同，如圖6所示.

多孔介質(zhì)的滲透率由下式確定：

圖6 人工隨機(jī)多孔介質(zhì)模型（孔隙率為0.577）Fig.6 Artificial random porous media model（porosity is 0.577）

圖7 表示穩(wěn)態(tài)流動時，流體平均速度和單位質(zhì)量體積力fb之間的關(guān)系.可以看出模擬結(jié)果呈現(xiàn)了很好的線性關(guān)系.根據(jù)圖7（a）、（b）中擬合直線的斜率和流體物性，由式（11）計算得到貝雷巖心模型和人工隨機(jī)多孔介質(zhì)的滲透率分別為k1＝2.74×10－11m2和k2＝4.64×10－10m2.這與通過達(dá)西定律所得到的預(yù)期結(jié)果相符.

圖7 平均流體速度和單位質(zhì)量體積力的關(guān)系Fig.7 Average flow velocity versus body force per unit mass

4 結(jié) 論

本文利用SPH方法模擬了經(jīng)典的Poiseuille流和Couette流，對所開發(fā)的SPH計算程序進(jìn)行了驗證，并分別利用SPH方法和FEM研究了流體在兩種不同的規(guī)則多孔介質(zhì)中的流動特性.在低速流動下，SPH方法計算結(jié)果和FEM的計算結(jié)果都非常一致，說明SPH模型可以在微觀尺度下直接模擬多孔介質(zhì)內(nèi)流體流動的詳細(xì)情況.最后，研究了流體在基于真實的貝雷巖心獲得的多孔介質(zhì)模型和以同尺度的圓柱為基元隨機(jī)構(gòu)造的多孔介質(zhì)模型內(nèi)的流體流動情況.結(jié)果表明，該SPH模型可以用于研究多孔介質(zhì)內(nèi)流體流動問題，能夠很好地模擬多孔介質(zhì)的微觀結(jié)構(gòu)對多孔介質(zhì)內(nèi)流體流動特性的影響，可以很容易地遷移應(yīng)用于任意形狀和大小的固體骨架結(jié)構(gòu)之中的流體流動問題.模擬得到的平均流體速度和單位質(zhì)量體積力有著很好的線性關(guān)系，符合通過達(dá)西定律所得到的預(yù)期結(jié)果.

［1］LIU Gui－rong，LIU Mou－bin.Smoothed Particle Hydrodynamics— A Meshfree Particle Method［M］.Singapore：World Scientific，2004.

［2］Gingold R A，Monaghan J J.Smoothed particle hydrodynamics — Theory and application to nonspherical stars ［J］.Monthly Notices of the Royal Astronomical Society，1977，181：375－389.

［3］Lucy L B.A numerical approach to the testing of the fission hypothesis ［J］.The Astronomical Journal，1977，82：1013－1024.

［4］Monaghan J J.Smoothed particle hydrodynamics［J］.Reports on Progress in Physics，2005，68（8）：1703.

［5］Morris J P，F(xiàn)ox P J，ZHU Yi.Modeling low Reynolds number incompressible flows using SPH［J］.Journal of Computational Physics，1997，136（1）：214－226.

［6］ZHU Yi， Fox P J.Smoothed particle hydrodynamics model for diffusion through porous media ［J］.Transport in Porous Media，2001，43（3）：441－471.

［7］ZHU Yi，F(xiàn)ox P J.Simulation of pore－scale dispersion in periodic porous media using smoothed particle hydrodynamics ［J ］.Journal of Computational Physics，2002，182（2）：622－645.

［8］ZHU Yi，F(xiàn)ox P J，Morris J P.A pore－scale numerical model for flow through porous media［J］.International Journal for Numerical and Analytical Methods in Geomechanics，1999，23（9）：881－904.

［9］胡曉燕.計算流體力學(xué)中的SPH方法和 MLSPH方法研究［D］.北京：中國工程物理研究院，2006.HU Xiao－Yan.Research on SPH and MLSPH method in computational fluid dynamics ［D］.Beijing：China Academy of Engineering Physics，2006.（in Chinese）

［10］Monaghan J J.SPH without a tensile instability［J］.Journal of Computational Physics，2000，159（2）：290－311.