李 亮， 石瑞芳， 林建忠

（1.浙江大學(xué) 航空航天學(xué)院 流體工程研究所，杭州 310027；2.廣東盈峰智能環(huán)衛(wèi)科技有限公司，廣東 佛山 528322）

引 言

氣固兩相流中，管道氣力輸送顆粒具有占地面積小、無(wú)污染、設(shè)置靈活、輸送距離長(zhǎng)等優(yōu)點(diǎn)，因而廣泛應(yīng)用于機(jī)械制造、冶金、發(fā)電、材料工程、制藥和食品生產(chǎn)等各個(gè)行業(yè).氣力輸送過(guò)程中，顆粒在管道中的沉積會(huì)導(dǎo)致管道的堵塞，降低系統(tǒng)的效率甚至影響安全生產(chǎn).可見(jiàn)，研究氣力輸送過(guò)程中顆粒在管道中的輸運(yùn)與沉積對(duì)工程應(yīng)用具有重要意義.顆粒在管道壁面上的沉積與熱泳力[1]、顆粒慣性力、重力[2]、Brown和湍流擴(kuò)散[3]有關(guān)，是個(gè)復(fù)雜的過(guò)程.在以往的研究中，已有一些內(nèi)容涉及顆粒在直管中的輸運(yùn)和沉積[4-8].然而，在實(shí)際應(yīng)用中彎管的情形是很普遍的，彎管中的流體在離心力的作用下產(chǎn)生二次流，將顆粒輸送至壁面附近的區(qū)域；顆粒本身也會(huì)在離心力作用下甩向彎道的外壁，這將導(dǎo)致顆粒在壁面沉積率的增加.因此，有必要研究顆粒流經(jīng)彎管湍流場(chǎng)時(shí)的沉積特性.

目前，已有一些顆粒流經(jīng)彎管湍流場(chǎng)的研究結(jié)果.Pui等[9]實(shí)驗(yàn)研究了Reynolds數(shù)100＜Re＜10 000、顆粒直徑0.1 μm≤dp≤10 μm的情況下，顆粒流經(jīng)彎管時(shí)的通過(guò)率，發(fā)現(xiàn)通過(guò)率既不依賴于彎管的曲率也不依賴于Re.Balásházy等[10]研究了顆粒的慣性碰撞和重力沉降，發(fā)現(xiàn)在吸入的情況下管道橫截面中的二次流對(duì)沉積率沒(méi)有顯著影響.Lee和Gieseke[11]指出，對(duì)于顆粒直徑0.035 μm≤dp≤1.3 μm、Reynolds數(shù)1 800＜Re＜15 600的情形而言，當(dāng)顆粒的慣性沉降和Brown擴(kuò)散同時(shí)起作用時(shí)，僅用以往的理論無(wú)法令人滿意地預(yù)測(cè)顆粒的最小沉積率.Sato等[12]發(fā)現(xiàn)顆粒的沉積率隨Stk和De的增加而提高，隨彎道曲率的減小而降低.Wang等[13]指出，當(dāng)Re較低時(shí)，直徑為5 nm≤dp≤15 nm顆粒擴(kuò)散損失的增加對(duì)彎管的取向比較敏感.Yook和Pui[14]發(fā)現(xiàn)在Dean數(shù)21≤De≤1 779的情況下，直徑為3 nm≤dp≤50 nm顆粒的通過(guò)率隨著粒徑和De的增加而增加，而當(dāng)De＞200時(shí), 曲率對(duì)通過(guò)率的影響可以忽略.Lin等[15]說(shuō)明當(dāng)Schmidt數(shù)Sc較小時(shí)，顆粒在管道中的分布主要由軸向速度決定，而當(dāng)Sc遠(yuǎn)大于1的數(shù)量級(jí)時(shí)，二次流將主導(dǎo)顆粒的分布，顆粒的沉積區(qū)域隨著De的增加而變得均勻.Wilson等[16]指出，在一定的Stk范圍內(nèi)，Re的增加并不會(huì)顯著改變顆粒沉積率的變化趨勢(shì)，而在0.1≤Stk≤0.4范圍內(nèi)，顆粒沉積率隨Re的增加而顯著增加；當(dāng)Stk=0.15時(shí)，Re從10 250增加到30 750會(huì)導(dǎo)致沉積率從0.14增加到0.36.Ghaffarpasand等[17]發(fā)現(xiàn)，當(dāng)1 426≤De≤2 885時(shí)，直徑為3 nm≤dp≤17 nm顆粒的通過(guò)率隨著曲率的增加而增加，而對(duì)Re的變化卻不敏感.Lin等[18]指出，對(duì)直徑為8 nm≤dp≤550 nm的顆粒，De對(duì)通過(guò)率的影響依賴于Sc；存在一個(gè)臨界De，超過(guò)該臨界De時(shí)，通過(guò)率由增加變?yōu)闇p少，而這個(gè)臨界值依賴于Sc；De越大，顆粒的通過(guò)率越高.

以上所述的顆粒都是圓球形狀，而在實(shí)際應(yīng)用中，非圓球顆粒例如圓柱狀或橢球狀顆粒還是很常見(jiàn)的.非圓球顆粒在氣流中的輸運(yùn)比圓球的情形復(fù)雜，因?yàn)轭w粒的旋轉(zhuǎn)及其取向與顆粒的平動(dòng)存在耦合，而迄今為止，鮮有對(duì)非圓球顆粒流經(jīng)彎管湍流場(chǎng)時(shí)顆粒取向和沉積特性的研究.圓柱狀顆粒是非圓球顆粒中最典型的一類，因此，本文首先數(shù)值求解流體的平均運(yùn)動(dòng)、湍動(dòng)能、耗散率和脈動(dòng)速度方程，然后數(shù)值求解圓柱狀顆粒的運(yùn)動(dòng)和取向方程，得到不同參數(shù)下顆粒取向沿流向不同截面和出口處的分布以及顆粒流經(jīng)彎管時(shí)的沉積率.

1 圓柱狀顆粒基本方程

圓柱狀顆粒流經(jīng)一個(gè)彎管的流場(chǎng)如圖1所示，圖中S是位于xOy平面的彎管的中線，r和θ是定義在橫截面上的極坐標(biāo)，u,v和w分別是r, θ和S方向的速度分量，a和R分別是彎管內(nèi)徑和曲率.

圖1 彎管流場(chǎng)和坐標(biāo)系Fig.1 The flow field of the curved tube and the coordinate system

在固定坐標(biāo)系Oxyz上的圓柱狀顆粒的取向可以由圖2中的 φ和ψ確定，其中坐標(biāo)系O123固定在顆粒上，方向3沿著顆粒的長(zhǎng)軸，方向1是顆粒長(zhǎng)軸在xOy平面上的投影方向，方向2是垂直于1O3平面的方向，φ是顆粒在xOy平面上的投影與x軸之間的夾角.本文假設(shè)顆粒為剛性圓柱狀顆粒，繞顆粒的流動(dòng)是Stokes流，顆粒對(duì)流體的影響忽略不計(jì).

圖2 兩個(gè)坐標(biāo)系

1.1 圓柱狀顆粒數(shù)值計(jì)算模型

圓柱狀顆粒在流體中運(yùn)動(dòng)的受力和速度模型采用Batchelor提出的細(xì)長(zhǎng)體理論[19]，即由圓柱狀顆粒引起的流場(chǎng)擾動(dòng)速度與由Stokeslet線分布引起的擾動(dòng)速度相同，Stokeslet是Stokes流中的奇點(diǎn)施加到流體上的力，基于對(duì)Stokes方程的求解，可以得到該力以及由該力誘導(dǎo)的流場(chǎng)擾動(dòng)速度.在細(xì)長(zhǎng)體理論中，圓柱狀顆粒被分為若干段，每一段施加到流體上的力用一個(gè)Stokeslet點(diǎn)力表示，沿著圓柱狀顆粒主軸長(zhǎng)度積分，便可以得到圓柱狀顆粒誘導(dǎo)的流場(chǎng)擾動(dòng)速度[20]：

式中u是流體的瞬時(shí)速度；xc是顆粒質(zhì)心的坐標(biāo)；vp和ω分別是顆粒的速度和角速度；l和p分別是顆粒長(zhǎng)度和取向的單位矢量；s是沿顆粒主軸的無(wú)量綱坐標(biāo)，-1≤s≤1；β是顆粒的長(zhǎng)徑比；b(s)是在s處顆粒橫截面的形狀因子(對(duì)于圓柱狀顆粒b(s)=1)；δ是單位矩陣；f(s)是流體施加在顆粒上的力.方程(1)左邊的項(xiàng)是顆粒和流體在流場(chǎng)某一點(diǎn)的速度差.首先根據(jù)Gauss積分點(diǎn)將圓柱狀顆粒劃分為M段，然后在每一段上把f(s)變換成f(si)，最后，將方程(1)右邊積分轉(zhuǎn)換為一系列f(si)的線性求和，這樣就可以將方程(1)轉(zhuǎn)換為3M個(gè)線性方程組：

式中l(wèi)m是圓柱狀顆粒第m段的長(zhǎng)度，f(si)由求解方程(2)得到.通過(guò)對(duì)以下方程數(shù)值積分，便可獲得流體施加在圓柱狀顆粒上的合力F與合力矩L：

1.2 圓柱狀顆粒運(yùn)動(dòng)方程

圓柱狀顆粒在流體拖曳力、離心力和隨機(jī)力作用下的運(yùn)動(dòng)方程如下所示[21]，其中隨機(jī)力基于Stokes-Einstein的色散理論：

式中m是顆粒質(zhì)量；Rb是均值和方差為零的Gauss隨機(jī)分布的隨機(jī)向量；dpe是圓柱狀顆粒的體積當(dāng)量直徑；kB是Boltzmann常數(shù)；T是溫度；μ是流體的動(dòng)力黏度；Δt是隨機(jī)力作用的時(shí)間間隔；Ji是顆粒的慣性矩；Li是力矩分量；ωi是顆粒角速度分量；下標(biāo)1、2和3表示坐標(biāo)系O123中的3個(gè)分量.

以下方程用于將方程(5)中的力矩從坐標(biāo)系O123轉(zhuǎn)換為坐標(biāo)系OrθS：

1.3 圓柱狀顆粒的取向分布函數(shù)

圓柱狀顆粒的取向分布函數(shù)Ψ由Fokker-Planck方程控制[22-23]：

式中ω是顆粒的角速度；D是旋轉(zhuǎn)擴(kuò)散系數(shù)；Rr是顆粒旋轉(zhuǎn)的阻力系數(shù)[24]；l是顆粒長(zhǎng)度；ρ和ν分別是流體的密度和運(yùn)動(dòng)黏性系數(shù).方程(7)右邊的第一項(xiàng)表示流體作用力的影響，第二項(xiàng)是隨機(jī)力影響.圓柱狀顆粒的取向由分布函數(shù)Ψ描述，它表示顆粒取向位于某一角度范圍內(nèi)的概率.

1.4 顆粒碰撞模型

間步長(zhǎng)為Δt= (u-vi)/， 這里u-是顆粒和流體的速度差顆粒在時(shí)間步長(zhǎng)i的加速度.在下一時(shí)間步

當(dāng)顆粒流過(guò)彎管時(shí)，顆粒會(huì)發(fā)生相互碰撞且有可能和管壁發(fā)生碰撞，這種碰撞會(huì)影響顆粒的平移和旋轉(zhuǎn)速度從而影響顆粒的分布和通過(guò)率.如圖3所示，假設(shè)兩個(gè)圓柱狀顆粒的碰撞為瞬時(shí)、非完全彈性碰撞，圖中兩顆粒接觸點(diǎn)及其法線方向n由兩顆粒的相對(duì)位置決定.本文圓柱狀顆粒基于細(xì)長(zhǎng)體理論建模，故一個(gè)顆粒被分成M段，每一段被視為一個(gè)單元，每個(gè)單元參與顆粒間碰撞的判斷.在對(duì)方程(4)、(5)的積分過(guò)程中，變時(shí)i+ 1，如果顆粒1的i段質(zhì)心與顆粒2的j段質(zhì)心之間的距離在dp-dp/10到dp+dp/10的范圍之內(nèi)，則認(rèn)為兩個(gè)顆粒發(fā)生碰撞.如果距離小于dp-dp/10，則原始時(shí)間步長(zhǎng)減少一半并重新計(jì)算.假設(shè)碰撞點(diǎn)O位于顆粒1的i段質(zhì)心和顆粒2的j段質(zhì)心連線的中點(diǎn)，l1和l2分別是從兩個(gè)顆粒質(zhì)心O1和O2到接觸點(diǎn)O的矢量，從而定義法向矢量n=l1+l2.碰撞時(shí)兩個(gè)顆粒沿法線方向獲得一個(gè)沖量，碰撞后兩個(gè)顆粒的平動(dòng)速度和角速度取決于該沖量：

圖3 兩顆粒碰撞示意圖Fig.3 Schematic of the collision between 2 particles

式中m和vp分別是顆粒質(zhì)量和速度；下標(biāo)1、2表示顆粒1和顆粒2；上標(biāo)“′”表示碰撞后.基于碰撞定律得

式中e是恢復(fù)系數(shù)；vp1o和vp2o是碰撞前接觸點(diǎn)處兩個(gè)顆粒沿法向的速度分量.作用在兩個(gè)顆粒上的力矩分別是In×l1和-In×l2，于是顆粒的旋轉(zhuǎn)方程為

那么沖量I可以寫成

當(dāng)顆粒1與壁面碰撞時(shí)，顆粒2可視為壁面，此時(shí)m2和vp2分別視為無(wú)窮大和零.

2 流場(chǎng)基本方程

要解方程(2) ～ (8)，首先要求解流場(chǎng)得到流場(chǎng)的瞬時(shí)速度u，該速度由平均速度U和脈動(dòng)速度u′組成.

2.1 平均運(yùn)動(dòng)方程

假設(shè)流動(dòng)為不可壓且充分發(fā)展的湍流，連續(xù)性方程和Reynolds平均運(yùn)動(dòng)方程為

式中Ui是平均速度，由3個(gè)分量Ur，Uθ和US構(gòu)成；P是平均壓力；ρ和ν分別是流體密度和運(yùn)動(dòng)黏性系數(shù)；是Reynolds應(yīng)力.定義平均軸向速度為

式中B=-?P/?S；μ是流體動(dòng)力黏性系數(shù)；a是管道半徑.定義無(wú)量綱參數(shù)如式(16)所示，用其對(duì)方程(14)無(wú)量綱化得方程(17)：

2.2 Reynolds應(yīng)力方程

方程(14)中的Reynolds應(yīng)力由以下方程描述：

式中p′是脈動(dòng)壓力.方程(18)右邊第一項(xiàng)的湍流擴(kuò)散項(xiàng)可以?；癁閇25]

以上方程中包含的k和ε的修正方程為

式中Cε1=1.44,Cε2=1.92, σk2=1.0, σε=1.3.

2.3 流體脈動(dòng)速度

湍流擴(kuò)散由流體脈動(dòng)速度u′體現(xiàn)，一個(gè)湍流場(chǎng)可以用Fourier模式、單點(diǎn)速度關(guān)聯(lián)表達(dá)式和標(biāo)量能譜表示.均勻各向同性湍流中的速度場(chǎng)可以模擬不同模態(tài)之間的動(dòng)力相互作用以及渦度場(chǎng)對(duì)流的運(yùn)動(dòng)學(xué)過(guò)程.Kolmogorov局部各向同性湍流的假設(shè)說(shuō)明，小尺度湍流運(yùn)動(dòng)在統(tǒng)計(jì)上是各向同性的，流體的脈動(dòng)速度u′與小尺度湍流運(yùn)動(dòng)有關(guān).因此，可以選擇均勻各向同性湍流模型來(lái)得到流體的脈動(dòng)速度.為了得到流體的脈動(dòng)速度u′，基于動(dòng)力學(xué)模擬掃掠模型[28-29]，用Fourier級(jí)數(shù)來(lái)表示u′：

式中ξ(n)和ζ(n)是滿足平均值為零、標(biāo)準(zhǔn)差為的Gauss分布的隨機(jī)矢量，其中N是常數(shù)，本文取N=100；頻率ω(n)是標(biāo)準(zhǔn)差為u′的均方根的Gauss隨機(jī)數(shù)；波數(shù)矢量k(n)是單位球面上各向同性的隨機(jī)矢量.ω(n)和半徑k(n)的分布函數(shù)由下式確定：

式中E(k)是 Karman-Pao湍流能譜；D(τ)是二階關(guān)聯(lián)函數(shù).那么ω(n)和k(n)服從Cauchy分布：

3 數(shù)值模擬

3.1 求解步驟

1)求解方程(13) ～ (24)得到Ui,k和ε；

2) 求解方程(25) ～ (28)得到u′；

3)對(duì)顆粒位置、取向、速度和角速度進(jìn)行初始化；

4)由方程(2)、(3)計(jì)算F和L；

5)由方程(6)將L從坐標(biāo)系Oxyz變換到坐標(biāo)系O123；

6)由方程(4)、(5)和方程(7)、(8)計(jì)算下一時(shí)間步的顆粒位置和取向；

7)返回步驟4)，直到顆粒流出管道或沉積在管壁上；

8)計(jì)算顆粒的沉積率Dep=(Nin-Nout)/Nin，這里Nout和Nin分別是管道出口和入口的顆粒數(shù).

3.2 數(shù)值模擬方法和相關(guān)參數(shù)

用有限體積法求解方程(7)、(8)和方程(13)、(14)，采用SIMPLE格式處理速度-壓力耦合項(xiàng)，對(duì)流項(xiàng)用冪律格式離散，該格式是一維對(duì)流-擴(kuò)散型方程精確解的分段近似，能給出物理真實(shí)解.模擬中采用交錯(cuò)網(wǎng)格，壓力、軸向速度US位于控制體中心，而橫截面上的速度分量Ur和Uθ位于控制體的邊界.壁面采用無(wú)滑移邊界條件，標(biāo)準(zhǔn)壁面函數(shù)用于近壁區(qū)域的計(jì)算，最靠近壁面單元的中心與壁面的距離為y+=30.用顯式Euler公式積分方程(4)、(5)，得到下一時(shí)間步顆粒的速度和角速度.數(shù)值模擬中的相關(guān)參數(shù)為ρ=1.205 kg/m3，μ=1.808 ×10-5Pa·s，T=293 K，kB=1.38 × 10-23J/K，e=0.37.20 000個(gè)圓柱狀顆粒初始均勻地分布于管道的入口，其取向沿各方向也是均勻分布，對(duì)中心線長(zhǎng)度Lb=πR/2的90°彎曲管進(jìn)行計(jì)算.

文中的Re，De和曲率的定義如式(16)所示，Stk定義為顆粒響應(yīng)時(shí)間與流動(dòng)特征時(shí)間之比：

式中ρp是顆粒密度；uτ是壁摩擦速度.

3.3 網(wǎng)格獨(dú)立性及數(shù)值模擬方法驗(yàn)證

網(wǎng)格單元的數(shù)量為128(r) × 128(θ) × 256(S)，在θ和S方向采用均勻網(wǎng)格、沿r方向的網(wǎng)格越靠近壁面越密.將r和θ方向的網(wǎng)格點(diǎn)從112變化為144、S方向網(wǎng)格點(diǎn)從240變化為272，進(jìn)行網(wǎng)格的獨(dú)立性驗(yàn)證，把每個(gè)物理量的所有歸一化剩余誤差小于10-4作為計(jì)算收斂的準(zhǔn)則.

為了驗(yàn)證數(shù)值模型和方法，將本文數(shù)值模擬得到的管道出口處流場(chǎng)截面上沿水平中線上的平均軸向速度分布和軸向脈動(dòng)速度的均方根與實(shí)驗(yàn)結(jié)果[31]進(jìn)行比較，結(jié)果如圖4所示，可見(jiàn)數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果基本相符.

圖4 平均軸向速度和軸向脈動(dòng)速度均方根的分布(Re=10 500, De=2 460)：（a）平均軸向速度；（b）軸向脈動(dòng)速度均方根Fig.4 Distributions of the mean axial velocity and the RMS value of the fluctuating axial velocity (Re=10 500, De=2 460): (a) the mean axial velocity;(b) the RMS value of the fluctuating axial velocity

4 結(jié)果與討論

4.1 顆粒取向在不同軸向位置的分布

不同軸向位置處橫截面上顆粒平均取向分布如圖5所示，圖中橫坐標(biāo) ?是顆粒的主軸與流場(chǎng)當(dāng)?shù)亓骶€的夾角，通過(guò)對(duì)固定軸向位置橫截面上所有顆粒的取向角進(jìn)行平均得到，縱坐標(biāo)是概率.流場(chǎng)中的速度矢量用速度的三個(gè)分量值計(jì)算得到，顆粒的取向矢量由求解取向分布函數(shù)方程(7)得到.由圖可見(jiàn)，顆粒的平均取向從入口的各向同性分布沿流動(dòng)方向變?yōu)榉歉飨蛲苑植?，顆粒的取向逐漸趨向于流動(dòng)方向.由方程(5) ～ (8)可知，顆粒的取向由流體施加在顆粒上的力矩和Brown擴(kuò)散所控制，后者使顆粒取向分布變得更加各向同性，那么顆粒取向趨向于流動(dòng)方向是由前者的力矩所造成，力矩一方面使顆粒繞渦量軸旋轉(zhuǎn)，另一方面使顆粒旋轉(zhuǎn)后的取向趨向于流動(dòng)方向.

圖5 橫截面上顆粒平均取向分布(Re=30 000, Stk = 1, De=2 200，β=8)Fig.5 Distributions of mean orientations of particles on the cross section (Re =30 000, Stk =1, De=2 200，β=8)

4.2 顆粒取向在出口截面處的分布

出口處橫截面上Stk和De對(duì)顆粒平均取向分布的影響如圖6、7所示.由圖6可見(jiàn)，隨著Stk的增加，顆粒取向趨向于流動(dòng)方向的現(xiàn)象更明顯.因?yàn)槿绶匠?7)、(8)所示，顆粒的取向由流體阻力和Brown隨機(jī)力控制，Brown隨機(jī)力使顆粒的取向趨向均勻分布，削弱其占優(yōu)取向的趨勢(shì)，而對(duì)于Stk較大的顆粒而言，Brown隨機(jī)力的作用相對(duì)較弱.

圖6 不同Stk時(shí)，顆粒平均取向分布(Re=30 000, De=2 200, β=8)Fig.6 Distributions of particle orientations for different Stk values(Re =30 000, De =2 200, β=8)

從圖7可以看出，隨著De的增加，顆粒取向趨向于流動(dòng)方向的現(xiàn)象變得不明顯，因?yàn)镈e與曲率κ成正比(見(jiàn)式(16))，而彎管內(nèi)二次流的強(qiáng)度與曲率κ成正比，大的De意味著大的二次流強(qiáng)度，二次流施加在顆粒上的力起著阻止顆粒取向趨向于流動(dòng)方向的作用.

圖7 不同De時(shí)，顆粒平均取向分布(Re =30 000, Stk=1, β=8)Fig.7 Distributions of particle orientations for different De values(Re=30 000, Stk=1, β=8)

出口橫截面上Re和顆粒長(zhǎng)徑比β對(duì)顆粒平均取向分布的影響如圖8、9所示.由圖8可見(jiàn)，隨著Re的增大，顆粒取向趨向于流動(dòng)方向的現(xiàn)象減弱.其原因有兩個(gè)方面，方程(2)、(4)表明，流體作用在顆粒上的力矩是由顆粒速度和流體瞬時(shí)速度之差產(chǎn)生的，而流體瞬時(shí)速度由平均速度U和脈動(dòng)速度u′組成，平均速度及其所誘導(dǎo)的剪切率使顆粒旋轉(zhuǎn)并導(dǎo)致取向趨向于流動(dòng)方向，而隨機(jī)的脈動(dòng)速度導(dǎo)致顆粒旋轉(zhuǎn)并使其取向趨向于各向同性分布.Reynolds數(shù)是表征湍流的重要參數(shù)，Reynolds數(shù)越大，最小渦的尺度越小，不同尺度的旋渦所包含的能譜分布越廣，脈動(dòng)速度的影響越大，于是顆粒取向趨向于流動(dòng)方向的現(xiàn)象越不明顯.另一方面，當(dāng)Re較大時(shí)，流動(dòng)中會(huì)出現(xiàn)較強(qiáng)的二次流，這也會(huì)阻礙顆粒的取向趨向于流動(dòng)方向.

圖8 不同Re時(shí)，顆粒取向分布(De=2 200, Stk=1, β=8)Fig.8 Distributions of particle orientations for different Re values(De=2 200, Stk=1, β=8)

在圖9中，隨著顆粒長(zhǎng)徑比β的增加，顆粒取向更趨向于流動(dòng)方向.流體作用在顆粒上的力矩與渦的尺度有關(guān)，即大尺度和小尺度的渦分別對(duì)應(yīng)于平均速度U和脈動(dòng)速度u′，長(zhǎng)徑比小的顆粒受脈動(dòng)速度u′的影響較大，于是顆粒取向趨向于流動(dòng)方向的現(xiàn)象不明顯.由圖還可見(jiàn)，在大長(zhǎng)徑比的情況下，β對(duì)顆粒平均取向分布的影響較小，這與Krushkal和Gallily[32]得出的結(jié)論一致，該結(jié)論認(rèn)為當(dāng)顆粒長(zhǎng)徑比大于5時(shí)，顆粒的取向分布對(duì)長(zhǎng)徑比的變化不敏感.

圖9 不同β時(shí)，顆粒取向分布(Re=30 000, De =2 200, Stk=1)Fig.9 Distributions of particle orientations for different β values(Re=30 000, De =2 200, Stk=1)

4.3 Stk和顆粒長(zhǎng)徑比對(duì)顆粒沉積率的影響

本文采用的細(xì)長(zhǎng)體理論將圓柱狀顆粒分成M段，每一段都被視為采樣單元作為與壁面碰撞的判斷依據(jù)，如果某一段的質(zhì)心與壁面之間的距離小于dp/2，則確定顆粒碰到壁面.當(dāng)撞擊壁面的顆粒碰撞速度較低時(shí)，顆粒將因不足以逃脫壁面吸引力而沉積在壁面上.模擬中當(dāng)顆粒碰撞后垂直于壁面的速度分量小于或等于零時(shí)，就認(rèn)為顆粒沉積在壁面上.如方程(4)、(5)所示，顆粒的沉積數(shù)量取決于顆粒軌跡，而該軌跡與作用在顆粒上的力、力矩和顆粒長(zhǎng)徑比有關(guān).二次流會(huì)增強(qiáng)顆粒向壁面的傳輸從而提高顆粒的沉積率.

為了驗(yàn)證數(shù)值模擬顆粒沉積率結(jié)果的可靠性，本文計(jì)算了當(dāng)顆粒長(zhǎng)徑比β=1時(shí)在不同Stk情況下顆粒流經(jīng)彎管的通過(guò)率Pe=Nout/Nin，這里Nout和Nin是管道出口和入口的顆粒數(shù)，通過(guò)率Pe和沉積率Dep的關(guān)系為Dep=1 -Pe.圖10是數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果[17]的比較，可見(jiàn)兩者吻合較好.

圖10 不同Stk下的顆粒通過(guò)率(De=1 862, Re=10 500, β=1)Fig.10 The pass ratio of particles at different Stk values (De=1 862,Re=10 500, β=1)

圖11是不同顆粒長(zhǎng)徑比(β)情況下沉積率與Stk的關(guān)系，可見(jiàn)沉積率并非Stk的單調(diào)函數(shù)，當(dāng)Stk較小時(shí)，隨機(jī)力對(duì)顆粒的影響大于慣性力的影響，近壁區(qū)的顆粒較容易被隨機(jī)力驅(qū)使到達(dá)壁面，導(dǎo)致較大的沉積率.隨著Stk的增加，隨機(jī)力的作用逐漸減弱，近壁區(qū)中的顆粒到達(dá)壁面的數(shù)量減少，沉積率也逐漸降低，這與以往的圓球顆粒的通過(guò)率隨著顆粒直徑的增加而增加的結(jié)論[14]一致，因?yàn)镾tk與顆粒直徑成正比.當(dāng)Stk進(jìn)一步增加到Stk≈7時(shí)，慣性力的作用大于隨機(jī)力的作用，較多的顆粒朝壁面附近運(yùn)動(dòng)并到達(dá)壁面，導(dǎo)致顆粒的沉積率有所回升，這與圓球情形的結(jié)論[12]是吻合的.

由圖11還可知，顆粒沉積率隨顆粒長(zhǎng)徑比β的增加而增加，其原因可分析如下：長(zhǎng)徑比越大的顆粒承受的慣性離心力越大，因而越容易在靠近彎管外側(cè)壁面的區(qū)域聚集，這為顆粒的沉積提供了條件.另一方面，具有較大長(zhǎng)徑比顆粒的旋轉(zhuǎn)導(dǎo)致顆粒更容易到達(dá)并沉積在壁面上.上述結(jié)論與直管情況下的結(jié)論一致，即剛性纖維的沉積率比圓球顆粒的沉積率高[33].

圖11 不同顆粒長(zhǎng)徑比時(shí)沉積率與Stk的關(guān)系(Re=30 000, De=2 200)Fig.11 The relationship between the deposition rate and Stk at different particle aspect ratios (Re=30 000, De=2 200)

4.4 Re和De對(duì)顆粒沉積率的影響

圖12和13分別給出了顆粒沉積率與Re和De的關(guān)系，可見(jiàn)沉積率隨著Re和De的增加而增加，這與圓球情形[33]相同.Re和De較大的情況下，顆粒更容易被較強(qiáng)的二次流帶向壁面，導(dǎo)致顆粒碰撞壁面從而沉積在壁面的可能性增大，這與式(16)中De與Re成正比的定義是吻合的.比較圖12和13可知，在本文計(jì)算的Re和De范圍內(nèi)，Re對(duì)顆粒沉積率的影響比De的影響明顯，在式(16)中，De與Re的關(guān)系相差一個(gè)曲率的因子κ，而κ小于1，所以這一結(jié)論也與式(16)中的定義吻合.

圖12 不同Re下沉積率與Stk的關(guān)系(De=2 200, β=8)Fig.12 The relationship between the deposition rate and Stk at different Re values (De =2 200, β=8)

圖13 不同De下沉積率與Stk的關(guān)系(Re =30 000, β=8)Fig.13 The relationship between the deposition rate and Stk at different De values (Re =30 000, β=8)

5 結(jié) 論

本文對(duì)圓柱狀顆粒流經(jīng)彎管湍流場(chǎng)時(shí)的取向與沉積特性進(jìn)行了研究，首先求解流體平均運(yùn)動(dòng)、湍動(dòng)能、耗散率和脈動(dòng)速度方程，然后求解顆粒運(yùn)動(dòng)方程和取向的Fokker-Planck方程，得到顆粒取向在不同軸向位置和出口截面處的分布以及顆粒流經(jīng)彎管時(shí)的沉積率，一些數(shù)值模擬結(jié)果與相關(guān)的實(shí)驗(yàn)進(jìn)行了比較，討論了各參數(shù)對(duì)顆粒取向以及沉積率的影響.研究結(jié)果表明，隨著Stk和顆粒長(zhǎng)徑比β的增加以及De和Re的減少，顆粒的主軸更趨向于流動(dòng)方向.顆粒流經(jīng)彎管時(shí)的沉積率隨著De，Re和顆粒長(zhǎng)徑比β的增加而增加，但隨Stk的變化呈現(xiàn)非單調(diào)趨勢(shì)，當(dāng)Stk＜7時(shí)，沉積率隨Stk的增加而減少，當(dāng)Stk＞7時(shí)則相反.

圓柱狀顆粒流經(jīng)彎管湍流場(chǎng)時(shí)的取向分布和沉積率是一個(gè)復(fù)雜問(wèn)題，本研究的數(shù)值模擬結(jié)果與圓球顆粒的結(jié)果進(jìn)行了一些比較，開(kāi)展圓柱狀顆粒流經(jīng)彎管湍流場(chǎng)的實(shí)驗(yàn)研究是將來(lái)需要進(jìn)行的工作.

參考文獻(xiàn)（ References ）：

［1］AKSHAT T M, MISRA S, GUDIYAWAR M Y, et al.Effect of electrospun nanofiber deposition on thermophysiology of functional clothing[J].Fibers and Polymers, 2019, 20（5）: 991-1002.

［2］TIAN L, AHMADI G, WANG Z C, et al.Transport and deposition of ellipsoidal particles in low Reynolds number flows[J].Journal of Aerosol Science, 2012, 45: 1-18.

［3］TU C X, YIN Z Q, LIN J Z, et al.A review of experimental techniques for measuring micro- to nano-particleladen gas flows[J].Applied Sciences, 2017, 7（2）.DOI: 10.3390/app7020120.

［4］SUN L, LIN J Z, BAO F B.Numerical simulation on the deposition of nanoparticles under laminar conditions[J].Journal of Hydrodynamics, 2006, 18（6）: 676-680.

［5］PHARES D J, SHARMA G.A DNS study of aerosol deposition in a turbulent square duct flow[J].Aerosol Science and Technology, 2006, 40（11）: 1016-1024.

［6］ARMAND P, BOULAUD D, POURPRIX M, et al.Two-fluid modeling of aerosol transport in laminar and turbulent flows[J].Journal of Aerosol Science, 1998, 29（8）: 961-983.

［7］羅萬(wàn)清, 李海燕, 梁劍寒.基于Euler-Lagrange模型的電弧風(fēng)洞噴管兩相流模擬[J].應(yīng)用數(shù)學(xué)和力學(xué), 2020, 41（1）:16-26.（LUO Wanqing, LI Haiyan, LIANG Jianhan.Simulation of 2-phase flow in the nozzle of the arc heated wind tunnel based on the Eulerian-Lagrange model[J].Applied Mathematics and Mechanics, 2020, 41（1）: 16-26.（in Chinese））

［8］LIN J Z, YIN Z Q, GAN F J, et al.Penetration efficiency and distribution of aerosol particles in turbulent pipe flow undergoing coagulation and breakage[J].International Journal of Multiphase Flow, 2014, 61: 28-36.

［9］PUI D Y H, ROMAY-NOVAS F, LIU B Y H.Experimental study of particle deposition in bends of circular cross section[J].Aerosol Science and Technology, 1987, 7（3）: 301-315.

［10］BALáSHáZY I, MARTONEN T B, HOFMANN W.Inertial impaction and gravitational deposition of aerosols in curved tubes and airway bifurcations[J].Aerosol Science and Technology, 1990, 13（3）: 308-321.

［11］LEE K W, GIESEKE J A.Deposition of particles in turbulent flow pipes[J].Journal of Aerosol Science, 2006,25（4）: 699-709.

［12］SATO S, CHEN D R, PUI D Y H.Particle transport at low pressure: deposition in bends of a circular cross-section[J].Aerosol Science and Technology, 2003, 37: 770-779.

［13］WANG J, FLAGAN R C, SEINFELD J H.Diffusional losses in particle sampling systems containing bends and elbows[J].Journal of Aerosol Science, 2002, 33（6）: 843-857.

［14］YOOK S J, PUI D Y H.Experimental study of nanoparticle penetration efficiency through coils of circular crosssections[J].Aerosol Science and Technology, 2006, 40（6）: 456-462.

［15］LIN J Z, LIN P F, CHEN H J.Research on the transport and deposition of nanoparticles in a rotating curved pipe[J].Physics of Fluid, 2009, 21（12）: 122001.

［16］WILSON S R, LIU Y A, MATIDA E A, et al.Aerosol deposition measurements as a function of Reynolds number for turbulent flow in a ninety-degree pipe bend[J].Aerosol Science and Technology, 2011, 45（3）: 364-375.

［17］GHAFFARPASAND O, DREWNICK F, HOSSEINIEBALAM F, et al.Penetration efficiency of nanometer-sized aerosol particles in tubes under turbulent flow conditions[J].Journal of Aerosol Science, 2012, 50: 11-25.

［18］LIN J Z, YIN Z Q, LIN P F, et al.Distribution and penetration efficiency of nanoparticles between 8～550 nm in pipe bends under laminar and turbulent flow conditions[J].International Journal of Heat and Mass Transfer,2015, 85: 61-70.

［19］BATCHELOR G K.Slender-body theory for particles of arbitrary cross-section in Stokes flow[J].Journal of Fluid Mechanics, 1970, 44（3）: 419-440.

［20］MACKAPLOW M B, SHAQFEH E S G.A numerical study of the sedimentation of fiber suspension[J].Journal of Fluid Mechanics, 1988, 376（1）: 149-182.

［21］MICHAELIDES E E.Brownian movement and thermophoresis of nanoparticles in liquids[J].International Journal of Heat and Mass Transfer, 2015, 81: 179-187.

［22］LEAL L G, HINCH E J.The effect of weak Brownian rotations on particles in shear flow[J].Journal of Fluid Mechanics, 1971, 46: 685-703.

［23］高振宇, 林建忠, 李俊.纖維懸浮剪切湍流中纖維旋轉(zhuǎn)擴(kuò)散系數(shù)的理論研究[J].應(yīng)用數(shù)學(xué)和力學(xué), 2007, 28（3）: 263-269.（GAO Zhenyu, LIN Jianzhong, LI Jun.Theoretical research on the rotational dispersion coefficient of fiber in the turbulent shear flow of fiber suspension[J].Applied Mathematics and Mechanics, 2007, 28（3）: 263-269.（in Chinese））

［24］CHEN S B, JIANG L.Orientation distribution in a dilute suspension of fibers subject to simple shear flow[J].Physics of Fluid, 1999, 11（10）: 2878-2890.

［25］LIEN F S, LESCHZINER M A.Assessment of turbulence-transport models including non-linear RNG eddy-viscosity formulation and second-moment closure for flow over a backward-facing step[J].Computers and Fluids,1994, 23（8）: 983-1004.

［26］LAUNDER B E.Second-moment closure and its use in modeling turbulent industrial flows[J].International Journal for Numerical Methods in Fluids, 1989, 9（8）: 963-985.

［27］LAUNDER B E.Second-moment closure: present and future[J].International Journal of Heat and Fluid Flow,1989, 10（4）: 282-300.

［28］FUNG J C H, HUNT J C R, MALIK N A, et al.Kinematic simulation of homogeneous turbulence by unsteady random Fourier modes[J].Journal of Fluid Mechanics, 1992, 236: 281-318.

［29］WANG L P, STOCK D E.Numerical simulation of heavy particle dispersion-scale ration and flow decay considerations[J].Journal of Fluids Engineering, 1994, 116（1）: 154-163.

［30］DONG S, FENG X, SALCUDEAN M, et al.Concentration of pulp fibers in 3D turbulent channel flow[J].International Journal of Multiphase Flow, 2003, 29（1）: 1-21.

［31］WEBSTER D R, HUMPHREY J A C.Experimental observations of flow instability in a helical coil[J].Journal of Fluids Engineering, 1993, 115（3）: 436-443.

［32］KRUSHKAL E M, GALLILY I.On the orientation distribution function of nonspherical aerosol particles in a general shear flow, Ⅱ: the turbulent case[J].Journal of Aerosol Science, 1988, 19（2）: 197-211.

［33］PODGóRSKI A, GRADO? L, GRZYBOWSKI P.Theoretical-study on deposition of flexible and stiff fibrous aerosol-particles on a cylindrical collector[J].Chemical Engineering Journal and the Biochemical Engineering Journal, 1995, 58（2）: 109-121.