基于多孔躍遷模型的流體阻力壓降特性研究

劉洋，趙立新,周龍大，季豪

(1.東北石油大學機械科學與工程學院，黑龍江大慶 163318；2.黑龍江省石油石化多相介質處理及污染防治重點實驗室，黑龍江大慶 163318)

0 前言

多孔介質作為流體滲流運移的載體，是自然界中廣泛存在的一種具有相互聯系的孔隙介質材料。從小型芯片材料到大型望遠鏡面板等，多孔介質的相關理論模型被廣泛應用于各個領域。BASTIAN將多孔介質定義為由持續(xù)存在的固體部分和相互連通的孔隙組成的區(qū)域?？紫吨锌梢蕴畛湟环N或多種流體，如油、水和氣體等；固相可以是石英砂、玻璃珠、硅藻土等非固結的填充顆粒集合體，也可以是相互連接的柔性或剛性骨架結構。

國內外研究者利用實驗測試和數值模擬方法，對非固結型填充床多孔介質中的流體流動阻力特性開展深入分析討論并取得一定的研究成果，而對固結型多孔介質區(qū)域內的阻力壓降和阻力系數的相關性研究較少。由于多孔介質復雜的內部結構，很難通過解析方法推導出阻力公式，因此，將多孔介質假設為一種虛擬的、均勻的連續(xù)介質，不同流速下的流體分子間以相互碰撞的形式實現動量交換，簡化流體在多孔介質內的流動問題。多孔介質中流體流動的壓降主要由黏性力和慣性力兩部分組成；低雷諾數時的慣性力影響可以忽略不計，而高雷諾數時的慣性力影響顯著。慣性力扭曲了流線，從而增加了速度梯度，導致壓降增加。當速度增加時，流動進入非線性層流狀態(tài)，此時多孔介質慣性效應不再可以忽略。

本文作者基于對多孔介質區(qū)域阻力壓降方程的分析，研究固結型金屬濾網內流體的阻力壓降特性。設計一種用于測量求解多孔介質區(qū)域內流體阻力系數的實驗裝置，得到不同孔徑條件下流體通過銅網的阻力壓降和流速的實驗值?；趯嶒灁祿\用二項式插值非線性擬合得到的阻力系數，采用多孔躍遷模型數值模擬多孔介質區(qū)域的阻力壓降特征，得到特定流速條件下的阻力壓降模擬值；比較分析阻力壓降的實驗值和模擬值，驗證了多孔介質阻力系數測定裝置和方法的準確性和便捷性；對銅網區(qū)域內流場的速度場、壓力場特性進行分析，驗證多孔躍遷模型的適用性和準確性，為研究多孔介質區(qū)域內流體流動的阻力壓降特性提供了實驗測量和數值模擬的新途徑。

1 阻力壓降模型方程

1.1 阻力源項組成

計算流體力學的控制方程通常由非穩(wěn)態(tài)項、對流項、擴散項、源項組成，如式(1)所示

(1)

其中，連續(xù)性方程、動量方程和能量方程的各系數變量如表1所示。其中：非穩(wěn)態(tài)項是由于時間的改變引起的變化；對流項是由于速度的改變引起的變化，即空間的一階導數引起的；擴散項則是由空間的二階導數引起的變化；而源項是發(fā)生在體積上的變化。

表1 控制方程的系數變量組成

多孔介質模型是計算流體動力學中一種常用的阻力介質模型，通常更多關注的是多孔介質區(qū)域對流體流動的影響。相對于標準流體流動方程，對于簡單多孔介質模型的建模，通過增加動量源項，實現將流體在區(qū)域中受到的阻力轉化為一種附加的分散阻力，將多孔介質區(qū)域簡化為增加了阻力源項的流體區(qū)域，利用多孔介質模型模擬流體流動在多孔介質區(qū)域內對動量方程的影響，如式(2)所示，動量源項作用于流體產生壓力梯度。

(2)

當在直角坐標系下時，動量源項表達式為

(3)

式中：Δ為多孔介質區(qū)域長度;為、、三個矢量方向的動量源項；和為系數矩陣。

當在柱坐標系下時，動量源項可以表示為

(4)

1.2 阻力壓降方程

多孔介質區(qū)域內流體流動過程中，跨層壓力梯度是系統幾何形狀、床層孔隙度、床層滲透率和流體物性的函數。DARCY(1956)基于在砂床上的體積流量和壓力差所進行的實驗測試，提出了一個估算體積流量的經驗方程，描述流體在多孔介質中的流量、黏度和壓力變化之間的比例關系。各向同性介質的達西定律表示為

(5)

對于各向異性的介質，達西定律可以表示為

(6)

在達西模型中，流體壓降與多孔介質中流體速度呈線性相關關系。已有多項研究證實并完善了多孔介質中流動的達西定律，并通過實驗發(fā)現，阻力壓降與流體速度的二次項有關。通常用Forchheimer方程來表示為

(7)

其中：Δ表示壓降，Pa；表示流動方向介質長度，m；表示流體密度，kg/m；表示流體動力黏度系數，kg/(m·s)；表示流體速度，m/s。通常將和分別稱為達西滲透系數和非達西滲透系數。為了描述多孔介質中的非線性流動，DUPUIT(1863)和FORCHHEIMER(1901)引入了具有、系數的二次項方程來描述流動方程。

(8)

Forchheimer方程是通過實驗獲得的經驗方程，僅給出了非線性方程的基本形式，而沒有討論方程中參數和的物理意義和影響因素，因此參數和僅能通過流體阻力的實驗數據獲得。

ERGUN(1952)結合Blake-Kozeny層流方程和Burke-Plummer湍流方程提出了如式(9)所示的Ergun方程，從流體速度、流體性質、孔隙率、流動方向、顆粒尺寸、顆粒形狀、顆粒狀固體表面等方面對該方程進行了檢驗。方程所適用的雷諾數范圍為0.4～10。

(9)

Ergun方程有效地解釋了同時發(fā)生的慣性和黏性能量損失，已獲得了各種顆粒形狀的大量壓降-流速數據的相關性。然而，有學者認為Ergun系數=150和=1.75不是常數，而是取決于雷諾數、孔隙度或顆粒形狀。因此，在研究多孔介質的流動特性時，需要對Ergun關系系數進行修正，以使其適用于不同系統，通常采用擬合阻力壓降和速度的實驗值方法獲得Ergun方程系數。

1.3 過濾阻力模型

如圖1所示，根據直角坐標系下達西定律和流體體積流量和過濾截面的關系式，建立直角坐標系下水平管段內流體過濾過程的數學模型：

圖1 直角坐標系下過濾介質內流體流動阻力模型

(10)

當過濾截面為圓形或橢圓形時，利用微元法選取長半軸和短半軸的微元半徑分別為d和d，則過濾截面面積表達式為

(11)

當長軸和短軸半徑分別為和時，將和分別關于取值范圍[0,]和[0,]積分，得到過濾截面面積表達式：

(12)

在處取厚度為d的薄層流體，過濾介質厚度(=0)=，(=)=0，過濾介質的阻力可表示為

(13)

將過濾截面積和過濾阻力代入達西公式，可以得到當流體通過的過濾介質厚度為時，流體在多孔介質區(qū)域內過濾產生的過濾阻力表示為

(14)

因此，流體在水平管段多孔介質區(qū)域內流動時，流動阻力與流體的體積流量成反比，與流體在多孔介質區(qū)域內產生的阻力壓降Δ成正比，此外，流體動力黏度系數越大，流動阻力越小，流動區(qū)域內截面尺寸越大，流體所受到的流動阻力越大。

1.4 阻力系數求解

在多孔介質模擬區(qū)域復雜性以及流體物理特性的綜合影響下，根據Buckingham提出的π定理法，假定阻力系數與液體動力黏度、液體密度、流量、孔隙率、孔隙大小、多孔介質區(qū)域水力直徑、多孔介質區(qū)域單位長度的壓降等有關。流體流動各參數間存在的關系為

(,,,,,,,Δ)=0

(15)

從上述8個物理量中選取3個基本變量：密度，流速，孔隙大小，根據π定理和量綱和諧原理，得到量綱為一關系方程：

(16)

(17)

(18)

=()()

(19)

其中:、為無量綱量;、、、為待定系數。因此，多孔介質的黏性阻力系數和慣性阻力系數模型與雷諾數和歐拉數相關，可根據實驗數據回歸得到量綱為一的量和待定系數。

準確的阻力系數對多孔介質模擬的設置至關重要。對于Ergun方程而言，通常采用公式法和二項式擬合法計算流體在多孔介質區(qū)域內的阻力系數。其中，公式法具體為利用Ergun方程計算得到阻力系數的解析解，如式(9)和式(20)；在沒有相關實驗數據的情況下,可以根據這些公式計算多孔介質的阻力系數，如式(21)和式(22)所示，分析流體流動的阻力特性。

(20)

(21)

(22)

二項式擬合法具體為利用流體在多孔介質區(qū)域內的多組流速和壓降實驗值，根據流體流動截面面積將流量轉換為表觀速度，根據實驗測量的多孔介質在流動方向上的長度將壓降轉換為壓力梯度，進行二項式插值擬合，得到擬合常數和，再計算出阻力系數。

Δ=+

(23)

如果在數值模擬前已開展實驗研究，利用對壓力降與流體流速的實驗數據進行二項式插值擬合的方法可確定以黏性阻力系數和慣性阻力系數為待求解參數，以流體在試驗管段內的流速為自變量，流體在試驗管段兩端所產生的壓力降值Δ為函數對阻力系數進行多項式擬合。

(24)

(25)

2 阻力壓降測量系統實驗

2.1 實驗原理和方法

阻力壓降實驗采用80目、120目、150目和200目的黃銅網作為固結型多孔介質區(qū)域，多孔介質的孔隙率被定義為孔隙體積與總體積的比值。圖2所示為黃銅網結構單元尺寸，濾網的目數與濾網絲徑和孔徑有關，根據假設所選黃銅網為絲徑和孔徑為均勻濾網的結構特點，濾網的孔隙率可由式(28)計算得到，各尺寸銅網孔隙率如表2所示。

圖2 濾網結構單元尺寸圖

表2 銅網尺寸參數和孔隙率

單位濾網的表面積和體積分別為式(26)和(27)所示

=π

(26)

(27)

其中:為單根金屬絲的直徑，又稱濾網厚度，m；為濾網單元的邊長，又稱濾網孔徑，m。濾網的目數可表示為在25.4 mm長度上排列著個孔隙，則該濾網為目，由此可知，濾網孔隙率為

(28)

其中:為整個濾網的孔隙體積，m；為整個濾網的體積，m。

根據固結型和非固結型多孔介質材料的結構特點，設計一種用于流體通過不同類型多孔介質時，測量計算阻力系數的實驗裝置，如圖3所示。該裝置特點在于可以測量固結型的金屬泡沫、濾網等，也可以測量非固結型填充床的阻力系數。

圖3 多孔介質阻力系數測量裝置

該裝置由儲液系統、動力系統、監(jiān)測系統、數據采集處理系統、測量管段系統等組成，選取金屬銅濾網作為固結型多孔介質材料進行實驗。水罐內安裝有加熱管和溫控熱電偶元件，實現對罐內液體的控溫調節(jié)。測量管段系統由材質、尺寸完全相同的測壓管段1和測壓管段2組成，管段2內通過填充顆粒和法蘭固定濾網的方式構建不同結構多孔介質區(qū)域。如圖4所示，兩組測量管段分別安裝連接壓差變送器的高壓隔膜端和低壓隔膜端，實現混合液通過管段1時產生的管段阻力壓降損失和混合液通過管段2內多孔介質區(qū)域時產生的阻力壓降損失的測量。

圖4 壓差變送器測壓端隔膜法蘭式連接裝置

根據流量計測得管段1、2內通過的體積流量和，可得到相應流速分別為和；根據壓降-流速的相關性對測量管段1所得的壓降損失和流速采用插值法進行二項式擬合，得到式(23)中擬合系數和;再將由測量管段2所得的流速代入擬合關系式，可得到管段1與管段2在相同流速條件下所產生的阻力壓降損失，因此，流體通過管段2內多孔介質區(qū)域時產生的阻力壓降損失為

Δ=-

(29)

將測壓管段2的壓降和流速實驗值代入式(20)中，擬合得到管段2內的相關阻力系數。

2.2 實驗結果及分析

20 ℃水分別通過80目、120目、150目、200目銅網時，得到的流速和壓降的實驗值如圖5所示。測量銅網壓降均隨入口速度的增大而增大，入口速度較小時，阻力壓降變化幅度較小，隨入口速度的增大，阻力壓降變化幅度逐漸增大。由孔隙率測量計算可知，隨著銅網從80目增大至200目，單位長度上的孔隙增多，單位孔徑尺寸減小，孔隙率減?。幌嗤俣葪l件下，濾網壓降隨著目數的增加而增大。利用不同孔徑銅網的壓降和速度實驗值得到的阻力系數擬合結果列于表3中。隨著銅網目數的增大，黏性阻力系數和慣性阻力系數均逐漸增大，由相關指數分別為99.823%、99.533%、99.356%、99.674%，說明壓降變化的差異中大于99%都是由速度變化引起的，方程擬合較好。

圖5 20 ℃水通過不同規(guī)格濾網阻力壓降變化

表3 黏性阻力系數和慣性阻力系數二項式插值非線性擬合結果

當單位長度銅網的孔隙數增大時，孔徑尺寸減小，孔隙率降低。流體通過濾網所受到的阻力增大，流體的滲透性降低。根據達西定律可知，流體的黏性阻力系數與滲透率成反比，因此，隨著濾網孔隙數的增大，黏性阻力系數增大。由于孔隙變化使得流體流動截面縮小，在濾網內產生速度梯度，加速運動的流體引起附加阻力，同時引起流體動能變化，附加阻力與加速度成正比，阻力系數即為慣性阻力系數，隨著銅網目數的增大，慣性阻力系數增大。

3 多孔躍遷模型數值模擬

通過對多孔介質模型中黏性阻力系數和慣性阻力系數的設置，模擬流體通過金屬濾網為多孔介質材料時的阻力壓降變化情況。阻力系數的準確性和多孔躍遷模型的有效性，都對數值模擬的計算結果產生較大的影響，通過設定實驗測量數據擬合得到阻力系數，模擬計算阻力壓降值，比較實驗值和模擬值，驗證多孔躍遷模型的可行性。

3.1 幾何模型和網格劃分

利用SolidWorks建模軟件建立圓管幾何模型，其中直徑為25 mm，長度為700 mm；利用ICEM CFD劃分結構網格，如圖6所示。

圖6 幾何模型和網格劃分

3.2 模擬條件和工況

多孔躍遷模型被用來模擬一個具有已知速度或壓降特性的面區(qū)域，而不是單元區(qū)域，本質上是對結構區(qū)域的多孔介質模型的一維簡化。多孔躍遷模型與完整的多孔介質模型相比更穩(wěn)健，能產生更好的收斂效果。

(30)

其中:是層流流體的黏度；是介質表面的滲透率；是壓力-躍遷系數；是多孔面的法線速度；Δ是介質的厚度。當多孔躍遷是將兩個流體區(qū)分開，或在一個流體區(qū)形成一個內部區(qū)域時，多孔躍遷定義的厚度和系數參數被用來計算跨越躍遷的壓力降，這個壓降被轉換為一個不包含黏性力的力矢量。

多孔躍遷模型需要設置多孔介質面區(qū)域的表面滲透率、多孔介質的有限厚度、壓力躍遷系數，其中，表面滲透率與黏性阻力系數成倒數關系，壓力躍遷系數等于慣性阻力系數，介質厚度為絲徑大小，具體流體及工況參數詳見表4。求解器采用壓力基耦合絕對速度Simple算法，二階迎風離散格式，穩(wěn)態(tài)單相流。采用Realizable-湍流模型，考慮流場的旋轉以及彎曲壁面流動，邊界條件采用速度入口，壓力出口，出口壓力=0，壁面無滑移，殘差收斂設置為10×10。

表4 流體及工況參數

3.3 網格無關性檢驗

數值模擬過程中，幾何模型的網格劃分精度決定了計算過程的精確性和準確度。通常劃分的網格數量越多，計算結果的精度越高，但是往往受限于計算過程的時間成本和計算設備的硬件條件限制，網格劃分的方法、數量影響著模擬方法的選擇。如表5所示，分別以18 795、38 520、54 912、72 225四種網格劃分數量為對象，討論當入口速度為0.65、0.7、0.75 m/s時的阻力壓降，如圖7所示。不同速度條件下，阻力壓降均呈現出先升高后降低至穩(wěn)定波動，確定選取網格數量為54 912開展數值模擬分析。

表5 幾何結構網格無關性檢驗參數

圖7 網格獨立性檢驗過程

3.4 結果分析和討論

圖8所示為200目銅網在不同入口速度0.2、0.6、1.0、1.4 m/s條件下的截面的阻力壓降云圖。結合幾何模型結構，由于在距離入口0.2 m處設置為多孔躍遷層，可以看出在圓管區(qū)域內阻力壓降發(fā)生明顯降低趨勢，隨著入口速度的增加，區(qū)域內的阻力壓降為1 672.135、5 709.251、10 623.79、16 417.1 Pa，呈逐漸增大趨勢。

圖8 200目時不同流速下yoz截面壓降云圖

如圖9所示，不同尺寸銅網的阻力壓降實驗值和模擬值比較，其中，圖9(a)(b)(c)(d)分別顯示了實驗值和模擬值間標準差的誤差線情況?？梢钥闯觯寒斎肟谒俣容^低時，誤差線較小；隨著入口速度的增大，誤差線有逐漸增大趨勢。說明隨著速度的增大，阻力壓降與流體速度間逐漸偏離達西公式的線性關系。

圖9(e)(f)(g)(h)分別顯示了當置信水平為95%時，阻力壓降的實驗值和模擬值的置信區(qū)間擬合線，實驗值和擬合值的預測關系為1∶1。線性擬合得到的比例關系分別為1.168 64、1.106 21、1.083 27、0.853 18，得到的相關指數分別為99.821%、99.592%、99.474%、99.65%，說明線性擬合效果較好，阻力壓降的實驗值和模擬值吻合較好，通過實驗值求解的黏性阻力系數和慣性阻力系數具有較高的準確性，所設計的阻力系數測量裝置具有較好的適用性，同時所采用的多孔躍遷模型能夠準確模擬阻力壓降的變化情況，體現出較好的運算精度和收斂性。

圖9 不同尺寸銅網的阻力壓降實驗值和模擬值比較

圖10所示為速度沿流體域的軸心線分布情況，可以看出：速度在流動方向上距離入口0.2 m附近先發(fā)生大幅度下降。結合幾何模型結構，此位置為多孔躍遷區(qū)域，分別設置為不同孔徑規(guī)格的80目、120目、150目和200目銅網，通過局部放大圖可以明顯看出初始入口速度越大，銅網附近速度下降幅度越大，在通過銅網孔隙的過程中流通截面減小，流出銅網后會產生流速增大的趨勢，所產生的阻力壓降一部分用來克服黏性阻力和慣性阻力，另一部分轉化為流體通過銅網后的動能，且初始入口速度越大，所產生的速度梯度越大，符合阻力壓降與流體速度的變化規(guī)律。

圖10 不同尺寸銅網在流動方向的軸心線上速度變化曲線

如圖11所示，分別顯示200目銅網作為數值模擬多孔躍遷區(qū)域，入口速度為0.2、0.6、1.0、1.4 m/s時，沿流動方向截面上的速度場變化云圖?？梢钥闯觯涸诰嚯x入口0.2 m處速度場梯度波動明顯。結合前述流動方向軸心線上速度變化分析，各速度場云圖的梯度變化趨勢基本一致，在0.2 m處為多孔躍遷模型位置，速度場從入口處先增大至多孔躍遷區(qū)域，在多孔躍遷區(qū)域內產生阻力壓降，之后速度場先短暫減小，經過多孔躍遷穩(wěn)定區(qū)域后，流體速度場逐漸增大直至流動區(qū)域出口，而由軸心至壁面處的速度場呈梯度降低趨勢。

圖11 不同流速下200目銅網yoz截面速度云圖

4 結論

(1)流體在水平管段多孔介質區(qū)域內的流動阻力與流體的體積流量成反比，與流體在多孔介質區(qū)域內產生的阻力壓降Δ成正比。流體動力黏度系數越大，流動阻力越小，流動區(qū)域內截面尺寸越大，流體所受到的流動阻力越大。

(2) 利用無量綱π定理分析法，多孔介質的黏性阻力系數和慣性阻力系數模型與雷諾數和歐拉數相關，根據實驗值回歸得到量綱為一的量和待定系數。

(3)隨著銅網從80目增大至200目，單位長度上的孔隙增多，單位孔徑尺寸減小，孔隙率減小。銅網阻力壓降均隨入口速度的增大而增大，入口速度較小時，阻力壓降變化幅度較小，隨入口速度的增大，阻力壓降變化幅度逐漸增大；相同速度條件下，銅網阻力壓降隨著目數的增加而增大。

(4)通過銅網內流體阻力系數測量裝置得到的不同孔徑條件下阻力壓降實驗值，與利用多孔躍遷模型得到的模擬值吻合較好。多孔躍遷模型區(qū)域內產生的阻力壓降部分用來克服黏性阻力和慣性阻力，部分用來轉化為影響流體通過多孔躍遷區(qū)域后速度場逐漸增加的動能。