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

田興旺，王平，徐士鳴

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

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

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

1 物理數(shù)學(xué)模型

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

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

冪律流體Darcy-Brinkman-Forchheimer流動(dòng)模型的動(dòng)量主控方程為

式中:u為多孔介質(zhì)達(dá)西流速，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)，是個(gè)定值。

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

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

對應(yīng)的無量綱邊界條件:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

隨著Da的增大，可以看出在邊界層的鄰近區(qū)域內(nèi)，2種模型與分析解存在的差異也變得明顯，但是拋物線模型顯然比線性模型所得的速度分布結(jié)果有一個(gè)更好的近似，說明了在非達(dá)西流區(qū)，邊界層內(nèi)速度分布采用拋物線模型得到的結(jié)果精度較高。值得注意的是拋物線模型在本質(zhì)上也更符合實(shí)際速度的分布形式，因?yàn)檫吔鐚雍穸忍幉豢赡艹霈F(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 計(jì)算結(jié)果分析

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

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

圖3 達(dá)西數(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可以看出，隨著達(dá)西數(shù)Da的減小，邊界層厚度δ變薄，說明了邊界粘性效應(yīng)被限制在貼近壁面處的一薄層內(nèi)。隨著冪律指數(shù)n的增大，邊界層厚度δ變厚，這一點(diǎn)不難理解，因?yàn)閮缏尚头桥ｎD流體在無多孔介質(zhì)的通道內(nèi)流動(dòng)時(shí)，剪切變稠流體(n=1.5)比剪切變稀流體(n=0.5)的速度剖面更像拋物線。此外，從圖3上還可以看出采用拋物線模型所得到的邊界層厚度要比線性模型厚一些。

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

3.2.1 達(dá)西數(shù)的影響

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

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

圖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及達(dá)西數(shù)Da=0.01時(shí)，采用拋物線模型所得的無量綱流核速度分布隨粘性比M和冪律指數(shù)n的變化情況。由圖6可知，隨著M的增大，流核速度U'略微增大，但變化很小;此外，隨著冪律指數(shù)n的增大，流核速度有所增加，但增加量也不大，綜上可以認(rè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 達(dá)西數(shù)和綜合慣性指數(shù)的影響

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

圖7 達(dá)西數(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時(shí)，冪律流體在多孔介質(zhì)中流動(dòng)無量綱軸向壓力降隨著孔隙率φ和粘性比M的變化關(guān)系，同時(shí)考慮了冪律指數(shù)n的影響。從圖上可知，在φ較小時(shí)，流體流動(dòng)阻力很大，隨著φ的增大，壓力降迅速降低，但減小變緩，特別是剪切變稀流體(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)推導(dǎo)出動(dòng)量主控方程及相應(yīng)邊界條件的無量綱表達(dá)式，結(jié)合假定的邊界層內(nèi)速度分布經(jīng)驗(yàn)關(guān)系式，得出邊界層厚度顯式表達(dá)式及無量綱的軸向流速和壓力降的計(jì)算表達(dá)式。

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

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

4)在Da較小時(shí)，流體流動(dòng)阻力很大，隨著Da的增大，壓力降先迅速降低，但之后降低減緩;隨著F'及冪律指數(shù)n的增大，流動(dòng)壓力降迅速上升;隨著孔隙率φ和粘性比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.