半圓形粗糙元壁面顆粒沉積數值模擬研究

黃秋焰，羅坤，樊建人

（浙江大學能源工程學院，浙江 杭州 310027）

隨著工業(yè)化的進步和人們生活水平的提高，粉煤燃燒發(fā)電工廠、汽車廢氣和建筑煙氣等帶來細微顆粒的排放，使大氣中的塵埃顆粒急劇增加，嚴重惡化的城市空氣質量危害人類健康。細顆粒不僅對人體呼吸系統(tǒng)、心血管系統(tǒng)及中樞神經系統(tǒng)造成傷害［1-3］，而且易在氣力輸送、熱交換器、粉末涂料設備等通風管道中沉積與懸?。?-5］，影響設備的正常工作與運轉。大多數通風設備的壁面都不光滑，存在大小不一、形式不同的粗糙元壁面。對粗糙元壁面通風管道顆粒沉積的研究可以為通風管道的設計提供指導。

與光滑元壁面相比，粗糙元壁面對顆粒沉積特性有加強作用。不同粗糙元壁面和粗糙元參數對顆粒沉積也有較大影響。國內外已經有大量相關研究，LAI等［6-7］通過實驗研究了方肋形粗糙元壁面的顆粒沉積特性，結果表明相比于光滑元壁面，方肋形粗糙元壁面的顆粒沉積率有較大提高。LU等［8］通過數值模擬研究了多組方肋粗糙元間距高度比與粗糙元高度管徑比參數下的顆粒沉積加強特性，發(fā)現(xiàn)當粗糙元間距高度比為2，高度管徑比為0.1時，顆粒沉積加強比可達100～1 000。HONG等［9］基于CFD-DEM方法分別對方肋形、尖角形、三角形粗糙元進行了三維數值模擬，結果發(fā)現(xiàn)，尖角形粗糙元的顆粒沉積率最大，而三角形顆粒沉積率最小。需要注意的是，尚未有學者對半圓形粗糙元壁面顆粒沉積進行研究?；诖?，本文研究不同半圓形粗糙元壁面參數（粗糙元間距高度比p/e、高度管徑比e/D）下顆粒沉積特性，并進一步與方形粗糙元壁面顆粒沉積進行了對比。

在顆粒沉積研究模型方面，眾多粗糙元壁面的顆粒沉積研究［8-12］忽略了顆粒在壁面的反彈，將壁面假設為附著壁面，即當顆粒與壁面接觸時認定顆粒沉積，這與實際情況不符。但也有學者在研究顆粒沉積時，考慮了顆粒的反彈、沉積、二次揚塵。如BRACH等［13］根據Johnson-Kendall-Roberts（JKR）理論，通過實驗總結了將半經驗反彈速度作為判斷顆粒沉積的條件，較全面地考慮了顆粒反彈及沉積現(xiàn)象。SOLTANI等［14］采用臨界動量理論判別了顆粒的二次揚塵情形，此方法被很多研究采納［15-17］。

在顆粒沉積研究方法方面，數值模擬已經成為一種重要的研究手段。相比于雷諾時均（RANS）方法，直接數值模擬（DNS）與大渦模擬（LES）在處理工程問題時往往計算成本偏高。因此，RANS方法是工程領域中一種重要的數值計算方法。TIAN等［10］研究了不同湍流模型對顆粒沉積特性的影響，結果發(fā)現(xiàn)，雷諾應力模型（RSM）可以很好地預測顆粒沉積，這是由于RSM考慮了湍流的各向異性，可獲得更準確的流場信息。RSM現(xiàn)已被廣泛運用于顆粒沉積的數值模擬研究［8-9，18-19］。因此，本文采用RSM求解流場。另外，考慮顆粒濃度偏低，屬于稀疏流，采用拉格朗日法對顆粒進行追蹤。

1 物理模型

考慮工業(yè)設備的通風管道大多為矩形管，因此，選取二維矩形通風管道作為模擬對象，管道長400 mm，高度20 mm，如圖1所示。管道下壁面分別用方形和半圓形粗糙元以不同e/D和p/e進行排列，如表1所示。管道上壁面為光滑壁面?？諝鈴墓艿雷髠冗M入，于右側流出。

表1 計算算例Table 1 Computation case

圖1 粗糙元壁面計算域示意Fig.1 The schematic of ribbed ducts

2 研究模型與數值方法

2.1 控制方程

假設空氣為常物性、不可壓縮，RSM控制方程為

其中，和分別為空氣平均速度和平均壓強。雷諾應力項′可用以下方程表述：

其中，νt為空氣湍流黏度，σk，C1，C2為經驗常數，取值分別為1.0，1.8，0.6。湍流耗散率ε的輸運方程為

其 中，σε，Cε1，Cε2為 經 驗 常 數，取 值 分 別 為1.30，1.44，1.92。

由于顆粒濃度偏低，可將顆粒視為離散相，采用拉格朗日法求解離散相的運動軌跡。假設顆粒為球形顆粒，不考慮顆粒的旋轉、顆粒間作用及顆粒對流場的耦合作用，僅考慮流體對顆粒的單向耦合。同時，主要考慮曳力、重力與浮力、薩夫曼升力對顆粒的影響，結合離散相模型（DPM），求解顆粒運動特性。離散相模型控制方程為

其中，Kc為系數，值為2.594，CD為曳力系數：

式中，Rep為顆粒雷諾數，Rep=為 空氣運動黏度。

2.2 顆粒沉積模型

BRACH等［13］在JKR理論的基礎上，通過實驗總結了半經驗反彈速度，并將其稱為顆粒沉積臨界速度Vcr，用以判斷顆粒是否沉積。當顆粒撞擊壁面的法向速度Vn小于臨界速度Vcr時，認為顆粒沉積；反之，認為顆粒與壁面作用后反彈，以一定速度匯入空氣中繼續(xù)運動。臨界速度Vcr表達式為

其中，有效剛度參數

Es和Ep分別是壁面與顆粒的楊氏模量，υs和υp分別為壁面與顆粒的泊松比，dp和ρp分別為顆粒的粒徑與密度，R為顆粒運動恢復系數，取R=0.9。

2.3 算例設置

以5 m·s-1的速度流入，自由流出；空氣在壁面運動采用無滑移邊界條件?？諝膺\動黏度ν=1.48×10-5m2·s-1，密 度ρ=1.225 kg·m-3。顆 粒 密 度 為ρp=2 450 kg·m-3，楊氏模量為192 GPa，泊松比為0.13。壁面材料用不銹鋼，其楊氏模量為215 GPa，泊松比為0.28。為了消除入口效應對顆粒運動的影響，保證湍流充分發(fā)展，在H=0.2 m處釋放45 000個顆粒進行追蹤。通過ANSYS Fluent 18.0自定義函數（UDF）實現(xiàn)顆粒沉積模型。

3 結果與討論

3.1 模型驗證

為驗證顆粒沉積模型在粗糙元壁面上運動的正確性，進行了流場與顆粒沉積特性兩方面驗證。選取X=150 mm的截面進行流場驗證，流場結果與KIM等［20］的DNS流 場數 據 及LU等［8］的RSM流 場數據吻合較好（圖2），獲取了準確的流場信息，包括流場無量綱速度和高度：

圖2 X=150 mm處流場驗證Fig.2 Air flow verification when X=150 mm

其中，U+為流場無量綱速度，y+為流場無量綱高度，u為當地流場速度，u*為摩擦速度，y為距離壁面的高度。

此外，分析了當e/D=0.1，p/e=8時，不同無量綱松弛時間τ+下的顆粒沉積速度Vd，并與LU等［8］的模擬結果進行對比，結果如圖3所示。由圖3可知，顆粒沉積速度整體趨勢一致，但當顆粒粒徑為1 μm時誤差較大，這由未對近壁面處的湍流脈動進行修正造成。當顆粒粒徑為2～10 μm時，本文結果與LU等［8］的結果較吻合，因此，可以用此積塵模型研究半圓形粗糙元壁面的顆粒沉積特性。τ+和Vd的表達式分別為

圖3 顆粒沉積模型驗證Fig.3 Particle deposition model verification

其中，Cc為科靈漢姆修正系數，S為顆粒與空氣密度之比，dp為顆粒粒徑，Ndep為顆粒沉積數量，Nin為進入計算域的顆粒數量。

3.2 空氣流場分析

本研究中的顆粒粒徑為1～10 μm，且僅考慮流場對顆粒的單向耦合作用。流場分析對于顆粒沉積研究起至關重要的作用。圖4為方形和半圓形2種粗糙元在e/D=0.1，p/e=8時的流場分析結果。X和Y分別為二維計算域的橫、縱坐標。由于粗糙元壁面流道變窄，空氣流場加速，且空氣在經過第一個粗糙元后，方形粗糙元壁面比半圓形粗糙元壁面加速明顯。與光滑的上壁面相比，2種粗糙元導致邊界層變厚，且下壁面的邊界層呈波浪式。

圖4 2種粗糙元壁面速度分布Fig.4 Velocity distribution in two ribbed ducts

圖5展示了2種粗糙元附近流場矢量圖與流線圖。當流體流經方形粗糙元壁面時，流場速度立刻增大，而當流體經過半圓形粗糙元壁面時，流場速度在經過半圓形粗糙元壁面頂點后才開始增大。這是由于空氣在經過半圓形粗糙元壁面時緊貼半圓表面流動，且流體流經的橫截面面積未發(fā)生突變。

圖5 2種粗糙元附近流場矢量圖與流線圖Fig.5 Vector and streamline diagram near two rough elements

空氣在繞過粗糙元后，2種粗糙元的背風面均形成一大一小2個回流區(qū)，同時粗糙元的迎風面前端也形成一個回流區(qū)。背風面幾何形狀使得半圓形粗糙元回流區(qū)2的區(qū)域面積略小于方形粗糙元。同樣，空氣的附壁流動導致半圓形粗糙元回流區(qū)1的區(qū)域面積小于方形粗糙元。2種粗糙元回流區(qū)3的區(qū)域面積大致相同?；亓鲄^(qū)不同不僅使顆粒停留時間出現(xiàn)差異，也導致湍流耗散、湍動能等產生差異。

圖6為7種算例下2種粗糙元壁面附近湍動能的分布。不同粗糙元參數下粗糙元壁面的湍動能差異較大。相比半圓形粗糙元壁面，方形粗糙元壁面的突變使得空氣擾動增強（圖5），因此同參數下的方形粗糙元附近空氣的湍動能更大。不同粗糙元高度會引起不同的空氣擾動，從而引起湍動能差異。在算例6中，當粗糙元的e/D=0.2時，其在所有算例中粗糙元高度最大，湍動能最大，壁面附近的最大湍動能達2.8 m2·s-2，而當e/D=0.02時，由于粗糙元高度較低，對流場擾動不明顯，未產生很大湍流效應，其附近湍動能最大僅為0.9 m2·s-2。另外，湍動能對顆粒沉積特性有促進作用，因此需要探究粗糙元參數對顆粒沉積特性的影響。

圖6 2種粗糙元附近湍動能分布（上方為方形粗糙元，下方為半圓形粗糙元）Fig.6 The distribution of turbulent kinetic energy near two rough elements（the upper is square，the lower is semicircle）

3.3 半圓形粗糙元與方形粗糙元顆粒沉積特性對比

圖7展示的為半圓形粗糙元與方形粗糙元在不同顆粒無量綱松弛時間τ+下顆粒沉積速度的對比?？芍?種粗糙元在不同高度與間距下顆粒沉積速度呈現(xiàn)的規(guī)律一致，在2種粗糙元壁面上，顆粒沉積速度總體上均隨顆粒無量綱松弛時間的增大而增大，而當顆粒無量綱時間大于10時，顆粒沉積速度略有下降。這是由于在顆粒無量綱松弛時間為0.1＜τ+＜10時，顆粒受湍流渦作用影響較大［21］，大量顆粒被夾帶進入湍流渦，又由于顆粒的慣性很小，難以逃離湍流渦的束縛，使得顆粒與壁面的碰撞大大增加，導致顆粒沉積現(xiàn)象嚴重。當τ+＞10時，顆粒屬于大粒徑［22］，慣性力占主導作用，顆粒沉積速度Vd保持穩(wěn)定。由于本研究中顆粒范圍較小，τ+＞10的顆粒較少，此處略有下降可以認為顆粒沉積速度處于穩(wěn)定狀態(tài)。

在同一粗糙元參數下，半圓形粗糙元顆粒沉積速度Vd均低于方形粗糙元。由于半圓形粗糙元的回流區(qū)較小，顆粒進入回流區(qū)的概率較低。此外，半圓形粗糙元的外形結構與來流的垂直性不及方形粗糙元，迎風面不能攔截更多的顆粒，所以相比方形粗糙元壁面，半圓形粗糙元壁面沒有出現(xiàn)很嚴重的顆粒沉積現(xiàn)象。

3.4 半圓形粗糙元間距與高度對顆粒沉積特性的影響

圖7（a）展示了當p/e=8時e/D對顆粒沉積特性的影響。當e/D=0.02時，顆粒沉積速度Vd最小且變化較明顯。這是因為粗糙元的高度較小，攔截顆粒運動的效率較低，迎風面與壁面上的顆粒沉積數量都較小，且小尺寸粗糙元對氣流的擾動較小，其附近不會引起較大的湍動能，導致受湍流影響較大的小顆粒跟隨空氣流出計算域，未在粗糙元壁面沉積。而當e/D=0.05，0.10，0.20時，顆粒沉積速度為0.04～0.10 m·s-1，變化不大。

圖7（b）展示了當e/D=0.05時不同p/e下顆粒沉積速度的變化。當p/e=3時，受湍流影響較大的小粒徑在單板壁面沉積較少，而大顆粒受慣性力沉積機制的影響，與其他粗糙元參數的沉積量級一致，使顆粒沉積速度范圍較大。在算例4中，p/e=3，與擁有3個回流區(qū)的其他算例不同的是，雖然算例4的粗糙元引起了較大的湍動能，但算例4中相鄰2個粗糙元的間隙較小，僅形成了一個大的回流區(qū)，顆粒難以向壁面方向運動。另外，本文其他半圓形粗糙元參數下顆粒沉積速度變化不大，這與LU等［8］研究同粒徑方形粗糙元顆粒沉積速度的沉積特性一致。

圖7 2種粗糙元在不同顆粒無量綱松弛時間下的沉積速度對比Fig.7 Velocity distribution with different particle dimensionless relaxation time

3.5 半圓形粗糙元迎風面上的顆粒沉積分布特性分析

圖8展示了半圓形粗糙元在5種不同粒徑下的顆粒沉積分布特性，橫坐標X表示計算域的位置，縱坐標表示顆粒沉積數量。由圖8可知，將顆粒從X=0.20 m處釋放，其在第一個粗糙元的迎風面大量沉積，出現(xiàn)顆粒沉積的峰值，且在下游的每個粗糙元迎風面上均出現(xiàn)峰值。這主要由粗糙元凸起的攔截作用造成，使大部分顆粒在粗糙元的迎風面沉積。而算例4與算例7顆粒沉積峰值不明顯的主要原因是粗糙元的e/D與p/e較小。算例7中，較小的粗糙元不僅顆粒攔截效率較低，而且對流場擾動較小，粗糙元附近湍動能不大，從而導致分峰不明顯。在算例4中，p/e=3，相鄰2個粗糙元之間的空隙僅為e。因此，距離較近導致分峰不明顯。

圖8 在不同粗糙元參數下半圓形粗糙元壁面的顆粒沉積分布Fig.8 Distribution of particles deposition on the wall with different roughness parameters

另外，e/D很好地描述了壁面的相對粗糙度。由圖7可知，粗糙元引起的流場擾動會影響湍流性質，從而對顆粒沉積產生影響。在算例7中，e/D=0.02，除了相對粗糙度較小引起顆粒的低攔截效率外，低粗糙度引起的低湍動能同樣對顆粒沉積產生影響。從圖6算例7的流場湍動能分布可以發(fā)現(xiàn)，流場湍動能的范圍均較小，顆粒沉積效率較低，導致相鄰2個粗糙元之間的顆粒沉積變化不大，形成的顆粒沉積峰不明顯。算例4雖然有較大的粗糙元高度，粗糙元壁面附近也引起了較大的湍動能，但并沒有導致較高的顆粒沉積數量。這是由于p/e較小，導致相鄰2個粗糙元之間只有一個回流區(qū)，顆粒難以進入，大多數顆粒隨空氣流出計算域，導致粗糙元間的顆粒沉積數量不高，從而出現(xiàn)分峰不明顯的現(xiàn)象。

顆粒在第1個粗糙元迎風面大量沉積后，在第2個粗糙元迎風面上的沉積卻大幅下降。其原因是，一方面，第2個粗糙元附近的流場受粗糙元的影響較?。ㄈ鐖D5所示），邊界層增長尚不大；另一方面，大部分顆粒被第1個粗糙元攔截沉積，導致第2個粗糙元的攔截沉積數量不大。同時，當X=0.3 m時，流場邊界層增長與湍動能開始趨于穩(wěn)定，之后顆粒沉積數量也達到穩(wěn)定狀態(tài)，基本不變。

4 結論

采用RSM與DPM模型，結合顆粒沉積的UDF程序，考慮顆粒的反彈與沉積，研究了二維半圓形粗糙元壁面顆粒沉積特性。得到以下結論：

（1）半圓形粗糙元壁面形成的回流區(qū)對顆粒沉積特性有很大影響，但回流區(qū)壁面通道較方形粗糙元小。

（2）半圓形粗糙元壁面顆粒沉積速度比同粗糙元參數下的方形粗糙元小，且顆粒粒徑越小，相差越大。

（3）半圓形粗糙元的高度與間距對顆粒沉積有一定影響。對于粒徑為1～10 μm的顆粒，當e/D=0.02，p/e=3時，顆粒沉積速度變化范圍較大；而在其他粗糙元參數下，顆粒沉積速度變化不大。

（4）半圓形粗糙元的迎風面是該管道顆粒沉積的主要位置，且顆粒沉積數量在第一個粗糙元迎風面上達到最大。

湍流對顆粒沉積的影響起至關重要的作用，將來可進一步對半圓形粗糙元壁面附近湍流脈動進行修正。本文僅討論了入口風速為5 m·s-1時的顆粒沉積特性，將來可進一步討論不同入口風速時的顆粒沉積特性。同時，為更加貼近工程實際，還需進一步開展更大粒徑范圍與三維結構通風管道數值計算的研究。