顆粒床內(nèi)固液過濾的三維CFD-DEM模擬

肖 桐， 王千紅， 沈盈鶯， 惠 卿， 陳 晗， 程振民

（華東理工大學化工學院，化學工程聯(lián)合國家重點實驗室，上海 200237）

顆粒床過濾是一種典型的深床過濾，與表面過濾方式不同，它主要將懸浮物攔截在過濾介質(zhì)的內(nèi)部。顆粒床過濾器結(jié)構(gòu)簡單、成本低廉、運行費用少、過濾效果好，在石油化工、生物制藥、污水處理、農(nóng)田灌溉等領(lǐng)域應用廣泛[1-4]。為深入了解深床過濾機理和規(guī)律，更好地指導過濾器的工藝設(shè)計和操作運行，研究者們構(gòu)建了不同的經(jīng)驗和理論模型，包括唯象模型[5]、跡線理論模型[6]、隨機模型[7]、逾滲理論模型[8]和網(wǎng)絡(luò)模型[9]等。但是這些經(jīng)驗或理論模型有許多缺陷：唯象模型是一種經(jīng)驗模型，有些參數(shù)并沒有實際的物理意義；跡線理論模型不能反映過濾過程中的水頭損失影響；隨機模型不能預測濾出液的顆粒濃度；逾滲理論模型不適合非截留機理主導過濾的情況；網(wǎng)絡(luò)模型假設(shè)的通道和實際的過濾床層孔道差別較大。因為這些數(shù)學模型的局限性，所以需要探索更為先進的研究方法。

隨著計算機技術(shù)的發(fā)展和離散單元法的提出，基于計算流體力學（CFD）和離散單元法（DEM）的數(shù)值模擬研究方法被越來越多的學者所采用。Zhao等[10]用耦合CFD-DEM的方法模擬與采礦和巖土工程相關(guān)應用中的流體和顆粒相互作用的行為。方乙琳等[11]運用CFD-DEM耦合的手段，結(jié)合標準k-ε模型（k為湍流動能，ε為湍流耗散率）、Eulerian多相流模型，模擬沸騰床外排廢催化劑的分選，得到催化劑分離狀態(tài)及流場變化，探究分選柱中催化劑的分選機理和顆粒運動特性。謝曉翠等[12]利用CFD軟件中的流體體積函數(shù)(VOF)模型捕捉氣液相界面流動，利用離散相模型（DPM）考察納米顆粒與水相的相互作用，進而模擬納米顆粒、水和空氣的流動。Tao等[13]采用三維多時間尺度CFD-DEM耦合方法研究不同類型纖維陣列對微米級顆粒的過濾，并比較了以平行和交錯方式均勻排列的纖維的過濾性能。Schilling等[14]用Euler-Lagrange方法和CFD-DEM耦合數(shù)值模擬濾料纖維上顆粒的遷移和沉積行為，并分析了過濾纖維上由沉積顆粒構(gòu)成的樹枝狀結(jié)構(gòu)的重新分配。Qian等[15]在研究纖維層過濾氣溶膠時，構(gòu)建出基于隨機算法的三維纖維床層，并將CFDDEM耦合模型應用到纖維床層過濾中，探究不同空隙率的纖維床層中雜質(zhì)顆粒的沉積形態(tài)，以及過濾效率和床層壓降隨表觀過濾速率、床層空隙率、雜質(zhì)顆粒大小的變化。Wang等[16]基于CFD模擬，對兩種典型填料顆粒過濾器的顆粒過濾特性進行了數(shù)值模擬，得到濾料顆粒表面上的流體流動、壓降、顆粒沉積特性以及顆粒沿過濾床的沉積分布。

本文采用CFD-DEM耦合方法，對高20 mm、顆粒直徑為1.0 mm的三維隨機堆積顆粒床進行固液過濾模擬研究。這種雙平臺耦合計算的方法，可以避免刻畫顆粒內(nèi)部的復雜流道，大大減少前期準備的工作量，為顆粒床過濾的模擬研究提供新思路。

1 模型描述和數(shù)值計算方法

1.1 數(shù)值計算模型

1.1.1 流體相的控制方程 在Eulerian模型中，考慮到雜質(zhì)顆粒沉積對濾料床層空隙率的影響，在Navier-Stokes方程中引入床層空隙率 ( φ) ，從而實現(xiàn)流體對顆粒及顆粒對流體的雙向耦合分析。流體相的連續(xù)性方程和動量方程分別表述如下：

其中：μ是流體黏度，Pa·s； u 是流體速度，m/s； ρf是流體密度，kg/m3； g 是重力加速度，m2/s；S是作用在網(wǎng)格單元體積的流體阻力 FD的總和，N； Δ V 是網(wǎng)格單元體積，m3。

對于填料床內(nèi)流體的流動狀態(tài)的判斷，景有海[17]依據(jù)前人的研究，提出Reynolds數(shù)計算公式（式（4）），只有當 Re 小于2時，流體流動狀態(tài)才是層流。

其中：dp為顆粒直徑； v 是濾層中水流平均速度； α 是濾料顆粒的面積形狀系數(shù)，取值為1。據(jù)此計算的層流臨界速率為0.006 8 m/s，而本實驗中的流速均大于0.13 m/s，所以，應該采用湍流模型計算。

應用SST k-ω湍流模型計算顆粒床層內(nèi)部的湍流流動。SST k-ω模型從k-ε和k-ω模型發(fā)展而來，可以在廣泛的領(lǐng)域中獨立于k-ε模型，使得在近壁自由流中k-ω模型有廣泛的應用范圍和精度。模型中的湍流動能k和比耗散率ω根據(jù)式（5）和式（6）計算：

其中： ui是流體速度的分量；Gk表示湍流的動能；Gω是由 ω 方程產(chǎn)生的項；Γk與 Γω分別代表k與 ω 的有效擴散項；Yk與Yω分別代表k與 ω 的發(fā)散項；Dω代表正交發(fā)散項；Sk與 Sω為自定義項。其他詳細的SST k-ω湍流計算模型可以參見文獻[18-19]。

1.1.2 固體相的離散模型 顆粒床過濾中包含兩種離散相：一是濾料床層的顆粒固相，二是流體中雜質(zhì)顆粒相。懸浮液中的雜質(zhì)顆粒在流體中運動，單一顆粒所受外力來自顆粒與顆粒之間作用力以及流體對顆粒的作用力。運用Lagrange方法研究雜質(zhì)顆粒在流體中的運動軌跡。由于模型中的雜質(zhì)顆粒足夠大，所以，本文只考慮流體對顆粒的曳力、顆粒自身的重力以及顆粒與顆粒的碰撞等因素，忽略范德華力、雙電層作用力和布朗運動等因素對顆粒的影響。顆粒在流場中的運動方程為：

其中：mp是顆粒質(zhì)量，kg； up是顆粒速度，m/s； Fg是顆粒的重力，N；fD(u-up) 是流體對顆粒的曳力，N。

fD可通過式（8）求得：

其 中非線性曳力系數(shù) (CD) 為：

顆粒與顆粒之間的接觸作用采用Hertz-Mindlin接觸模型[20]進行計算，這種接觸作用可用非線性的軟球模型表示。顆粒與顆粒之間的作用力和力矩表示如下：

其中： IP是顆粒的轉(zhuǎn)動慣量，kg·m2； wp是顆粒的轉(zhuǎn)動速度，rad/s； Mi是顆粒碰撞扭矩，N·m。詳細的顆粒與顆粒的作用力和力矩計算參見文獻[21]。

1.1.3 耦合計算方法 CFD與DEM的耦合采用Eulerian耦合法，除考慮液相和固相之間的動量交換外，還考慮固體顆粒對流體相的影響。CFD平臺采用Fluent軟件，DEM平臺采用EDEM軟件，兩個平臺之間的耦合程序采用自編UDF（User-Defined Function，用戶自定義函數(shù)）。如圖1所示，首先在Fluent中進行液相流場初始化，在一次時間步長內(nèi)，流場計算至收斂；然后在EDEM中確定顆粒位置。考慮到流體對顆粒的曳力、顆粒自身以及顆粒之間的作用等，計算顆粒碰撞方程，更新顆粒位置信息；將顆粒對流體的作用力、顆粒位置等信息，傳遞給Fluent進行下一時間步長內(nèi)液相流場的計算。

圖 1 CFD-DEM耦合計算流程Fig. 1 Flow chart of the CFD-DEM coupling

1.2 三維過濾模型

圖2所示為隨機堆積的顆粒床層。通過EDEM中圓柱形“顆粒工廠”快速生成隨機位置的顆粒，在重力作用下于過濾床內(nèi)堆積，進而形成隨機堆積的三維過濾模型。圖3所示為數(shù)值計算模型中的顆粒過濾床，直徑（R）為10 mm，高（H）為20 mm，其中球形顆粒的直徑（D）為1.0 mm，床層空隙率（φ）為0.43。額外的入口和出口延伸區(qū)域用以維持顆粒均勻的入口速度并防止其從出口回流。流體介質(zhì)為常溫水，密度約為1 000 kg/m3，黏度約為0.001 Pa·s。具體幾何參數(shù)和表觀過濾速率（uA）等模擬條件見表1，顆粒相關(guān)參數(shù)的設(shè)置見表2。網(wǎng)格單元采用結(jié)構(gòu)性網(wǎng)格劃分，網(wǎng)格數(shù)為2 261，最小體積為6.5×10-10m3，最大體積為3.7×10-9m3。Fluent中的時間步長設(shè)置為2×10-3s，為了準確捕獲顆粒接觸行為，EDEM的時間步長比Fluent的時間步長通常要小很多，這里兩者比值設(shè)置為1∶100。所有的模擬計算均在Intel?Xeon?E5645 2.4 GHz（24核，48 GB內(nèi)存）服務(wù)器上進行，計算的平均時長約為48 h。

圖 2 隨機堆積的顆粒床層Fig. 2 Randomly packed granular layer

2 數(shù)值模型的實驗驗證

在顆粒床過濾實驗中，過濾效率 ( η ) 定義為：

其中： c1為母液質(zhì)量濃度，mg/L； c2為濾出液質(zhì)量濃度，mg/L。假設(shè)每個懸浮顆粒質(zhì)量相等，采取計數(shù)的方法代替質(zhì)量濃度的測量。因此，過濾效率等于母液和濾出液中顆粒數(shù)量之差與母液中顆粒數(shù)量的比值，初始過濾效率（η0）定義為：在零時刻，濾層處于干凈狀態(tài)下的過濾效率。

圖 3 DEM計算區(qū)域和邊界條件Fig. 3 Illustration of the DEM solution domain and boundary conditions

表 1 幾何參數(shù)和模擬條件Table 1 Geometric parameters and the simulated conditions

表 2 模擬計算參數(shù)Table 2 Parameters in the simulated calculation

為了驗證三維模型的CFD-DEM耦合計算方法的可行性，設(shè)計了小型過濾床實驗，如圖4所示。過濾床內(nèi)的濾料是直徑為1.0 mm的氧化鋁顆粒，過濾床高為20 mm，雜質(zhì)粒徑分別選取100 μm和200 μm。不同表觀過濾速率下過濾床對雜質(zhì)顆粒的過濾效率的實驗值與計算值比較如圖5所示。由圖5可以看出模型計算值和實驗測量值基本吻合。

3 結(jié)果與討論

3.1 初始過濾效率和床層壓降

圖6所示為雙因素（雜質(zhì)顆粒大小、表觀過濾速率）下初始過濾效率的預測值和實驗值的比較，以及由計算值擬合得出的多項式曲面過濾效率的預測。由圖6可以看到，過濾床層對較大雜質(zhì)顆粒的攔截效率高，表觀過濾速率的增大會減小過濾效率。由于過濾效率與雜質(zhì)顆粒直徑（dp）、表觀過濾速率有關(guān)，初始過濾效率可以寫成：

圖 4 過濾床示意圖Fig. 4 Schematic diagram of the filter bed

圖 5 過濾效率的實驗值與模擬計算值比較Fig. 5 Comparison of the filtration efficiencies between experimental and simulated data

圖 6 雙因素下的初始過濾效率： (a) 擬合曲面；(b) 預測值與實驗值Fig. 6 Initial filtration efficiency under two factors: (a) Fitting surface; (b) Predicted value and experimental value

利用Matlab的多項式函數(shù)進行非線性擬合，在130 mm/s＜uA＜330 mm/s，100 μm ＜ dp＜ 200 μm范圍內(nèi)得到初始過濾效率的計算式，如式（14）所示。需要注意 的是，式（14）中 uA、 dp的單位分別是mm/s、mm。

繪制出計算值的擬合曲面，并與100 μm和200 μm時的實驗值進行對比，結(jié)果如圖6（a）所示。

由回歸分析得到的擬合曲面和多項式函數(shù)關(guān)系，進行初始過濾效率的預測。如圖6（b）所示，100 μm和200 μm雜質(zhì)顆粒的過濾效率，其預測值和實驗值的偏差均在10%以內(nèi)。Ergun方程[22]可以用來計算填料床層的壓降值，如式（15）所示。

其中： ψ 為濾料顆粒的球形度，取值為1； de為濾料顆粒的體積等效直徑，與 dp相等。

圖7所示為過濾床的初始壓降。由圖7可得干凈床層的初始壓降的模型計算值與Ergun方程所得壓降值的比較。Ergun方程的計算誤差在±25%以內(nèi)有效，雖然模型計算值相比Ergun方程計算值偏小，但誤差都在25%以內(nèi)。

3.2 過濾速率對顆粒沉積的影響

雜質(zhì)顆粒在過濾中的一個重要行為是聚集。當雜質(zhì)顆粒聚集在一起形成不同的簇（聚集體）時，則發(fā)生聚集現(xiàn)象。新形成的簇中的每個單獨的顆粒將失去動力學獨立性。較大的顆粒（聚集體）在過濾期間更可能沉積。當已經(jīng)沉積的顆粒受到流體的沖刷、顆粒的碰撞等外力作用時，有可能再次剝離沉積表面，這些現(xiàn)象顯然與過濾速率有關(guān)。

圖 7 過濾床的初始壓降Fig. 7 Initial pressure drop of the filter bed

圖8示出了過濾一段時間（10 s）時，不同表觀過濾速率下同一縱截面的顆粒沉積情況（藍色代表顆粒速度接近零，紅色代表顆粒速度最大，綠色代表顆粒介于兩者之間，全文同），流體走向為從上至下?？梢钥吹?，當表觀過濾速率為130 mm/s時，懸浮顆粒在濾料間隙成團的比例較大；當表觀過濾速率增大到230 mm/s時，“團簇”的顆粒數(shù)量減少；當表觀過濾速度增大到330 mm/s時，幾乎看不到成團，且顆粒沉積量明顯減少。

圖 8 不同表觀過濾速率下同一縱截面的顆粒沉積Fig. 8 Particle deposition in the same longitudinal section at different superficial velocities

在EDEM的后處理界面，使用“Clipping”功能將過濾床模型等分為20層，并逐個統(tǒng)計每層中雜質(zhì)顆粒數(shù)隨時間的變化關(guān)系，然后單獨統(tǒng)計某一時刻不同層的顆粒數(shù)量。圖9所示為相同時刻（10 s）、不同表觀過濾速率下過濾床不同深度處的雜質(zhì)顆粒數(shù)量。可以看出，隨著表觀過濾速率的增大，每層的雜質(zhì)顆粒數(shù)量有減少的趨勢，總的沉積顆粒數(shù)量相應減少。所以，表觀過濾速率會影響雜質(zhì)顆粒在濾層中的成團趨勢，會使得已經(jīng)沉積的顆粒再次返回流體中，進而影響濾層的過濾效率。

圖 9 不同表觀過濾速率下不同床層深度處的雜質(zhì)顆粒數(shù)量Fig. 9 Numbers of particles at different bed depths and superficial velocities

3.3 過濾時間對顆粒沉積的影響

圖10示出了過濾床中不同深度的雜質(zhì)顆粒隨時間的沉積變化情況。將同一時刻（分別選取1，5、10 s）、不同顆粒層內(nèi)的雜質(zhì)顆粒數(shù)統(tǒng)計到一起，構(gòu)成雜質(zhì)顆粒數(shù)隨過濾床深度的變化圖，如圖11所示。圖10和圖11均在表觀過濾速率為130 mm/s，雜質(zhì)顆粒直徑為100 μm的條件下進行，其他條件參見表1。其中圖10是分別選取距濾層頂部0、10 mm和19 mm截面處的雜質(zhì)顆粒沉積情況。

從圖10和圖11可以看到，沿著過濾床從上到下，雜質(zhì)顆粒數(shù)量呈現(xiàn)減少趨勢；靠近入口處的雜質(zhì)顆粒成團趨勢明顯；在靠近出口處，速度較大的雜質(zhì)顆粒數(shù)量比例較大。再對比不同時間的顆粒沉積情況，取前3層作為表層，后3層作為底層。圖11（a）中柱狀分布較“尖銳”，表層顆粒數(shù)超過45%，而底層顆粒數(shù)僅為4%；圖11（b）、（c）中柱狀分布趨向平緩，表層顆粒數(shù)分別為29%、25%，底層分別為14%、17%，表層和底層顆粒數(shù)的差值進一步減少。從而說明，在過濾前期，靠近入口處的雜質(zhì)顆粒數(shù)量變化較大；當過濾逐漸加深，靠近入口處的雜質(zhì)顆粒數(shù)量變化減弱，較深的床層處也逐漸積累更多的雜質(zhì)顆粒，最后床層積累的雜質(zhì)顆粒數(shù)量越來越均勻。

圖 11 過濾床內(nèi)不同時刻不同深度處的雜質(zhì)顆粒數(shù)占所有沉積顆粒數(shù)的比例Fig. 11 Ratio of impurity particles to all deposited particles at different depths of the filter bed at different time

3.4 顆粒沉積對過濾效率和壓降的影響

隨著過濾時間的進行，雜質(zhì)顆粒的沉積會影響過濾床內(nèi)的流體流動、過濾效率和床層壓降。雜質(zhì)顆粒隨流體進入過濾床層，受到濾料顆粒的攔截作用，一部分雜質(zhì)顆粒沉積在濾料顆粒表面，另一部分受前面顆粒架橋作用的影響而被攔截，但仍有被流體沖刷而脫離的可能性。其他的顆粒隨流體被帶出過濾床層。沉積在床層的顆粒會改變后來流體和顆粒的流道，影響其流動速率和被攔截的概率，進而影響過濾床層的過濾效率和壓降。

圖12（a）、12（b）分別為在表觀過濾速率為230 mm/s、連續(xù)過濾條件下，過濾床的過濾效率和床層壓降隨單位體積沉積量增長的變化情況。可以看到，對于不同大小的雜質(zhì)顆粒，過濾效率均有所提高，但是床層壓降也隨之增大。雜質(zhì)顆粒的沉積會改變?yōu)V料床層的流體通道大小，使得通道尺寸小于雜質(zhì)顆粒尺寸的可能性更大，因此床層攔截的雜質(zhì)顆粒更多，過濾效率更高。小顆粒之間的縫隙比大顆粒之間的縫隙小，流體流過小顆粒之間的縫隙時的速率更大，并且床層空隙率也隨顆粒沉積的增多而變小。床層壓降與流速和空隙率是正相關(guān)關(guān)系，所以小顆粒雜質(zhì)沉積產(chǎn)生的壓降就較高。

圖 12 過濾效率和床層壓降隨雜質(zhì)顆粒單位體積沉積量的變化Fig. 12 Filtration efficiencies and pressure drop variations with the particle deposition mass per unit volume

4 結(jié) 論

采用CFD-DEM耦合方法，對高20 mm、顆粒直徑為1.0 mm的三維隨機堆積顆粒床層進行固液過濾研究。主要結(jié)論如下：

（1）隨著床層的加深，壓降值逐漸增大；流體因為顆粒之間的縫隙減小而加快速度。與Ergun方程相比，不同表觀過濾速率下的三維模型壓降值的計算偏差在25%以內(nèi)。

（2）過濾床層對較大雜質(zhì)顆粒的攔截效率高，表觀過濾速率的增大會減小過濾效率。通過回歸分析得到擬合函數(shù)，并進行初始過濾效率的預測，對于100 μm和200 μm的雜質(zhì)顆粒的過濾效率，其預測值和實驗值的偏差均在10%以內(nèi)。

（3）沿著過濾床從上到下，雜質(zhì)顆粒數(shù)量呈現(xiàn)減少趨勢；靠近入口處的雜質(zhì)顆粒成團趨勢明顯；在靠近出口處，速度較大的雜質(zhì)顆粒數(shù)量較多。在過濾前期，靠近入口處的雜質(zhì)顆粒數(shù)量變化較大；當過濾逐漸加深，靠近入口處的雜質(zhì)顆粒數(shù)量變化減弱，床層積累的雜質(zhì)顆粒數(shù)量越來越均勻。

（4）隨著過濾時間的增加，雜質(zhì)顆粒的沉積會影響過濾床內(nèi)的流體流動，進而影響過濾效率和床層壓降。對于不同大小的雜質(zhì)顆粒，隨著時間進行，其過濾效率均有所提高，但是引起的床層壓降也隨之增大。