徐芳芳，付海明，雷澤明，甘 靈

（1.東華大學 環(huán)境科學與工程學院，上海 201620；2.中建三局建設工程股份有限公司，湖北 武漢 430000）

壓力損失是纖維過濾介質的重要參數(shù)之一.當纖維多孔介質內的流體流速較低或者雷諾數(shù)（Re）較小時，一般假設其流場服從經(jīng)典的達西定律，即壓力梯度和平均速度存在線性關系.隨著Re增加到一極限值時，其不再服從達西定律，多孔介質內的流體流動由于慣性及黏性等作用呈現(xiàn)非線性趨勢[1].VAFAI[2]是較早研究達西定律擴展現(xiàn)象的研究者之一，其使用宏觀平均技術來研究達西定律中的摩擦阻力和慣性阻力，認為這些非線性的影響使得壓力損失比達西定律得到的要高.基于剛性多孔介質的一維穩(wěn)態(tài)不可壓縮層流來求解Navier-Stokes（N-S）方程，文獻[3]提出一個包括非達西影響的修正表達式.DEIBER等[4]考慮了黏性和慣性影響，采用渦量流方程，研究了高Re下多孔介質的自然對流.文獻[5]通過成堆馬鈴薯的壓力損失－速度關系試驗，驗證了Forchheimer方程的有效性，并利用有限元素 方 法 對 三 維 的 Darcy-Forchheimer-Brinkman（DFB）模型進行數(shù)值求解，結果表明該模型的計算結果與壓力損失的試驗數(shù)據(jù)一致.文獻[6]基于達西拓展方程，使用有限體積方法，模擬了多孔介質表面存在液滴時的影響.YANG等[7]基于DFB方法，模擬了飽和多孔介質內的流體流動，得出Re與速度一樣，都是影響壓力損失的重要因素.但有些研究者認為，F(xiàn)orchheimer方程僅在高Re時有效，而在達到完全發(fā)展的紊流階段前，還應該存在過渡流.MEI等[8]通過嚴密的理論推導，給出了各向同性介質壓力損失的三次方表達式.ANDRADE等[9]通過對試驗數(shù)據(jù)的整理，采用速度的三次方對Forchheimer方程進行修正，用以描述非達西現(xiàn)象.CAI等[10]利用波爾茲曼格子法，給出了二維簡化多孔介質的三次方及過渡區(qū)的表達式.雖然這些研究結果對于非達西現(xiàn)象的研究很有應用價值，但其都是建立在一維或二維的簡單多孔介質模型上的，這對真實的纖維多孔介質而言，并不適用.因為纖維多孔介質的流場發(fā)生在三維空間，且纖維一般是隨機排列的.另外，F(xiàn)orchheimer方程中的系數(shù)理論上難以給出，一般需通過試驗數(shù)據(jù)得到，且給出值通常隨不同研究者或不同的多孔介質模型而不同[11].隨著計算機技術的發(fā)展與計算能力的提高，非達西流的數(shù)值研究得到較大的發(fā)展.本文基于VBA程序語言，創(chuàng)建了一系列三維隨機排列纖維的多孔介質模型，采用計算流體動力學（CFD）軟件模擬了孔隙尺度下纖維多孔介質模型內部的氣相流場，并將模擬數(shù)據(jù)與一些經(jīng)驗關系式得到的結果進行比較，最后通過數(shù)據(jù)擬合得到一個可以描述非達西現(xiàn)象的表達式.

1 多孔介質的經(jīng)驗模型

多孔介質的不可壓縮流可以通過兩種方法來描述：一種是孔隙尺度下求解N-S方程；另一種是直接求解平均方程或已存在的經(jīng)驗方程式[10].由于多孔介質幾何的復雜性和不規(guī)律性，直接在孔隙尺度下求解N-S方程式，計算量大且很難得到精確的計算結果.通常對多孔介質進行簡化，忽略其復雜性，得出相應的經(jīng)驗公式并應用于工程方面.但這些簡化模型往往是基于經(jīng)驗假設，并不適合用來研究多孔介質內復雜的流動情況.事實上，這兩種方法之間并不是相互獨立的.

由于Re與速度都是影響壓力損失的重要參數(shù)[10].故可將廣泛使用的Re作為標準，直觀地研究非達西現(xiàn)象中的流動情況，其中，Re為慣性力與黏性力的比值.

其中：ρ為流體密度；為流體平均速度；μ為流體動力學黏度；Df為纖維直徑，在本文中作為特征長度.假設纖維多孔介質流動分為線性達西流區(qū)以及非線性流區(qū)兩個不同的區(qū)域，分別研究其流動特性.

（1）線性達西流區(qū)：在這一區(qū)域內，流場流動符合達西定律，其壓力梯度－速度表達式如式（2）所示.

其中：Δp為流體的壓力損失；K為滲透率；L為過濾介質厚度.為了直觀地表述慣性的影響，轉變式（2）為式（3）所示.

其中：Ke為有效滲透率；c為常數(shù).在線性達西流區(qū)域內，黏性力占主導地位.

（2）非線性流區(qū)：在這一區(qū)域內，隨著Re不斷增加，多孔介質的流動可用Forchheimer方程進行描述[3].

其中：b1為非達西流系數(shù)，其值與孔隙尺寸、形狀和孔隙率有關；方程右邊的第一、二項分別為黏性力和慣性力.對式（4）進行變換得如式（5）所示的Forchheimer方程.

其中：b′1＝b1Df，為一無量綱參數(shù).在非線性流區(qū)域內，慣性力影響比黏性力影響大，流體行為主要受慣性力支配.

2 數(shù)值模擬計算

2.1 三維隨機纖維過濾介質模型的建立

本文中基于VBA語言開發(fā)技術和AutoCAD平臺，采用隨機算法生成三維纖維多孔介質.在模型生成過程中，長方體代表纖維多孔介質材料，細圓柱體代表纖維.通過輸入恰當?shù)膮?shù)（如纖維的長、寬及纖維的直徑、數(shù)量），即可得到相應孔隙率的模型.需要注意的是，為了生成高質量的網(wǎng)格，模型中不允許纖維間相互接觸.因此，在模型生成過程中，忽略纖維的彎曲、壓延及生產(chǎn)過程中造成的變形.

2.2 模型的網(wǎng)格劃分及邊界條件

圖1給出了本文的計算模型及相應的邊界條件.假設氣體通過速度入口進入計算區(qū)域，然后從壓力出口流出計算區(qū)域.在本文中，為了消除纖維面較強流速或壓力梯度的影響，入口邊界和出口邊界分別位于距離過濾器上游10倍于纖維直徑和下游5倍于纖維直徑處[12].由于流體流動主要在纖維厚度方向，側邊的邊界條件對模擬結果的影響小，因此，計算域的側邊采用對稱性邊界條件.另外，對于纖維表面的流體，本文假設為無滑移邊界.

圖1 數(shù)值模擬區(qū)域及邊界條件Fig.1 Numerical simulation area and boundary conditions

在給出壓力損失的模擬結果前，需確保模擬結果獨立于所劃分的網(wǎng)格.為此，本文計算了纖維多孔介質的填充率（α）為4.51%的模型在不同的纖維邊界網(wǎng)格尺度數(shù)及總網(wǎng)格數(shù)量下的壓力損失.改變纖維周邊的網(wǎng)格尺度數(shù)，將其由9開始減小至2，得網(wǎng)格尺度數(shù)對壓力損失的影響結果如圖2所示.由圖2可以看出，降低網(wǎng)格尺度數(shù)開始會導致壓力損失增加，至一定程度（纖維網(wǎng)格邊界尺度數(shù)小于2.5）后，壓力損失幾乎沒有變化.

總網(wǎng)格數(shù)對壓力損失的影響如圖3所示.由圖3可以看出，增加總網(wǎng)格數(shù)量，開始時壓力損失逐步增加，但當總網(wǎng)格數(shù)超過2×106時，壓力損失幾乎

圖2 纖維邊界網(wǎng)格尺度對壓力損失的影響Fig.2 Effect of fiber boundary mesh size on pressure drop

沒有變化.因此，在本文的數(shù)值模擬計算中，所有模型的計算工況均選擇網(wǎng)格尺度數(shù)為2.5和總網(wǎng)格數(shù)為2×106，以保證計算結果的準確性.另外，本文涉及的所有多孔介質模型中的纖維直徑為16μm，過濾介質厚度為5mm；纖維介質尺寸為500μm×500 μm×200μm.

圖3 纖維邊界網(wǎng)格數(shù)量對壓力損失的影響Fig.3 Effect of fiber boundary mesh numbers on pressure drop

3 模擬結果與討論

3.1 流體流動分析

根據(jù)模擬計算，得出不同α下隨機排列的三維纖維多孔介質的壓力損失與流體平均速度的關系如圖4所示.同樣，為了更好地觀察到慣性的影響，將壓力梯度－流體平均速度關系轉化為有效滲透率與Re的關系（如式（5）），并繪出D2f／Ke與Re的關系曲線如圖5所示.從圖4和5可以清楚地看出，隨著流體平均速度或Re的增加，纖維過濾多孔介質內的壓力損失逐漸增加，但這種增加是非線性的，因此，單純用達西定律來描述該模型下的流場是不合適的.值和由式（5）得出的K的計算值，發(fā)現(xiàn)所得值有些不同，但是兩者之間的差異非常小，最大相對誤差不超過4%.因此，式（5）可以精確地描述非線性流區(qū)的流動規(guī)律.

選擇其中一種模型（α＝6.76%）的模擬結果（如圖6所示），具體分析纖維過濾介質有效滲透率與Re的關系.由圖6（a）可以看出，當Re≤0.33（對應的流體平均速度為0.3m／s）時，有效滲透率與Re的比值幾乎為一常數(shù).數(shù)值模擬結果與達西定律表達式非常吻合，該模擬區(qū)域為線性達西流區(qū).由圖6（b）可以看出，Re足夠大時，其與纖維過濾介質有效滲透率之間的比值不再是常數(shù)，即數(shù)值模擬結果不再符合達西定律.將數(shù)值模擬結果與Forchheimer方程對比發(fā)現(xiàn)，模擬結果可以用Forchheimer方程來預測，該模擬區(qū)域為非線性流區(qū).

分別比較由式（3）推算出的不同α下滲透率K

3.2 壓力損失－速度關系式

分別模擬計算5種纖維多孔介質填充率（α分別為2.81%，3.36%，3.82%，4.51%，6.76%）和7種迎面風速下（分別為0.008 3，0.050 0，0.100 0，0.150 0，0.200 0，0.250 0，0.300 0m／s）的流場分布情況，并依據(jù)這35組數(shù)據(jù)結果進行分析，以方程K＝D2f／[aα3／2（1＋bα3）]為基礎，a和b為待定系數(shù)，可得出任意分布纖維多孔介質滲透率表達式如式（6）所示[13].

結合數(shù)據(jù)擬合軟件，得出系數(shù)b1的表達式如式（7）所示.

將現(xiàn)有模擬結果與文獻[14-15]的結果進行比較，結果如圖7所示（本文的Df與文獻[14-15]的Df統(tǒng)一為定值）.由圖7可以看出，文獻[14]提出的表達式往往過高地估計系數(shù)b1值，而文獻[15]給出的關系式往往低估了系數(shù)b1值.

圖7 b1Df與ε的關系曲線Fig.7 Relationship between b1Dfandε

綜上所述，隨著Re的增大，達西定律出現(xiàn)偏差，這主要是由于慣性力的影響，使得隨機纖維多孔介質流動出現(xiàn)不同的流動區(qū)域.本文最后得到的適合非達西現(xiàn)象的壓力損失拓展表達式如下所述.

線性流區(qū)（≤0.3m／s）：

4 結 語

本文通過創(chuàng)建隨機排列的三維纖維過濾介質微觀結構模型，采用計算流體動力學軟件模擬計算不同流體平均速度及多孔介質填充率條件下三維纖維過濾介質內部的氣相流場，得出如下主要結論.

（1）在本文建立的模型中，當Re≤0.33或當≤0.3m／s時，達西定律可以精確地描述纖維過濾多孔介質的流場，即壓力損失與速度呈線性關系；而隨著Re或速度逐漸增大，壓力損失隨著其增大呈現(xiàn)出非線性增加，則達西定律不適用于描述此時的流體流動.

參 考 文 獻

[1]HAO Z Y，CHENG N S，TAN S K.Investigation of inertial effect in simplified porous media flow[J].Advances in Water Resources and Hydraulic Engineering，2009，1：160-165.

[2]VAFAI K.Convection flow and heat transfer in variable porous media[J].Journal of Fluid Mechanics，1984，147：233-259.

[3]FORCHHEIMER P.Wasserbewegung durch Boden[J].Zeitschrift Verein Deutscher Ingenieure，1901，45：1782-1788.

[4]DEIBER J，BORTOLOZZI R.A two-field model for natural convection in a porous annulus at high Rayleigh numbers[J].Chemical Engineering Science，1998，53（8）：1505-1516.

[5]VAN DER SMAN R G M.Prediction of airflow through a vented box by the Darcy-Forchheimer equation[J].Journal of Food Engineering，2002，55（1）：49-57.

[6]JR REIS N C，GRIFFITHS R F，SANTORS J M.Numerical simulation of the impact of liquid droplets on porous surfaces[J].Journal of Computational Physics，2004，198（2）：747-770.

[7]YANG D X，YANG Y P，COSTA V A F.Numerical simulation of non-Darcian flow through a porous medium [J].Particuology，2009，7（3）：193-198.

[8]MEI C C，AURIAULT J L.The effect of weak inertia on flow through a porous medium[J].Journal of Fluid Mechanics，1991，222：647-663.

[9]ANDRADE J S，COSTA U M S，ALMEIDA M P，et al.Inertial effects on fluid flow through disordered porous media[J].Physical Review Letters，1999，82（26）：5249-5252.

[10]CAI Z H，SHI B C，LU J H，et al.Non-Darcy flow in disordered porous media：A lattice Boltzmann study[J].Computers ＆ Fluids，2010，39（10）：2069-2077.

[11]TENG H，ZHAO T S.An extension of Darcy＇s law to non-Stokes flow in porous media[J].Chemical Engineering Science，2000，55（14）：2727-2735.

[12]HOSSEINI S A，VAHEDI TAFRESHI H.3-D simulation of particle filtration in electrospun nanofibrous filters[J].Powder Technology，2010，201（2）：153-160.

[13]李艷艷.纖維過濾材料結構參數(shù)優(yōu)化數(shù)值模擬[D].上海：東華大學環(huán)境科學與工程學院，2011：20-22.

[14]ERGUN S.Fluid flow through packed columns[J].Chemical and Engineering Progress，1952，48（2）：89-94.

[15]LEE S L，YANG J H.Modelling of Darcy-Forchheimer drag for fluid flow across a bank of circular cylinders [J].International Journal of Heat and Mass Transfer，1997，40（13）：3149-3155.