基于CFD-DEM 的固液分級過濾模擬

趙鐘杰， 張建鵬， 唐艷玲， 肖 桐， 黃子賓， 程振民

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

顆粒層過濾作為一種低成本的分離方式被人們廣泛關注，并應用于煙氣除塵、水處理、油漿凈化等領域[1-3]。研究者們總結了包括慣性碰撞、攔截和擴散等在內的基本顆粒層過濾機理，得出了過濾性能與流速、粒徑、床層深度等因素有關的結論[4-6]。Zamani 等[7]歸納了顆粒層過濾的微觀和宏觀模型，這些模型的應用范圍是非普適性的，均存在各自的缺點：如宏觀經驗模型無法揭示內部過濾機理；隨機模型無法計算過濾效率；跡線模型無法計算過濾壓降等。

由于上述模型應用的局限性，計算流體力學(Computational Fluid Dynamics，CFD) 和離散單元法(Discrete Element Method，DEM) 耦合模擬的方法被許多研究者所重視，并被應用于過濾研究。Qian 等[8]利用隨機算法構建了三維纖維過濾器模型，并采用CFD-DEM 耦合方法研究不同氣速和空隙率下的氣溶膠過濾性能。Yue 等[9]基于CFD-DEM 的耦合模型對含塵氣體的過濾進行了模擬，全面分析了顆粒的運動以及顆粒所受接觸力對過濾性能的影響。肖桐等[10]采用CFD-DEM 耦合方法，構建三維隨機堆積的顆粒床模型，研究不同時刻和表觀過濾速率下顆粒沉積的變化，并考慮了顆粒沉積對于過濾效率和壓降的影響。

顆粒層過濾器中雜質顆粒的沉積分布通常是不均勻的，當濾料粒徑較小時，顆粒更容易沉積在過濾器前端[11]，床層深處的濾料幾乎沒有被利用。濾料粒徑越小，過濾效率越高，然而這也意味著過濾器容垢能力的降低。為了同時提高過濾效率和容垢量，分級過濾器的概念被提出。Yang 等[12]進行了分級過濾器的除塵實驗，上、下層分別填裝輕而大和重而小的濾料。實驗結果表明，在相同的床層壓降下，分級過濾器的容塵能力是單層過濾器的10 倍。Shi 等[13]通過實驗研究了上層膨脹珍珠巖和下層細砂組成的分級過濾器的過濾性能，結果發現分級過濾器的壓降增長率與單層粗濾料過濾器基本相同，且遠低于單層細濾料過濾器的壓降增長率。

綜上所述，可發現前人對于分級過濾的研究主要集中在實驗層面，而雜質顆粒在過濾器內的沉積分布無法通過實驗方法準確獲得。本文利用CFDDEM 耦合方法，構建不同濾層結構的三維過濾器模型，獲取顆粒沉積分布結果以輔助解釋實驗結果。

1 實驗描述

圖1 為顆粒層過濾實驗裝置簡圖，包括原料液供應和控制系統、輸送原料液的泵系統、顆粒過濾系統。充分攪拌的原料液通過泵系統輸送至濾層總高為40 mm 的過濾器的入口，隨后原料液流經床層，未被收集的顆粒從過濾器出口逃逸。

實驗中使用的濾料為α-Al2O3，并采用分級法進行填裝。上層填充的粗濾料(記為D1)粒徑為4 mm，下層填充的細濾料(記為D2)粒徑為2 mm。圖2 所示為5 種不同濾層結構的顆粒層過濾器，其中Ⅰ、Ⅴ為單層過濾器，Ⅱ、Ⅲ、Ⅳ為分級過濾器。利用流量計測量表觀過濾速率(u)，并分別控制表觀過濾速率為50、100、150、200、250 mm/s。原料罐內料液的密度近似為1 000 kg/m3，動力黏度為0.001 Pa·s。常規的固液顆粒層過濾器中顆粒的體積分數在0.01%以下[14]，本文中顆粒體積分數為0.005 6%，質量濃度控制在180 mg/L。進入過濾器的雜質顆粒通過細篩分獲得，直徑約為0.4 mm。每次實驗在20 ℃的室溫下重復3 次，以盡量消除偶然性。

圖2 不同濾料填裝方式的過濾器Fig.2 Filters with different granular packing methods

2 數值方法和模型描述

2.1 數值方法

2.1.1 流體相控制方程 選取Eulerian 模型描述流體相的行為變化，并在控制方程中引入了過濾器的床層空隙率ε，實現流體相和離散相的雙向耦合，以考察顆粒沉積對于過濾的影響。流體相的控制方程如下：

其中：ρf是流體密度；g是重力加速度；μ是流體黏度；S是作用在網格單元體積內的阻力FD的總和，ΔV表示網格單元體積。

流體在顆粒床層內的流動狀態主要取決于床層雷諾數Re[15]：

其中：ρ為顆粒密度，D是顆粒直徑，v是水流平均速率，α為濾料形狀系數。當Re>2 時，流體流動狀態為湍流。

根據式(4)，若視顆粒為球體，則層流臨界速率(Re=2 時的水流平均速率)為0.003 43 m/s。由于在文中最小表觀過濾速率為50 mm/s，因此選擇湍流模型。采用SSTk-ω模型描述介質空隙間湍流的影響，湍流動能k和比耗散率ω的方程分別為[16-17]：

其中：xi、xj分別為i、j方向上的坐標值；Γk和Γω為湍流模型的有效擴散項；Gk為湍流動能；Gω代表ω方程；Yk和Yω表示湍流模型的發散項；Dω表示正交散度項；Sk和Sω由用戶自定義。

2.1.2 固體離散相模型 在過濾過程中，由于顆粒與顆粒、濾料和流體相間作用力的存在，使得顆粒發生遷移，遷移軌跡通過Lagrange 方法求解。通常微米級顆粒才需要考慮非接觸力，而本文中模擬涉及的顆粒粒徑遠大于微米級，所以不考慮非接觸力的影響。顆粒間的接觸力則利用Hertz-Mindlin 模型[18]計算，并表示為非線性的軟球模型。在流體相中，顆粒的遷移運動方程如式(7)、式(8)所示：

其中：mp、up、ωp、fn、ft和Tp分別表示顆粒的質量(g)、速度(m/s)、角速度(rad/s)、法向碰撞力(N)、切向碰撞力(N)、碰撞力矩(N·m)；Ip是顆粒的慣性矩(m4)；Fg和fD(u-up)分別是顆粒的重力(N)和流體對于顆粒的阻力(N)，fD表示為：

其中：ρp是顆粒密度；CD是非線性阻力系數，表達式如下：

有關顆粒之間力的詳細計算，可參閱文獻[19]。

2.1.3 耦合計算方法 在CFD-DEM 雙平臺數值模擬過程中采用Eulerian 耦合法，其不僅考慮液固兩相間的動量交互過程，同時也考慮顆粒對流場的影響[10]。DEM 平臺選用EDEM (V 2018)，CFD 平臺選用Fluent (V 16.0)，耦合流程由編譯的UDF (User-Defined Function，用戶定義函數) 加載。

圖3 為CFD-DEM 耦合流程圖，首先在Fluent 中初始化流場并進行零時刻流場數據的計算，隨后將流場數據傳輸至EDEM 中，并利用Hertz-Mindlin 模型計算顆粒受到的曳力等所有接觸力，由此計算顆粒速度等信息；將作用力等信息通過耦合程序傳輸至Fluent 中進行下一時間步長流場數據的計算，直至達到目標模擬時間。計算流場的時間步長取為3×10-4s，為準確獲取顆粒接觸力的信息，EDEM 時間步長設置為3×10-6s。本文的計算在Intel(R) Core(TM) i7-9750H CPU @ 2.60 GHz (6 核心，16 GB 內存)上進行，計算平均時長為60 h。

圖3 CFD-DEM 耦合流程圖Fig.3 Flow chart of the CFD-DEM coupling

2.2 過濾模型

圖4 所示為過濾器的幾何模型(以過濾器Ⅳ為例)及其邊界條件，模型整體是半徑為25 mm 的三維圓柱，DEM 計算域的高度為40 mm。濾料由EDEM中的“顆粒工廠”在DEM 計算域內隨機生成并快速堆積。計算域底部的交織濾網在支撐濾料的同時保證雜質顆粒都能通過。流體的入口和出口分別采用速度進口和壓力出口，進口速率與實驗保持一致。為了保持進口速度的穩定，并避免流體從出口回流，分別對模型的進口和出口區域進行了加長處理。顆粒由顆粒入射面進入，并從壓力出口邊界離開模擬區域。模型的壁面設置為靜止且無滑移，模擬參數如濾料粒徑和雜質顆粒粒徑等均與實驗條件相同。表1 所示為模擬所用物性參數。結構化網格被應用于劃分Fluent 中過濾器模型的網格，網格數量為2 261。為了保證數值結果的準確性，EDEM 中過濾器模型的網格數量分別依次劃分為982 600、1 706 256、3 342 336 和7 860 800。模擬結果顯示，網格數量為3342 336 和7 860 800 的模型所得到的初始壓降值和過濾效率均非常接近。因此，從數值結果的精度和計算機資源兩方面考慮，EDEM 中網格數量劃分選取3 342 336。

表1 模擬物性參數Table1 Physical parameters in the simulation

圖4 過濾器模型及其邊界條件Fig.4 Filter model and its boundary conditions

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

其中：Cin為進口顆粒質量濃度，mg/L；Cout為出口顆粒質量濃度，mg/L。假設所有顆粒的物性參數相同并視為球體，則可采用顆粒數目代替顆粒質量濃度的方法計算過濾效率。而初始過濾效率定義為濾料處于干凈狀態下的過濾效率。

Ergun 方程[20]可用于預測顆粒床的初始壓降（ Δp）：

其中： 為顆粒球形度；de為顆粒體積等效直徑，m。

雖然不存在計算容垢量的公式，但很顯然，顆粒沉積均勻度(λ)越大，則相同質量的顆粒沉積引起的壓降增量越小，過濾床的容垢能力就越高。因此，本文采用沉積均勻度來評價過濾器的容垢能力。

圖5 所示為顆粒層過濾器分節示意圖。將過濾器沿著過濾方向分為20 節，顆粒沉積均勻度由式(14)定義：

圖5 顆粒層過濾器分節示意圖Fig.5 Schematic diagram of granular bed filter segmentation

其中：n為床層節數，Ni為每一節床層中沉積的顆粒數目，Nave為沉積顆粒數目的算術平均數。顆粒沉積越均勻，λ越趨向于1。每一節床層中的顆粒沉積分數( ? )定義為：

3 結果與討論

3.1 初始壓降與過濾效率

圖6(a)所示為過濾器Ⅰ、Ⅳ和Ⅴ的初始壓降結果。從圖6(a)可知，壓降模擬結果與實驗結果基本一致，但后者數值略大，這可能是因為在實驗過程中，流體流過過濾器壁面時，由于摩擦而產生了微小的壓降。在相同過濾速率下，過濾器Ⅳ的壓降僅為過濾器Ⅰ和Ⅴ壓降之和的一半左右。

圖6(b) 所示為表觀過濾速率為50 mm/s 時，不同床層深度下過濾器V 的初始壓降模擬結果和Ergun 方程結果之間的比較。Ergun 方程允許的誤差范圍為±25%，雖然模擬計算得到的壓降值稍大，但兩者之間的偏差仍處于有效范圍內。

圖6 (a)不同過濾器初始壓降的實驗值和計算值比較；(b)不同床層深度下過濾器Ⅴ的初始壓降Fig.6 (a) Comparison of initial pressure drop of different filters between experiment and simulation; (b) Initial pressure drop under different bed depths of filter with case Ⅴ

圖7(a) 所示為不同表觀過濾速率下過濾器Ⅰ、Ⅱ和Ⅴ的初始過濾效率變化，黑點和紅點分別代表實驗和計算結果。由圖可得，模擬結果與實驗結果具有良好的一致性。過濾器Ⅰ的過濾效率明顯小于其他濾層結構的過濾器。如當表觀過濾速率為50 mm/s時，過濾器Ⅰ的過濾效率約為70%，而其他濾層結構過濾器的過濾效率均約為90%，說明采用分級填料可以大幅提高過濾效率。

很顯然，通過模擬分析可便捷地得到過濾效率隨床層深度的連續變化。圖7(b)所示為表觀過濾速率為50 mm/s 時，不同床層結構和不同床層深度下的過濾效率模擬值?？梢钥闯?，當過濾器Ⅴ的床層深度約為24 mm(稱為“轉折點”)時，其過濾效率達到最大值并保持穩定。過濾器Ⅲ和Ⅳ的轉折點分別約為32 mm 和28 mm；而過濾器Ⅰ和Ⅱ的轉折點大于40 mm。綜上分析可知，下層細濾料層深度的增加使轉折點提前出現。

圖7 (a)不同過濾器過濾效率的實驗值和模擬值比較；(b)不同過濾器不同床層深度下的過濾效率模擬值Fig.7 (a) Comparison of filtration efficiency of different filters between experiment and simulation; (b) Simulation results of filtration efficiency under different bed depths of different filters

過濾器的綜合性能由過濾效率和壓降共同決定，不同濾層結構過濾器的綜合評價指標(Y) 根據式(16)定義[21]，計算得到的Y值如圖8 所示，可見在不同表觀過濾速率下，分級過濾器相比于單層過濾器具有更高的Y值，即分級過濾的綜合過濾性能最佳。

圖8 不同床層結構過濾器的Y 值Fig.8 Y Value of filters with different bed structures

3.2 顆粒沉積形貌

雜質顆粒沉積形貌的結果更加直觀，可加深對過濾現象的理解。圖9 所示為當表觀過濾速率為200 mm/s、過濾時間為10 s 時，過濾器Ⅱ和Ⅴ的細濾料層表面的顆粒沉積形貌，可發現兩者的顆粒沉積圖是近似的，過濾器Ⅴ的細濾料層表面沉積的雜質顆粒數目大于分級過濾器Ⅱ的細濾料層表面雜質顆粒數。由此可以推斷，分級過濾可以使雜質顆粒部分容納在粗濾料層中，提高了床層容垢量。

圖9 細濾料層表面的顆粒沉積形貌Fig.9 Deposition morphology on the surface of fine granular layer

圖10 所示為當表觀過濾速率為50 mm/s，過濾時間為10 s 時，模擬所得同一縱截面下過濾器Ⅰ、Ⅱ和Ⅴ的顆粒沉積圖?？梢钥闯?，過濾器Ⅴ近入口處及過濾器Ⅱ上、下層濾料的分界處形成明顯的團聚沉積結構。這可以從慣性機制[6]的角度進行解釋，當濾料粒徑為2 mm、過濾速率為50 mm/s 時，斯托克斯數(St)約為1.5，顆粒的慣性發揮作用，使得顆粒容易偏離流體流向，其運動軌跡與濾料表面相交的幾率增大，顆粒發生高頻率碰撞并被快速捕集。較高的捕集幾率促成顆粒間的架橋現象，從而形成團聚結構，其可視為新的“過濾器”來截留顆粒，因而大幅降低了濾層表面的空隙率，也解釋了圖6(a)中過濾器Ⅴ的壓降遠高于其他過濾器的原因。

由圖10 可知，過濾器Ⅰ及過濾器Ⅱ的上層濾料表面并未出現明顯的團聚結構。相比于過濾器Ⅴ，過濾器Ⅰ較大的濾料粒徑使得慣性作用明顯減弱，顆粒更容易隨著流體流線運動而被夾帶至床層的深處，使得顆粒的分布較為均勻。此時可考慮攔截機制[22]，而只有當顆粒的軌跡恰好落在濾料表面到濾料半徑的流線范圍內時，攔截機制才對顆粒遷移產生影響，顆粒才會被濾料表面截留，因此難以形成團聚結構。

圖10 不同過濾器同一縱截面的顆粒沉積形貌Fig.10 Particle deposition morphology in the same longitudinal section of different filters

3.3 顆粒沉積均勻度

從雜質顆粒沉積形貌來分析顆粒沉積分布僅是定性分析，通過EDEM 的“group bin”后處理功能可得到不同床層節數上濾料表面及間隙所沉積顆粒的數量，從而定量分析顆粒沉積均勻度。當表觀過濾速率為200 mm/s、過濾時間為10 s 時，不同床層節數上的雜質顆粒沉積分數如圖11(a)所示。對于過濾器Ⅴ，雜質顆粒的沉積分數隨著床層節數的增大而明顯減小?？梢钥闯?，前5 節床層的顆粒沉積分數之和達到了71.2%；當床層節數達到10 左右時，顆粒沉積分數就幾乎降低為0，后半段床層的顆粒沉積分數之和僅為8.4%，這意味著過濾器Ⅴ收集的顆粒幾乎都容納在過濾器的前端部分。過濾器Ⅰ中每一節床層的顆粒沉積分數基本都為0.05～1.00。對于過濾器Ⅱ，當床層節數小于16 時，每一節床層的顆粒沉積分數相差不大，提供較高容垢量；而當在床層節數為17 時，顆粒沉積分數顯著增大，表明顆粒收集量快速增加，保證了分級過濾器的過濾效率。

圖11 顆粒的沉積分布Fig.11 Particle deposition distribution

圖11(b)定量分析了不同濾層結構下顆粒沉積均勻度的變化。結果表明，沉積均勻度隨表觀過濾速率的增加而增加。如當過濾速率從50 mm/s 增加至250 mm/s 時，過濾器II 的沉積均勻度由0.73 增加到0.84。這是由于過濾速率的增大使濾料表面已沉積顆粒的剝離行為更顯著。過濾器V 的沉積均勻度明顯低于其他4 種濾層結構過濾器。當表觀過濾速率為200 mm/s時，過濾器II 的沉積均勻度比過濾器V 高59.4%；當濾層結構由V 變化至IV，沉積均勻度增加0.217，而從IV 變化到I，沉積均勻度僅增加0.103。

引入分級過濾器的平均濾料直徑(Dave)，并定義為：

式中：D1、D2分別為上、下層濾料直徑；L1、L2分別為上、下層濾料層厚度。

對影響沉積均勻度的各個參數進行量綱分析，通過白金漢定理和量綱一致性，沉積均勻度可表示為：

利用1st Opt 軟件進行多元非線性回歸分析，得到在50 mm/su＜ 250 mm/s, 2 mmDave＜4.0 mm 范圍內，沉積均勻度的關系式如下：

圖12 對比了沉積均勻度的預測結果（ λpre:）與計算結果（ λcal:）。由圖可知，沉積均勻度的模擬結果與回歸方程值吻合良好，偏差小于5%，說明擬合具有足夠的精度。

圖12 沉積均勻度的預測結果和計算結果的對比Fig.12 Comparison of deposition uniformity between predicted and calculated results

4 結 論

采用CFD-DEM 耦合方法，對不同濾層結構的顆粒層過濾器進行固液分級過濾的數值模擬研究，主要結論如下：

(1) 過濾效率的模擬計算值與實驗值吻合良好，壓降值的偏差在Ergun 方程允許誤差范圍內。過濾效率和壓降隨床層深度的增加而增加，當床層深度增大至轉折點時，過濾效率基本不變，且隨著下層細濾料層深度的增加，轉折點提前出現。

(2) 過濾器的容垢能力由顆粒沉積均勻度表示。顆粒沉積均勻度隨表觀過濾速率的增大而增大，濾層的分級填裝使得雜質顆粒在床層中的沉積均勻度大幅提升。通過擬合回歸得到不同濾層結構過濾器沉積均勻度的關聯式，且偏差小于5%。

(3) 雜質顆粒的沉積分布表明：單層細濾料過濾器由于較強的顆粒慣性作用而形成顆粒團聚結構，顆粒沉積主要發生在近入口處，容垢量較低；分級過濾器的上層粗濾料截留了部分的雜質顆粒，在保證過濾效率的前提下，降低了過濾器的壓降，同時提高了容垢量。