郜凱凱，陳寶明,2,3*，姜昊，蘆凱

(1. 山東建筑大學(xué) 熱能工程學(xué)院，山東 濟南 250101；2. 可再生能源建筑利用技術(shù)省部共建教育部重點實驗室，

?

基于LBM的通道內(nèi)多孔介質(zhì)交界面滑移效應(yīng)研究

郜凱凱1，陳寶明1,2,3*，姜昊1，蘆凱1

(1. 山東建筑大學(xué) 熱能工程學(xué)院，山東 濟南 250101；2. 可再生能源建筑利用技術(shù)省部共建教育部重點實驗室，

山東 濟南 250101；3. 山東省可再生能源建筑應(yīng)用技術(shù)重點實驗室，山東 濟南 250101)

多孔介質(zhì)與流體所構(gòu)成復(fù)合區(qū)域內(nèi)流體流動現(xiàn)象在自然界及社會許多行業(yè)之中廣泛存在，研究含多孔介質(zhì)通道內(nèi)流體流動的問題具有重要意義。文章基于格子Boltzmann方法對局部填充多孔介質(zhì)的通道內(nèi)流體流動進行模擬，研究了不同工況條件下對多孔介質(zhì)區(qū)域與純流體交界面處的滑移效應(yīng)變化規(guī)律，并采用編程進行數(shù)值模擬及結(jié)果驗證，分析了不同Re數(shù)和孔隙率對多孔介質(zhì)區(qū)域與純流體交界面處的滑移效應(yīng)。結(jié)果表明：速度滑移系數(shù)α始終為正，應(yīng)力滑移系數(shù)β始終為負值；速度滑移系數(shù)α和應(yīng)力滑移系數(shù)β的數(shù)值隨Re數(shù)和孔隙率ε增大時的變化趨勢不同，但滑移效應(yīng)的變化趨勢相同，即速度滑移效應(yīng)和應(yīng)力滑移效應(yīng)都隨Re數(shù)的增大而增強，隨孔隙率ε的增大而減弱。

通道；多孔介質(zhì)；格子Boltzmann方法；交界面；滑移效應(yīng)；

0　引言

多孔介質(zhì)與流體構(gòu)成復(fù)合區(qū)域內(nèi)流體流動現(xiàn)象廣泛存在于自然界和人類社會各行各業(yè)生產(chǎn)等工程均有極其重要的利用價值和實際意義。通道內(nèi)含多孔介質(zhì)的流體流動，通常涉及到其內(nèi)部的動量和質(zhì)量的傳遞現(xiàn)象，這些流動都是由外部因素所引起的。由于細通道的尺寸較小。其內(nèi)部有無多孔介質(zhì)及多孔介質(zhì)的情況，對流體的流動特性影響十分顯著。研究多孔介質(zhì)與純流體交界面處的流動變化規(guī)律變得十分重要。流體與多孔介質(zhì)交界處的速度滑移和應(yīng)力滑移一直以來是國內(nèi)外學(xué)者研究的重點。Beavers等提出在可滲透多孔介質(zhì)表面與流體空間交界面處存在速度不連續(xù)現(xiàn)象并通過實驗研究提出了速度滑移系數(shù)的半經(jīng)驗公式[1]，奠定了界面滑移效應(yīng)的研究基礎(chǔ)。Ochoa等利用修正Brinkman方程提出了交界面處應(yīng)力滑移條件[2-3]。Valdes等從解析解的角度分析了滲透率和孔隙率與界面應(yīng)力滑移相關(guān)[4]。伴隨對多孔介質(zhì)內(nèi)流體流動研究的深入，Baytas等對多孔介質(zhì)和流體交界面為水平及有臺階兩種情形分別進行了數(shù)值模擬分析[5]。Chen等和Liu等對方腔內(nèi)局部填充多孔介質(zhì)自然對流模型的交界面處滑移系數(shù)進行較為詳細分析和研究[6-7]。隨著數(shù)值模擬方法的不斷更新，格子Boltzmann方法(LBM)在計算流體力學(xué)領(lǐng)域得到迅速發(fā)展。而與傳統(tǒng)CFD模擬方法求解時間推進動量方程不同，LBM是基于粒子動力學(xué)展開的，利用介觀動理學(xué)方程—Boltzmann方程，通過Chapman-Enskog展開獲得宏觀運動方程。由于其計算簡單，并行性好等優(yōu)點，LBM被越來越多的被用來微小尺度的介質(zhì)進行模擬研究[8]。

文章采用格子Boltzmann方法對內(nèi)部局部填充多孔介質(zhì)通道的流體流動進行模擬。獲取了局部填充多孔介質(zhì)的通道內(nèi)流體流動變化的基本特征，探討不同Re數(shù)、孔隙率ε對通道內(nèi)流體流動的影響規(guī)律，重點研究不同Re數(shù)和孔隙率ε條件對多孔介質(zhì)與純流體交界面處滑移效應(yīng)的影響及其變化規(guī)律。

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

1.1物理模型的建立

建立的物理模型如圖1所示，在通道中間部分部分填充多孔介質(zhì)，上下邊界為固體壁面，左側(cè)為入口，入口速度為u，右側(cè)為出口。通過改變多孔介質(zhì)部分孔隙率ε(改變圓柱的直徑、位置)及流體的Re數(shù)獲得不同的工況。

圖1　問題的物理模型圖

1.2基本假設(shè)及數(shù)學(xué)模型

為了便于研究Re數(shù)、孔隙率ε等對部分填充多孔介質(zhì)通道內(nèi)流體及溫度場的影響，現(xiàn)作出以下4條基本假設(shè)：

(1) 流體區(qū)域為單一的純流體介質(zhì)，且區(qū)域內(nèi)部不存在任何障礙物。

(2) 流體中的粘性耗散項忽略不計且流體流動為二維流動。

(3) 通道內(nèi)流體為不可壓縮流體。

(4) 通道內(nèi)流體與多孔介質(zhì)骨架無任何化學(xué)反應(yīng)。

基于上述假設(shè)，文章采用格子Boltzmann方法模擬，LBM是應(yīng)用格子演化變量求解偏微分方程的一種十分活躍而極具發(fā)展前景的模擬方法，可以從描述微流體流動的連續(xù)Boltz-mann 方程得到。文章中粒子演化方程采用D2Q9格子Boltzmann模型(其離散速度如圖2所示)，模擬過程中使用格子單位(以下物理量符號沒有單位)進行模擬研究，其對應(yīng)的速度演化方程由式(1)[8-10]表示為

(1)

(2)

(3)

(4)

宏觀的流體密度ρ和速度u可通過式(5)、(6)分別計算得到，其中fα為粒子分布函數(shù), eα為方向α的格子速度。

(5)

(6)

平衡態(tài)函數(shù)可通過Chapman-Enskog展開方法得到D2Q9模型對應(yīng)的宏觀方程[11]，分別由式(7)、(8)表示為

(7)

=-p+·[ρμ(u+(u)T)]

(8)

1.3數(shù)學(xué)模型的邊界設(shè)置

該物理模型中邊界條件涉及通道壁面、出入口及多孔介質(zhì)壁面處速度邊界。其中，上下邊界為固體壁面，其速度邊界采用無速度滑移邊界條件，模擬中采用的是反彈與鏡面反射混合格式；多孔介質(zhì)壁面處采用反彈格式(bounce-back)，多孔介質(zhì)中ux=uy=0；出口處采用Zou-He邊界條件[12]；入口速度ux=0.1、uy=0(后續(xù)文章中u不加x、y下標(biāo)的，默認為ux)，計算區(qū)域格子解析度采用x、y方向上格子數(shù)分別為1400、200。

1.4滑移系數(shù)

多孔介質(zhì)骨架層與純流體之間存在交界面區(qū)域。而在該區(qū)域向多孔介質(zhì)區(qū)域過渡的微小空間，結(jié)構(gòu)較薄、影響機理復(fù)雜，因此兩區(qū)域交界面處動量傳遞較為復(fù)雜,如圖3所示。

隨著學(xué)者對該區(qū)域的研究深入，正確認識多孔介質(zhì)區(qū)域與流體區(qū)域交界面處的速度滑移與應(yīng)力滑移效應(yīng)成為探討流體在部分填充多孔介質(zhì)通道內(nèi)流動的關(guān)鍵。在通道的純流體區(qū)域下側(cè)，水平邊界層內(nèi)主要以水平流動為主。與水平速度相比，豎直方向速度較小，故在交界面處利用速度滑移和應(yīng)力滑移條件進行分析，并建立在交界面區(qū)域速度滑移系數(shù)α、應(yīng)力滑移系數(shù)β與速度u和剪切應(yīng)力du/dy之間的表達式[13]分別由式(9)、(10)表示為

圖3　多孔介質(zhì)骨架層與純流體構(gòu)成的復(fù)合域圖

(9)

(10)

式(9)、(10)經(jīng)變形，無量綱化后，獲得速度滑移系數(shù)α和應(yīng)力滑移系數(shù)β分別由式(11)、(12)表示為

(11)

(12)

2　數(shù)值模擬結(jié)果驗證

驗證文章數(shù)學(xué)模型及程序的正確性，這里選取雷諾數(shù)Re=10、20、40時的二維靜止圓柱繞流問題進行數(shù)值模擬，計算區(qū)域格子解析度采用x方向1200格子，y方向600格子，圓柱直徑d=30，上下為固體邊界，左側(cè)入口流速u=0.1，右側(cè)為出口，給出了計算結(jié)果，測量回流區(qū)的長度L′，計算阻力系數(shù)CD并與Nieuwstadt和He的結(jié)果分別進行了比較，其中,Re=ud/ν[14-15]。結(jié)果見表1，文章模型的與參考模型的計算很是接近(誤差率均在5%范圍內(nèi))，可知文章模型的正確性。

表1　穩(wěn)定時的流動參數(shù)對比(r為圓柱半徑)

3　數(shù)值模擬主要參數(shù)影響分析

3.1雷諾數(shù)Re的影響

Re數(shù)的變化對通道內(nèi)流體流動的影響十分明顯，它表征了粘滯力的強度。層流中，相比慣性力，粘滯力對流體流動的影響較大。不同的Re數(shù)條件對多孔介質(zhì)交界面處的速度滑移系數(shù)α與應(yīng)力滑移系數(shù)β的影響也十分重要。

圖4　不同Re數(shù)下滑移系數(shù)沿x方向的變化圖(孔隙率ε=0.8) (a)速度滑移系數(shù)α ；(b) 應(yīng)力滑移系數(shù)β

3.2 孔隙率ε的影響

孔隙率是多孔介質(zhì)的重要參數(shù)，孔隙率的變化直接影響多孔介質(zhì)內(nèi)部的流場的變化，進而影響純流體與多孔介質(zhì)交界面處的滑移現(xiàn)象。對比通道內(nèi)部不同孔隙率的多孔介質(zhì)區(qū)域與純流體區(qū)域交界面處的速度及應(yīng)力滑移系數(shù)，分析多孔介質(zhì)孔隙率對界面處的影響效應(yīng)。其中三種多孔介質(zhì)(ε=0.7、0.8和0.9)為不同直徑的圓柱排列而成，ε=0.9對應(yīng)的無量綱圓柱直徑為6×10-2，ε=0.8對應(yīng)的無量綱圓柱直徑為8.5×10-2，ε=0.7對應(yīng)的無量綱圓柱直徑為10.5×10-2。多孔介質(zhì)部分的滲透率k和Da數(shù)由式(13)表示[16]，其中L為特征長度(取1)。

(13)

由于采用圓柱構(gòu)造的多孔介質(zhì)區(qū)域，其Da數(shù)隨孔隙率ε的增大而增大，經(jīng)計算結(jié)果表明：在相同Re數(shù)條件下，速度滑移系數(shù)α隨孔隙率ε的增大而增大(如圖7(a)所示)，對應(yīng)的應(yīng)力滑移系數(shù)β隨孔隙率ε增大而減小(如圖7(b)所示)。

4　結(jié)論

純流體區(qū)域交界面處的滑移現(xiàn)象進行模擬分析，探討在不同Re數(shù)和不同多孔介質(zhì)的孔隙率ε下，交界面處滑移系數(shù)的變化規(guī)律，計算結(jié)果表明：

(1) 速度滑移系數(shù)α總為正值而應(yīng)力滑移系數(shù)β總為負值。

(2) 速度滑移效應(yīng)和應(yīng)力滑移效應(yīng)都隨Re數(shù)的增大而增強，隨孔隙率ε的增大而減弱。

(3) 當(dāng)多孔介質(zhì)區(qū)域處的空隙率ε一定時，交界面處的速度滑移系數(shù)α隨流體的Re數(shù)增大而增大，應(yīng)力滑移系數(shù)β隨著Re數(shù)的增大而減小。

(4) 當(dāng)流體的Re數(shù)不變時，鑒于增大圓柱構(gòu)造多孔介質(zhì)部分的孔隙率ε使得對應(yīng)的Da數(shù)也相應(yīng)的增大，經(jīng)計算可知：交界面處速度滑移系數(shù)α隨空隙率ε的增大而增大，應(yīng)力滑移系數(shù)β隨孔隙率ε的增大而減小。

圖5　不同孔隙率ε下交界面處當(dāng)量系數(shù)沿x方向的分布圖(Re=100)(a)系數(shù)α/Da0.5 ；(b)系數(shù)β/Da0.5

圖6　不同孔隙率ε下無量綱水平速度ux沿高度y的分布圖(Re=100)

圖7　不同孔隙率ε下滑移系數(shù)沿x方向的分布圖(Re=100)(a)速度滑移系數(shù)α ；(b)應(yīng)力滑移系數(shù)β

[1]Beavers G. S., Joseph D. D.. Boundary Conditions at a Naturally Permeable Wall[J]. Journal of Fluid Mechanics, 1967, 30(1):197-207.

[2]Ochoa-Tapia J. A., Whitaker S.. Momentum transfer at the boundary between a porous medium and a homogeneous fluid-I: theoretical development [J]. International Journal of Heat and Mass Transfer, 1995, 38(4): 2635-2646.

[3]Ochoa-Tapia J. A., Whitaker S. Momentum transfer at the boundary between a porous medium and a homogeneous fluid—II. Comp-arison with experiment[J]. International Journal of Heat & Mass Transfer, 1995, 38(14):2647-2655.

[4]Valdés-Parada F. J.,Ochoa-Tapia J. A., Aguilar-Madera C. G.,etal. Velocity and stress jump conditions between a porous medium and a fluid[J]. Advances in Water Resources, 2013(62): 327-339.

[5]Baytas A. C., Baytas A. F., Ingham D. B.,etal.. Double diffusive natural convection in an enclosure filled with a step type porous layer: Non-darcy flow [J]. International Journal of Thermal Sciences, 2009, 48(4):665-673.

[6]Liu F., Chen B., Wang L.. Experimental and numerical estimation of slip coefficient in a partially porous cavity[J]. Experimental Thermal & Fluid Science, 2013, 44(1):431-438.

[7]Chen B., Liu F.. Numerical study of slip effect at the porous media/liquid interface in an enclosure partially filled with a porous medium[C].Penyazkov:Heat Transfer Research. New York: Begell House Inc, 2016.

[8]陶文銓.傳熱與流動問題的多尺度數(shù)值模擬方法與應(yīng)用[M].北京：科學(xué)出版社,2008.

[9]何雅玲,王勇,李慶.格子Boltzmann方法的理論及應(yīng)用[M].北京：科學(xué)出版社,2009.

[10]Guo Z., Shu C.. Lattice Boltzmann Method and Its Applications in Engineering[M]. Singapore: World Scientific, 2013.

[11]Meng X., Guo Z.. Localized lattice Boltzmann equation model for simulating miscible viscous displacement in porous media[J]. International Journal of Heat and Mass Transfer, 2016(100): 767-778.

[12]Nazari M., Mohebbi R., Kayhani M. H.. Power-law fluid flow and heat transfer in a channel with a built-in porous square cylinder: Lattice Boltzmann simulation[J]. Journal of Non-Newtonian Fluid Mechanics, 2014(204):38-49.

[13] Loussif N., Orfi J.. Slip flow heat transfer in micro-tubes with viscous dissipation[J]. Desalination and Water Treatment, 2015, 53(5): 1263-1274.

[13]Nieuwstadt F., Keller H. B.. Viscous flow past circular cylinders[J]. Computers & Fluids, 1972, 1(1):59-71.

[14]He X., Doolen G.. Lattice Boltzmann Method on curvilinear coordinates system: flow around a circular cylinder [J]. Journal of Computational Physics, 1997, 134(2):306-315.

[15]Nield D. A., Bejan A.. Convection in Porous Media[M]. New York: Springer, 2013.

(學(xué)科責(zé)編：吳芹)

Study on slippage effect of porous media interface within the channel by Lattice Boltzmann Method

Gao Kaikai1, Chen Baoming1,2,3*, Jiang Hao1,etal.

(1. School of Thermal Energy Engineering, Shandong Jianzhu University, Jinan 250101, China; 2. Key Laboratory of Renewable Energy Utilization Technologies in Building, Ministry of Education, Jinan 250101, China; 3. Shandong Key Laboratory of Renewable Energy Applied Technologies, Jinan 250101, China)

Fluid flow phenomenon on the composite area of porous media and fluid is so widespread in nature and human society. The study on the flow in porous media channels has important practical significance. The model on channel partially filled with structured porous layer was simulated numerically by lattice Boltzmann method. Study on the change rule of porous media area and pure fluid interface slippage effect under different conditions, numerical simulation by programming, analyzing principally the impact of different Re number and the porosity ε of the porous medium on the slip effect at the porous/fluid interface. The numerical results show that the velocity slip coefficient is always positive. On the contrary, the stress slip coefficient is always negative; the change trends of velocity slip and stress slip coefficient with Re and the porosity ε are different, instead, the tendency of the slip effect is the same, which means that the velocity slip effect and the stress slip effect are enhanced with the increase of Re number, and decreases with the increase of the porosity ε.

channel; porous media; Lattice Boltzmann Method; porous/fluid interface; slippage effect;

2016-03-04

國家自然科學(xué)基金資助項目(51076086)

郜凱凱(1990-)，男，在讀碩士，主要從事多孔介質(zhì)內(nèi)流動換熱等方面的研究.E-mail: wsgaokai0917@163.com

*:陳寶明( 1963-) ，男，教授，博士，主要從事復(fù)雜體系中的傳熱傳質(zhì)等方面的研究.E-mail: chenbm@sdjzu.edu.cn

1673-7644(2016)03-0237-07

TU996

A