冪律流體飽和多孔介質平板通道流動特性分析

田興旺，王平，徐士鳴

(1.大連理工大學能源與動力學院，遼寧 大連116024;2.大連海洋大學海洋與土木工程學院，遼寧大連116023)

冪律型非牛頓流體在多孔介質中的流動有著廣泛而重要的工程應用背景［1-4］，如化工工程中填充塔的性能優化、資源工程中的石油熱采等。由于冪律流體在多孔介質中的非達西流流動方程存在著很強的非線性效應，很難得到精確解，除了采用數值模擬方法［5-8］外，有學者借助積分的近似方法來求解［9-12］。如Nakayama等［9］對冪律流體飽和多孔介質的強迫對流問題進行了研究，運用邊界層動量積分方法求解了Darcy-Brinkman-Forchheimer流動模型的動量主控方程，但積分采用了簡單的線性速度分布模型。在文獻［9］基礎上，Nield等［10］采用同樣積分方法進一步研究了常壁溫條件下兩平板構成的冪律流體飽和多孔介質通道內強迫對流的熱發展情況。Chen等［11］通過邊界層動量積分方法研究了冪律流體在多孔介質通道中的強迫對流換熱過程，與文獻［9］不同的是積分采用拋物線速度分布模型，得到了含有無限級數的邊界層厚度的隱式表達式，但求解過程較復雜。而Thayalan等［12］也采用邊界層拋物線速度分布模型，運用邊界層動量積分方法，對流動充分發展的冪律流體動量主控方程進行求解，得到了邊界層厚度的顯式表達式及壓力梯度的表達式，但沒有考慮非達西流區的慣性效應。

綜上，本文在前人研究的基礎上，采用應用更廣泛的Darcy-Brinkman-Forchheimer流動模型，通過比較邊界層內線性和拋物線2種速度分布經驗關系式，運用邊界層動量積分方法更全面的分析了冪律流體飽和多孔介質通道的流動特性問題。得到了邊界層厚度顯式表達式及軸向速度和壓力降的無量綱表達式，并進一步討論了一些重要的無量綱參數，如粘性比、達西數、綜合慣性參數、孔隙率等對不同流變指數流體流動速度分布及壓力降的影響。

1 物理數學模型

如圖1所示，研究對象為一個填充各向同性多孔介質的板間距為2H的平行板通道流動問題。考慮二維流動充分發展的穩態情況，假設流體為單相、不可壓縮的冪律流體，除了流體粘性外其他流體參數物性是定值。

圖1 多孔介質平板通道流動模型示意圖Fig.1 Schematic diagram of flow model in a parallel plate channel filled with porous media

冪律流體Darcy-Brinkman-Forchheimer流動模型的動量主控方程為

式中:u為多孔介質達西流速，m/s;μeff為冪律流體的有效稠度系數，N·s/m2;μ為冪律流體的稠度系數，N·s/m2;φ為多孔介質孔隙率，n為冪律指數，K為修正的多孔介質滲透率，mn+1;b為Forchheimer慣性系數，m－1;dp/dx為壓力梯度，kg/(m2·s2)，是個定值。

定義無量綱的變量如下:無量綱坐標Y=y/H，無量綱粘度比為M=μeff/μ，特征參考速度為ud=(Kdp/(μdx))1/n，基于特征參考速度的無量綱流速為U=u/ud，基于特征參考速度修正的Forchheimer慣性參數為F=ρbK(ud)2－n/μ，修正的達西數為Da=(K/φn)2/(n+1)/H2。

把上述無量綱參數帶入方程式(1)得

對應的無量綱邊界條件:

按照文獻［12］的分析，粘性效應被限制在貼近板壁處的薄邊界層δ內，其速度梯度滿足方程(2)，而邊界層外為流核區，其無量綱的流核速度Uc滿足方程式(4)，從而使得流動呈現結構流的特性。

對式(2)在半個通道斷面上采用邊界層動量積分得

相應的邊界條件方程式(3)變為

對式(5)封閉求解，需要在邊界層中增加一個假想的速度分布經驗關系式。

假定邊界層內速度線性分布，則無量綱速度分布表達式為

將速度線性分布式(7)代入到式(5)，考慮邊界條件(6)，得到邊界層無量綱厚度值為

假定邊界層內速度二次拋物線分布，則無量綱速度分布表達式為

把拋物線速度分布表達式(9)帶入到方程(5)，考慮邊界條件(6)，得到邊界層無量綱厚度值為

其中，Γ為伽馬函數，給定一個冪律指數n，即可求出Γ(n)。

在這里定義體積平均速度um和無量綱體積平均速度Um如下:

當把式(11)代入式(7)和(8)得

將式(11)代入式(9)、(10)得

其中，F'=ρbK(um)2－n/μ ，是基于體積平均速度修正的Forchheimer慣性參數，也稱為綜合慣性參數，U'c是基于體積平均速度的無量綱流核速度。

根據式(1)可得到冪律型流體在多孔介質中流動的壓降關系式為

定義無量綱的范寧(Fanning)摩擦系數和修正的無量綱雷諾數Re分別為

將式(17)、(18)代入到式(16)中，可得到無量綱軸向壓力降方程式為

綜上，已知冪律指數n、孔隙度φ，達西數Da及Forchheimer參數F'，聯立式(12)、(13)或者(14)、(15)，采用迭代求解方法就能得到邊界層厚度值δ和無量綱流核速度值U'c。然后代入到式(12)或(14)及(19)即可得到平板通道斷面相應的流速分布曲線和壓力降的具體結果。同時可以看出，邊界層厚度值一方面與流體的冪律指數n和多孔介質的達西數Da密切相關，另一方面也隨著流動的綜合慣性參數F'的大小而變化。

對于滲透率較低的多孔介質，在流速較低時慣性效應很小，近似F'=0，此時，速度線性分布模型和拋物線分布模型所對應的邊界層厚度分別為

而對于滲透率較高的多孔材料，在流動慣性效應較大時，其邊界層厚度和壓力降與達西數Da和綜合慣性參數F'均有關。

2 拋物線模型與線性模型的比較

為了證實二次拋物線速度分布和線性速度分布的有效性，將這2個模型所得的結果分別與分析解作比較。在流體為牛頓流體(n=1)、無量綱粘度比M=1和綜合慣性參數F'=0的情況下，對式(2)化簡得

基于特征速度和平均速度的分析解表達式分別為

圖2描述了牛頓流體(n=1)在粘度比M=1、綜合慣性參數F'=0時，比較了不同達西數Da情況下，拋物線和線性2種模型得出的速度分布與分析解接近程度。由圖中可知，在Da較小時，拋物線和線性2種模型所對應的速度分布結果與分析解吻合度較好，這是由于較小的達西數使得邊界粘性效應降低，流核區域變大，速度剖面曲線變得更平坦。這說明了在達西流區(一般可認為Da＜0.001)，邊界層內假定流速分布關系式的選取對計算結果影響很小。

隨著Da的增大，可以看出在邊界層的鄰近區域內，2種模型與分析解存在的差異也變得明顯，但是拋物線模型顯然比線性模型所得的速度分布結果有一個更好的近似，說明了在非達西流區，邊界層內速度分布采用拋物線模型得到的結果精度較高。值得注意的是拋物線模型在本質上也更符合實際速度的分布形式，因為邊界層厚度處不可能出現速度的突變。

圖2 拋物線、線性模型與分析解的比較(n=1，M=1，F'=0)Fig.2 The comparison of the parabolic and linear velocity distribution models with the analysis solution for the case(n=1，M=1，F'=0)

3 計算結果分析

3.1 達西數對流動邊界層厚度的影響

圖3描述了在粘性比M=1及綜合慣性參數F'=100的條件下，邊界層厚度δ在不同冪律指數n下隨達西數Da的變化情況。

圖3 達西數對邊界層厚度的影響(M=1，F'=100)Fig.3 The effect of the Darcy number on the boundary layer thickness for the case(M=1，F'=100)

從圖3可以看出，隨著達西數Da的減小，邊界層厚度δ變薄，說明了邊界粘性效應被限制在貼近壁面處的一薄層內。隨著冪律指數n的增大，邊界層厚度δ變厚，這一點不難理解，因為冪律型非牛頓流體在無多孔介質的通道內流動時，剪切變稠流體(n=1.5)比剪切變稀流體(n=0.5)的速度剖面更像拋物線。此外，從圖3上還可以看出采用拋物線模型所得到的邊界層厚度要比線性模型厚一些。

3.2 無量綱軸向流速分布的主要影響因素

3.2.1 達西數的影響

圖4描述了在粘性比M=1及綜合慣性參數F'=100的條件下，采用拋物線模型所得的無量綱軸向速度分布隨達西數Da的變化情況，同時考慮冪律指數n的影響。從圖4可知Da的變化對剪切變稠流體(n=1.5)的影響比剪切變稀流體(n=0.5)大，即隨著Da的減小，剪切變稠流體的流核區域變化較明顯，速度剖面變得平坦，壁面處邊界層流動特征更加顯著。這是因為多孔介質的固體基質對流體流動的阻擋作用對剪切變稠流體的影響更明顯。

圖4 達西數對無量綱軸向流速分布的影響(M=1，F'=100)Fig.4 The effect of the Darcy number on the dimensionless axial velocity distribution for the case(M=1，F'=100)

3.2.2 綜合慣性參數的影響

圖5描述了在粘性比M=1及達西數Da=0.01時，采用拋物線模型所得的無量綱軸向速度分布隨F'數的變化情況，同時考慮冪律指數n的影響。從圖5可知，綜合慣性參數F'數的變化對剪切變稀流體(n=0.5)的影響比剪切變稠流體(n=1.5)大，即隨著F'數的增大，剪切變稀流體的邊界層更薄，邊界層內速度梯度更大，速度剖面也變得更平坦，使得流核區域變化較明顯，壁面處邊界層流動特征更顯著。

圖5 綜合慣性參數對無量綱軸向流速分布的影響(M=1，Da=0.01)Fig.5 The effect of the integrated inertial parameter on the dimensionless axial velocity distribution for the case(M=1，Da=0.01)

3.2.3 粘性比的影響

圖6描述了在綜合慣性參數F'=100及達西數Da=0.01時，采用拋物線模型所得的無量綱流核速度分布隨粘性比M和冪律指數n的變化情況。由圖6可知，隨著M的增大，流核速度U'略微增大，但變化很小;此外，隨著冪律指數n的增大，流核速度有所增加，但增加量也不大，綜上可以認為無量綱粘性比對冪律型非牛頓的流核速度的影響較小，即對無量綱軸向流速分布的影響可以忽略。

圖6 粘性比對無量綱軸向流速分布的影響(F'=100，Da=0.01)Fig.6 The effect of the viscosity ratio on the dimensionless axial velocity distribution for the case(F'=100，Da=0.01)

3.3 無量綱軸向壓力降的主要影響因素

3.3.1 達西數和綜合慣性指數的影響

圖7是在M=1，φ =0.5時，冪律流體在多孔介質中流動無量綱軸向壓力降隨著達西數Da和綜合慣性參數F'的變化關系，同時考慮了冪律指數n的影響。從圖7可知，在F'數較小時，流體流動阻力趨近于Darcy流模型(F'=0)所預測出的軸向壓力降方程為常數的情況;隨著F'的增大，流體流動從Darcy流轉變成Forchheimer慣性流，壓力降迅速上升;在慣性效應較強時，壓力降隨綜合慣性參數F'成線性增長。同時可以看出，多孔介質的特性和流體的流變特性對壓力降有著較大的影響，即隨著滲透率的減小(Da為0.01～0.000 1)和冪律指數n的增大，壓力降也是顯著增大的。

圖7 達西數和綜合慣性指數對無量綱軸向壓力降的影響(M=1，φ =0.5)Fig.7 The effect of the Darcy number and integrated inertial parameter on the dimensionless pressure drop for the case(M=1，φ =0.5)

3.3.2 孔隙率和粘性比對無量綱軸向壓力降的影響

圖8是在F'=100，Da=0.01時，冪律流體在多孔介質中流動無量綱軸向壓力降隨著孔隙率φ和粘性比M的變化關系，同時考慮了冪律指數n的影響。從圖上可知，在φ較小時，流體流動阻力很大，隨著φ的增大，壓力降迅速降低，但減小變緩，特別是剪切變稀流體(n=0.5)，隨著φ的增大，壓力降變化不明顯。還可看出隨著M數的增大壓力降也是逐漸下降的，其中剪切變稠流體(n=1.5)較其他2種流體隨著φ的增大，壓力降降低更為明顯。

圖8 孔隙率和粘性比對無量綱軸向壓力降的影響(F'=100，Da=0.01)Fig.8 The effect of the porosity and viscosity ratio on the dimensionless pressure drop for the case(F'=100，Da=0.01)

4 結論

1)推導出動量主控方程及相應邊界條件的無量綱表達式，結合假定的邊界層內速度分布經驗關系式，得出邊界層厚度顯式表達式及無量綱的軸向流速和壓力降的計算表達式。

2)在Da較小時，拋物線和線性2種模型所對應的速度分布結果與分析解吻合度較好，邊界層內假定流速分布關系式的選取對計算結果影響很小。隨著Da的增大，拋物線模型比線性模型所得的速度分布結果有一個更好的近似，邊界層內速度分布采用拋物線模型得到的結果精度較高。

3)流動邊界層隨著達西數Da的減小而變薄，隨著冪律指數n的增大而變厚;達西數Da的變化對剪切變稠流體的影響較大，而綜合慣性指數F'的變化對剪切變稀流體的影響較大;無量綱粘性比M對冪律型非牛頓的流速分布的影響較小。

4)在Da較小時，流體流動阻力很大，隨著Da的增大，壓力降先迅速降低，但之后降低減緩;隨著F'及冪律指數n的增大，流動壓力降迅速上升;隨著孔隙率φ和粘性比M數的增大，壓力降逐漸降低。

［1］HARRIS S D，INGHAM D B，POP I.Mixed convection boundary-layer flow near the stagnation point on a vertical surface in a porous medium:Brinkman model with slip［J］.Transport in Porous Media，2009，77:267-285.

［2］CHEN H T，CHEN C K.Natural convection of a non-Newtonian fluid about a horizontal cylinder and sphere in a porous medium［J］.International Communications in Heat and Mass Transfer，1988，15:605-614.

［3］SHENOY A V.Darcy-Forchheimer natural，forced and mixed convection heat transfer in non-Newtonian power-law fluid-saturated porous media［J］.Transport in Porous Media，1993，11(3):219-241.

［4］CHENG C Y.Soret and Dufour effects on free convection boundary layers of non-Newtonian power law fluids with yield stress in porous media over a vertical plate with variable wall heat and mass fluxes［J］.International Communications in Heat and Mass Transfer，2011，38:615-619.

［5］PANTOKRATORAS A，MAGYARI E.Forced convection flow of power-law fluids over a flat plate embedded in a Darcy-Brinkman porous medium［J］.Transport in Porous Media，2010，85(1):143-155.

［6］HADIM H.Non-Darcy natural convection of a non-Newtonian fluid in a porous cavity［J］.International Communications in Heat and Mass Transfer，2006，33(10):1179-1189.

［7］CHEN G M，TSO C P.Effects of viscous dissipation on forced convective heat transfer in a channel embedded in a power-law fluid saturated porous medium［J］.International Communications in Heat and Mass Transfer，2011，38(1):57-62.

［8］EL-AMIN M F，SALAMA A，et al.A conditionally stable scheme for a transient flow of a non-newtonian fluid Saturating a porous medium［J］.Procedia Computer Science，2012，9:651-660.

［9］NAKAYAMA A，SHENOY A V.Non-Darcy forced convective heat transfer in achannel embedded in a non-Newtonian inelastic fluid-saturated porous medium［J］.Can J Chem Eng，1993，71:168-173.

［10］NIELD D A，KUZNETSOV A V.Thermally developing forced convection in a channel occupied by a porous medium saturated by a non-Newtonian fluid［J］.International Journal of Heat and Mass Transfer，2005，48(6):1214-1218.

［11］CHEN G，HADIM H A.Forced convection of a power-law fluid in a porous channel—integral solutions［J］.Journal of Porous Media，1999，2(1):59-70.

［12］THAYALAN N，HUNG Y M.Momentum integral method for forced convection in thermal nonequilibrium power-law fluid-saturated porous media［J］.Chemical Engineering Communications，2013，200(2):269-288.