顆粒床內固液過濾的三維CFD-DEM模擬

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

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

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

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

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

1 模型描述和數值計算方法

1.1 數值計算模型

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

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

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

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

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

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

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

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

fD可通過式（8）求得：

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

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

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

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

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

1.2 三維過濾模型

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

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

2 數值模型的實驗驗證

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

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

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

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

表 2 模擬計算參數Table 2 Parameters in the simulated calculation

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

3 結果與討論

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

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

圖 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的多項式函數進行非線性擬合，在130 mm/s＜uA＜330 mm/s，100 μm ＜ dp＜ 200 μm范圍內得到初始過濾效率的計算式，如式（14）所示。需要注意 的是，式（14）中 uA、 dp的單位分別是mm/s、mm。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結 論

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

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

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

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

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