黃志江 呂孝飛 趙會(huì)軍 宋琳琳 楊草來 張建龍

(常州大學(xué)石油工程學(xué)院，江蘇省油氣儲(chǔ)運(yùn)技術(shù)重點(diǎn)實(shí)驗(yàn)室 江蘇常州 213164)

0 引言

油品在土壤、砂礫等多孔介質(zhì)中的流動(dòng)阻力特性研究在實(shí)際工程應(yīng)用中具有重要意義，如研究埋地油品管道泄漏后的遷移擴(kuò)散規(guī)律、清除泄漏后的油污區(qū)域等問題時(shí)，油品流經(jīng)土壤、砂礫時(shí)的阻力壓降是進(jìn)行相關(guān)研究的重要參數(shù)[1-2]。目前，國(guó)內(nèi)外學(xué)者借助實(shí)驗(yàn)和數(shù)值模擬方法，對(duì)多孔介質(zhì)中的流體流動(dòng)阻力特性進(jìn)行了研究并有了一定的認(rèn)識(shí)。PASCAL J P等[3]研究表明，當(dāng)流體以低雷諾數(shù)流動(dòng)時(shí)，壓力梯度主要用于克服黏性阻力，當(dāng)流速增大到一定限度，慣性阻力的影響增大，壓力梯度需同時(shí)克服黏性阻力和慣性阻力；JAN H[4]研究了單相水在各種均勻分級(jí)顆粒多孔介質(zhì)中的流動(dòng)行為，發(fā)現(xiàn)多孔介質(zhì)中孔隙結(jié)構(gòu)的幾何形狀對(duì)流動(dòng)阻力影響很大；SYLWIA R[5]發(fā)現(xiàn)聚合物溶液通過多孔介質(zhì)時(shí)，會(huì)受到剪切和拉伸致使流動(dòng)阻力大幅增加；張震等[6]對(duì)多孔介質(zhì)內(nèi)單相流的絕熱流動(dòng)阻力進(jìn)行了數(shù)值模擬，分析了慣性損失系數(shù)對(duì)流動(dòng)阻力的影響；劉學(xué)強(qiáng)等[7]對(duì)紊流流態(tài)區(qū)間下多孔介質(zhì)中單相流的阻力特性進(jìn)行了實(shí)驗(yàn)研究，得出相應(yīng)的數(shù)學(xué)模型；李振鵬等[8]采用玻璃球填充的多孔介質(zhì)通道對(duì)單相水的流動(dòng)阻力特性進(jìn)行實(shí)驗(yàn)研究，得出相應(yīng)阻力壓降計(jì)算公式；于立章等[9-10]應(yīng)用相似理論，合理地縮減多孔介質(zhì)通道的尺寸，建立了多孔介質(zhì)通道中單相水絕熱流動(dòng)壓降預(yù)測(cè)模型；吳國(guó)忠、齊晗兵等[11-12]測(cè)量了油水混合液在顆粒多孔介質(zhì)通道內(nèi)的流動(dòng)阻力，擬合了其壓降-速度表達(dá)式，并得到其黏性和慣性阻力系數(shù)；張淑淑等[13]研究了儲(chǔ)罐內(nèi)填充不同的多孔材料厚度下，液化天然氣在儲(chǔ)罐內(nèi)的流動(dòng)情況。以上這些研究結(jié)果使人們對(duì)多孔介質(zhì)中單相流的流動(dòng)規(guī)律有了比較深入的認(rèn)識(shí)，但對(duì)兩相流及多相流在顆粒多孔介質(zhì)中的流動(dòng)特性的認(rèn)識(shí)較少。本文考慮到油水兩相流比單相流的流動(dòng)形態(tài)更為復(fù)雜，油水相比可能對(duì)阻力壓降的影響顯著，多孔介質(zhì)中孔隙分布具有隨機(jī)性，流經(jīng)不同長(zhǎng)度多孔介質(zhì)通道時(shí)產(chǎn)生的阻力壓降也應(yīng)有較大差異，故對(duì)油水兩相的相比和多孔介質(zhì)通道長(zhǎng)度進(jìn)行對(duì)比分析，得到它們對(duì)流動(dòng)阻力壓降的影響程度，對(duì)實(shí)際工程應(yīng)用有很好的指導(dǎo)意義。

本文構(gòu)建小球顆粒填充多孔介質(zhì)通道幾何模型并結(jié)合CFD數(shù)值模擬計(jì)算，探討不同油水比工況和不同多孔介質(zhì)通道長(zhǎng)度工況下，流動(dòng)阻力壓降隨入口流速的變化規(guī)律，為輸油管道泄漏污染物在多孔介質(zhì)中的流動(dòng)特性提供可靠的依據(jù)。

1 數(shù)值模擬

1.1 幾何模型及網(wǎng)格劃分

為保證多孔介質(zhì)通道幾何模型的可行性，采用參考文獻(xiàn)[6]中的實(shí)驗(yàn)參數(shù)，即內(nèi)徑為0.05 m，管長(zhǎng)為1 m的有機(jī)玻璃管內(nèi)緊密填充直徑8 mm的均勻粒徑玻璃小球顆粒，其孔隙率為0.408。模型尺寸以實(shí)驗(yàn)所用多孔介質(zhì)通道尺寸為依據(jù)，考慮到網(wǎng)格數(shù)量及網(wǎng)格質(zhì)量的影響后，按照2.5∶1的比例繪制上述多孔介質(zhì)通道簡(jiǎn)化的幾何模型。為簡(jiǎn)化網(wǎng)格數(shù)量和減少計(jì)算量，根據(jù)圓管模型的對(duì)稱性，可選取1/4圓管通道作為模擬流域，小球顆粒以四棱錐結(jié)構(gòu)為基本單元進(jìn)行填充，在近壁面間隙較大部分，使用不同體積的半球進(jìn)行填充，如圖1所示。

(a)俯視圖(局部) (b)左視圖(局部)

(c)主視圖 (d)斜視圖圖1 多孔介質(zhì)通道小球排列示意

在管道無障礙空間內(nèi)作穩(wěn)定流動(dòng)的流體流入多孔介質(zhì)通道時(shí)，會(huì)經(jīng)過一個(gè)過渡段(過渡段長(zhǎng)度l與管道內(nèi)徑d相等)才能成為速度均勻發(fā)展的穩(wěn)定流動(dòng)，如圖2所示[14]。可以適當(dāng)縮小幾何模型中多孔介質(zhì)通道穩(wěn)定流動(dòng)段長(zhǎng)度來合理簡(jiǎn)化模型。根據(jù)以上分析，確定模型入口空管段長(zhǎng)度為0.02 m，小球填充多孔介質(zhì)段長(zhǎng)度為0.045 m，出口空管段長(zhǎng)度為0.02 m，如圖3所示。

圖2 流體流入多孔介質(zhì)示意

圖3 流動(dòng)通道示意

為了保證網(wǎng)格正常生成，在構(gòu)建幾何模型時(shí)，需使小球顆粒與小球顆粒之間留有微小的縫隙。因此，本文對(duì)徑向的小球分布采取相互分離的排列方式，對(duì)軸向的每層小球采取鑲嵌的排列方式。通過這種正負(fù)誤差抵消的方法，可以降低整個(gè)通道內(nèi)的誤差[9]。

由于模型不規(guī)則程度高，不易生成計(jì)算量大的六面體網(wǎng)格，故選用非結(jié)構(gòu)四面體網(wǎng)格進(jìn)行劃分，為保證計(jì)算結(jié)果的精確，對(duì)小球顆粒間的縫隙區(qū)域生成的網(wǎng)格進(jìn)行加密處理，如圖4所示。經(jīng)網(wǎng)格無關(guān)驗(yàn)證，確定網(wǎng)格數(shù)為4 326 709個(gè)，其中95%的網(wǎng)格的偏斜度為0到0.25，網(wǎng)格質(zhì)量滿足計(jì)算要求。

圖4 顆粒間網(wǎng)格示意

1.2 數(shù)值模型

1.2.1 歐拉模型

歐拉方法對(duì)每一相求解動(dòng)量方程和連續(xù)性方程，并通過相間作用力來實(shí)現(xiàn)相間耦合[15]，能夠有效地模擬相間混合，對(duì)于不可壓縮流體來說，每一相的平均質(zhì)量和動(dòng)量方程表示如下：

(1)連續(xù)性方程

(1)

式中，α、ρ和u分別表示體積分?jǐn)?shù)、密度和速度；下標(biāo)c表示連續(xù)相的油或水。

(2)體積分?jǐn)?shù)守恒方程

多孔介質(zhì)通道中的油水兩相流，只存在油相和水相，因此各自體積分?jǐn)?shù)αs，αl任意時(shí)刻都服從：

αs+αl=1

(2)

(3)動(dòng)量方程

(3)

式中，G為重力；FM是兩相間的轉(zhuǎn)換動(dòng)量或兩相之間的作用力，包括曳力、升力、虛擬質(zhì)量力和湍流耗散力；τ為應(yīng)力張量；上標(biāo)l和t分別表示層流和湍流。

1.2.2 湍流模型

Realizablek-ε湍流模型將流場(chǎng)的旋轉(zhuǎn)以及彎曲壁面流動(dòng)考慮在內(nèi)，可用于中等強(qiáng)度旋流的流場(chǎng)。故本文選用Realizablek-ε湍流模型，其表達(dá)式為

(4)

(5)

式中，t表示時(shí)間，ρ表示密度，Gk與Gb均表示湍動(dòng)能，YM表示湍流的脈動(dòng)分量在運(yùn)動(dòng)過程中對(duì)耗算率的影響，ut表示湍流粘性系數(shù)。Clε與C2ε是常量，σk為湍動(dòng)能的Prandtl數(shù)，σε是耗算率的Prandtl數(shù)。Clε=1.44，C2ε=1.29，σk=1.0，σε=1.2。

1.2.3 邊界條件及求解器設(shè)置

邊界條件設(shè)置如下：入口條件，設(shè)為速度邊界條件，分別設(shè)置5%、30%、60%、90% 4組油水體積分?jǐn)?shù)。當(dāng)油水體積分?jǐn)?shù)大于50%時(shí)，主相為油相，密度為730 kg/m3，動(dòng)力粘度為2.4×10-3Pa·s；當(dāng)油水體積分?jǐn)?shù)小于50%時(shí)，主相為水相，密度為998.2 kg/m3，動(dòng)力粘度為1×10-3Pa·s；出口條件，設(shè)為自由出流；壓力條件，僅對(duì)油水兩相進(jìn)行數(shù)值模擬研究，不考慮氣相，因此操作壓力設(shè)為101.325 kPa，壓力參考點(diǎn)設(shè)置在出口邊界面，防止結(jié)果出現(xiàn)負(fù)壓的情況；圓管對(duì)稱邊界面設(shè)置為symmetry邊界條件；壁面無滑移，近壁面采用標(biāo)準(zhǔn)壁面函數(shù)處理。

采用pressure-based求解器以及適于壓力和速度耦合的SIMPLE算法進(jìn)行求解，各控制方程采用二階迎風(fēng)離散格式；為使計(jì)算結(jié)果收斂更快可適當(dāng)調(diào)低松弛因子；監(jiān)測(cè)殘差設(shè)置為10-5。

2 結(jié)果與分析

2.1 模擬結(jié)果

圖5為4種不同油水比工況下，等差值多孔介質(zhì)長(zhǎng)度下阻力壓降的模擬值?？梢钥闯觯S著入口流速的增大，不同油水比之間的阻力壓降差也變得更加顯著，每段多孔介質(zhì)中阻力壓降隨入口流速的增大呈指數(shù)式上升趨勢(shì)，在入口流速較低時(shí)，即層流流態(tài)，油水混合液主要受黏性阻力，流速-阻力壓降曲線呈一次函數(shù)曲線，當(dāng)流速越大時(shí)，雷諾數(shù)超過一定臨界值時(shí)，慣性阻力項(xiàng)主導(dǎo)流動(dòng)阻力，流速-阻力壓降曲線呈二次函數(shù)曲線。

(a)5%油水比

(b)30%油水比

(c)60%油水比

(d)90%油水比圖5 不同油水比下每段多孔介質(zhì)的壓降

2.2 油水比對(duì)多孔介質(zhì)阻力壓降的影響

圖6是等體積分?jǐn)?shù)差的油水混合液下，阻力壓降隨入口流速變化的模擬值。由圖可以看出，入口流速較小時(shí)，不同油水比下的壓降變化均較小，壓降與流速均成一次函數(shù)關(guān)系，但隨著入口流速的增大，阻力壓降變化幅度也逐漸增大。按等體積分?jǐn)?shù)改變油水比，入口流速為0.05 m/s時(shí)，阻力壓降改變較小，5%油水比的最小壓降與90%油水比的最小壓降僅相差4 kPa，而隨著入口流速增大，在0.5 m/s時(shí)，阻力壓降產(chǎn)生巨大改變，5%油水比的最大壓降與90%油水比的最大壓降相差接近147 kPa。同時(shí)可以看出，油水比一定時(shí)，隨入口流速的增大，流速-阻力壓降曲線呈現(xiàn)二次函數(shù)曲線。

圖6 不同油水比下入口流速與壓降的關(guān)系

研究表明，相鄰體積分?jǐn)?shù)差值的油水比下，相同入口流速對(duì)應(yīng)的阻力壓降之間差距大致是相等的，說明在不同流速時(shí)，油水比的變化不會(huì)改變阻力壓降上升的變化趨勢(shì)，但油水比越小時(shí)，阻力壓降值隨流速變化得越快。

2.3 多孔介質(zhì)通道長(zhǎng)度對(duì)阻力壓降的影響

圖7是以5%油水比工況為例，按等差值截取的多孔介質(zhì)通道長(zhǎng)度下，阻力壓降隨入口流速變化的模擬值。由圖可以看出，入口流速較小時(shí)，不同多孔介質(zhì)長(zhǎng)度下的壓降變化均較小，流速-阻力壓降曲線均呈一次函數(shù)曲線關(guān)系。多孔介質(zhì)長(zhǎng)度較長(zhǎng)時(shí)，隨著入口流速的增大，阻力壓降變化也逐漸增大。按等差值截取多孔介質(zhì)長(zhǎng)度，入口流速為0.05 m/s時(shí)，阻力壓降改變較小，9 mm長(zhǎng)度的最小壓降與45 mm長(zhǎng)度的最小壓降僅相差3.4 kPa，而隨著入口流速增大，在0.5 m/s時(shí)，阻力壓降產(chǎn)生更大改變，9 mm長(zhǎng)度的最大壓降與45 mm長(zhǎng)度的最大壓降相差接近192 kPa。同時(shí)可以看出，多孔介質(zhì)長(zhǎng)度一定時(shí)，隨入口流速的增大，阻力壓降-流速曲線的斜率越大。這是由于流動(dòng)阻力項(xiàng)中的黏性阻力項(xiàng)和慣性阻力項(xiàng)會(huì)隨層流區(qū)和紊流區(qū)的流態(tài)變化，依次成為主導(dǎo)項(xiàng)。同時(shí)，從圖7可以看出，相鄰差值的多孔介質(zhì)長(zhǎng)度下，相同入口流速對(duì)應(yīng)的阻力壓降之間差距大致是相等的，說明在不同流速時(shí)，多孔介質(zhì)長(zhǎng)度的變化不能改變阻力壓降上升的變化趨勢(shì)。

圖7 不同多孔介質(zhì)長(zhǎng)度下入口流速與壓降的關(guān)系

3 結(jié)論

通過油水兩相在顆粒多孔介質(zhì)中的流動(dòng)阻力特性三維數(shù)值模擬研究，得到不同油水比和等差值多孔介質(zhì)長(zhǎng)度工況下，阻力壓降的變化規(guī)律。

(1)不同入口流速下，油水比越大，多孔介質(zhì)進(jìn)出口阻力壓降反而越小；當(dāng)入口流速較低時(shí)，不同油水比下的進(jìn)出口阻力壓降相差較小，流速-阻力壓降曲線呈一次函數(shù)曲線；隨入口流速逐漸增大，不同油水比下的進(jìn)出口阻力壓降差值越顯著，流速-阻力壓降曲線呈二次函數(shù)曲線。

(2)不同入口流速下，多孔介質(zhì)長(zhǎng)度越長(zhǎng)，多孔介質(zhì)進(jìn)出口阻力壓降越大；當(dāng)入口流速較低時(shí)，不同多孔介質(zhì)長(zhǎng)度下的進(jìn)出口阻力壓降相差較小，流速-阻力壓降曲線呈一次函數(shù)曲線；隨入口流速逐漸增大，不同多孔介質(zhì)長(zhǎng)度下的進(jìn)出口阻力壓降相差越大，流速-阻力壓降曲線呈二次函數(shù)曲線。

(3)按等體積分?jǐn)?shù)差改變油水比以及等差值截取多孔介質(zhì)長(zhǎng)度，在相同入口流速時(shí)，相鄰工況的阻力壓降差值大致相等。這說明入口流速一定時(shí)，油水比以及多孔介質(zhì)長(zhǎng)度的等差值改變并不能改變阻力壓降的上升趨勢(shì)。