多孔介質(zhì)孔隙尺度下不可壓縮流體流動(dòng)特性SPH模擬

李 維 仲， 趙 月 帥， 宋 永 臣

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

動(dòng)量方程：

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

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

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

粒子按照以下形式運(yùn)動(dòng)：

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

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

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

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

用經(jīng)典算例Poiseuille流和Couette流對(duì)計(jì)算方法進(jìn)行驗(yàn)證，并將模擬結(jié)果和解析解進(jìn)行對(duì)比.本文中，計(jì)算域?yàn)?.001m×0.001m，流體為水，密度ρ＝1×103kg／m3，運(yùn)動(dòng)黏度ν＝1×10－6m2／s，為了減少計(jì)算誤差，使流體流動(dòng)的特征速度接近mm／s的數(shù)量級(jí)，Poiseuille流中的水平方向單位質(zhì)量體積力fb＝0.02N／kg，Couette流中的頂邊運(yùn)動(dòng)速度v0＝2.5×10－3m／s，粒子數(shù)量為40×40，光滑長(zhǎng)度取定值，為粒子初始間距的1.1倍，時(shí)間步長(zhǎng)取為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步計(jì)算后，模擬達(dá)到穩(wěn)定狀態(tài).圖1和2分別給出了Poiseuille流與Couette流，由SPH方法和級(jí)數(shù)求解法（解析法）得到的在t＝0.01、0.1s和最終狀態(tài)t＝∞時(shí)的速度分布圖，可見SPH方法所得到的結(jié)果和解析法所得到的結(jié)果相當(dāng)吻合，相對(duì)誤差在1%以內(nèi).

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

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

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

同時(shí)采用SPH方法和FEM對(duì)穩(wěn)態(tài)不可壓縮流體流動(dòng)進(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ì)中的流體流動(dòng)

計(jì)算域?yàn)?.001m×0.001m，初始粒子分布為64×64，x和y方向均采用周期性邊界條件，流體為水，運(yùn)動(dòng)黏度和密度分別為ν＝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ù)二值化處理后得到的圖像.考慮計(jì)算效率，只選取其中的局部區(qū)域作為研究對(duì)象，考察其孔隙結(jié)構(gòu)內(nèi)部的流體流動(dòng)情況，所選取的巖心模型如圖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)流動(dòng)時(shí)，流體平均速度和單位質(zhì)量體積力fb之間的關(guān)系.可以看出模擬結(jié)果呈現(xiàn)了很好的線性關(guān)系.根據(jù)圖7（a）、（b）中擬合直線的斜率和流體物性，由式（11）計(jì)算得到貝雷巖心模型和人工隨機(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流，對(duì)所開發(fā)的SPH計(jì)算程序進(jìn)行了驗(yàn)證，并分別利用SPH方法和FEM研究了流體在兩種不同的規(guī)則多孔介質(zhì)中的流動(dòng)特性.在低速流動(dòng)下，SPH方法計(jì)算結(jié)果和FEM的計(jì)算結(jié)果都非常一致，說明SPH模型可以在微觀尺度下直接模擬多孔介質(zhì)內(nèi)流體流動(dòng)的詳細(xì)情況.最后，研究了流體在基于真實(shí)的貝雷巖心獲得的多孔介質(zhì)模型和以同尺度的圓柱為基元隨機(jī)構(gòu)造的多孔介質(zhì)模型內(nèi)的流體流動(dòng)情況.結(jié)果表明，該SPH模型可以用于研究多孔介質(zhì)內(nèi)流體流動(dòng)問題，能夠很好地模擬多孔介質(zhì)的微觀結(jié)構(gòu)對(duì)多孔介質(zhì)內(nèi)流體流動(dòng)特性的影響，可以很容易地遷移應(yīng)用于任意形狀和大小的固體骨架結(jié)構(gòu)之中的流體流動(dòng)問題.模擬得到的平均流體速度和單位質(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］胡曉燕.計(jì)算流體力學(xué)中的SPH方法和 MLSPH方法研究［D］.北京：中國(guó)工程物理研究院，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.