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

田興旺，王平，徐士鳴

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

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

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

1 物理數(shù)學模型

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

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

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

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

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

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

對應的無量綱邊界條件:

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

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

相應的邊界條件方程式(3)變?yōu)?/p>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 拋物線、線性模型與分析解的比較(n=1，M=1，F(xiàn)'=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(xiàn)'=0)

3 計算結(jié)果分析

3.1 達西數(shù)對流動邊界層厚度的影響

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

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

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

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

3.2.1 達西數(shù)的影響

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

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

3.2.2 綜合慣性參數(shù)的影響

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

圖5 綜合慣性參數(shù)對無量綱軸向流速分布的影響(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描述了在綜合慣性參數(shù)F'=100及達西數(shù)Da=0.01時，采用拋物線模型所得的無量綱流核速度分布隨粘性比M和冪律指數(shù)n的變化情況。由圖6可知，隨著M的增大，流核速度U'略微增大，但變化很小;此外，隨著冪律指數(shù)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 達西數(shù)和綜合慣性指數(shù)的影響

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

圖7 達西數(shù)和綜合慣性指數(shù)對無量綱軸向壓力降的影響(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時，冪律流體在多孔介質(zhì)中流動無量綱軸向壓力降隨著孔隙率φ和粘性比M的變化關系，同時考慮了冪律指數(shù)n的影響。從圖上可知，在φ較小時，流體流動阻力很大，隨著φ的增大，壓力降迅速降低，但減小變緩，特別是剪切變稀流體(n=0.5)，隨著φ的增大，壓力降變化不明顯。還可看出隨著M數(shù)的增大壓力降也是逐漸下降的，其中剪切變稠流體(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 結(jié)論

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

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

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

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

［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.