基于分形理論的駐波聲場中顆粒團運動特性數值預測?

賈文龍 凡鳳仙 蘇明旭

(上海理工大學能源與動力工程學院 上海 200093)

0 引言

懸浮在大氣中空氣動力學直徑小于或等于2.5μm的顆粒物被稱為PM2.5，其極易富集重金屬和多環芳烴等強致癌污染物，且能夠經呼吸系統進入人體，對人群健康帶來嚴重危害；同時，PM2.5影響大氣環境，加劇酸雨、光化學煙霧以及霧霾的形成，威脅到人類賴以生存的環境[1?2]。文獻報道顯示，燃煤電站、交通車輛、工業過程的排放是PM2.5的主要來源[3?5]。由于PM2.5粒徑細微，傳統的除塵方式，如慣性除塵、濕式除塵、過濾式除塵、靜電除塵等對PM2.5的清除均難以奏效，使得大量PM2.5排放到大氣中。隨著環境空氣質量標準的日趨嚴格，PM2.5排放控制技術研究顯得尤為重要且迫切。為了控制PM2.5的排放，一些學者提出了PM2.5凝并預處理技術，以期通過物理或化學作用促進PM2.5長大為粒徑較大的顆粒，從而提高傳統除塵裝置的效率。聲凝并是一種重要的PM2.5預處理技術，其原理是利用外加聲場作用促進PM2.5發生相對運動、碰撞接觸進而凝并為粒徑較大的顆粒物[6?7]。該項技術因裝置簡單、適應性強，而備受研究者關注[8?11]。

為探討聲凝并的內在機理和動力學過程，一些學者對聲場中單顆粒的運動[12?13]、顆粒間相互作用[14?17]、顆粒碰撞[18]和凝并[9,19?23]開展了一系列理論和數值模擬研究。然而，這些研究主要針對球形顆粒，甚至在對聲凝并進行建模時將固體顆粒聲凝并后生成的顆粒團也視為球形[9,20,22?23]。事實上，聲凝并形成的顆粒團普遍具有分形結構[19,21,23?24]。雖然研究者利用分形維數描述顆粒團形狀[19,21]，并將分形維數引入顆粒凝并核函數以建立聲凝并模型，模擬聲凝并過程中顆粒粒徑分布的演變，但模型中未能反映出顆粒團分形結構對其運動特性的影響。趙兵等[24]對顆粒團在駐波聲場中的運動軌跡進行可視化實驗并對顆粒團弛豫時間進行理論分析，但缺少對聲場中顆粒團運動規律的深入探究。楊旭峰等[25]雖然對顆粒團在聲場中的動力學特性進行過探討，但所針對的是直鏈顆粒團。在聲凝并的實際應用中，一旦顆粒凝并形成顆粒團，顆粒團的運動特性將對聲凝并過程起決定作用。因此，欲全面掌握聲凝并過程中顆粒動力學行為的細節信息，很有必要對聲場中顆粒團的運動特性進行研究。

本文將基于分形理論，利用顆粒團的質量-半徑關系式和無量綱流體動力學半徑-分形維數關系式，求解顆粒團的流體動力學半徑和分形維數，進而建立駐波聲場中顆粒團動力學模型，并將數值模擬結果與文獻中的實驗結果進行對比，驗證模型的準確性。在此基礎上，利用數值模擬方法，預測組成顆粒團的原生顆粒半徑、數目和排列情況對顆粒團在聲場中運動特性的影響，以揭示顆粒團的行為規律，并為建立更精確的聲凝并模型提供理論基礎。

1 數學模型與數值計算方法

1.1 分形理論對顆粒團的描述

圖1為球形顆粒組成的顆粒團的結構及特征參數示意圖，圖中ai為組成顆粒團的第i個原生顆粒的半徑，R為顆粒團半徑，即能夠包含顆粒團的最小球體的半徑；Rh為顆粒團的流體動力學半徑，即與顆粒團質量和動力學特性相同的非滲透性球形顆粒的半徑。

圖1 顆粒團結構及特征參數示意圖Fig.1 Schematic diagram of agglomerate structure and characteristic parameters

組成顆粒團的原生顆粒數目N、原生顆粒半徑a與顆粒團的流體動力學半徑Rh的函數關系可由質量-半徑關系式給出[26?27]：

式(1)中，Df為顆粒團的分形維數，表示顆粒團的空間填充致密程度，取值范圍為1

當原生顆粒半徑不相等時，式(1)中的原生顆粒半徑a可采用體積平均半徑aave代替，即

式(2)中，Vt是組成顆粒團的原生顆粒的總體積，

顆粒團的流體動力學半徑Rh與顆粒團半徑R的比值可寫成分形維數Df的函數[27?28]，即

聯立式(1)、式(4)，通過逐次逼近的迭代算法，即可獲得顆粒團的流體動力學半徑和分形維數。已有研究表明，顆粒團的質量-半徑關系式和無量綱流體動力學半徑-分形維數關系式，適用于描述原生顆粒數目為2～100的顆粒團[28?30]。

1.2 聲波的波動方程

將水平方向設為x向，由Navier-Stokes方程可推導出一維水平平面駐波聲場的波動方程為

式(5)中，u為聲波引起的氣體振動速度；u為速度振幅；k為波數，k=ω/c，c為聲速，ω=2πf，f為聲波頻率；x為聲波波動方向位置坐標；t為時間。

通常采用聲壓級來描述聲場的強度，聲壓級的表達式為

式(6)中，L為聲壓級；ρg為氣體密度；Pr為參考聲壓，Pr=2×10?5Pa。

1.3 顆粒團的動力學模型

為著重探討水平駐波聲場中單個顆粒團的動力學特性，忽略顆粒團的破碎和重組。聲場中顆粒團所受作用力包括重力、浮力、Stokes力和非穩定力(Basset力、虛擬質量力、壓力梯度力等)，對于氣相中的固體顆粒，與Stokes力相比，非穩定力可以忽略不計[13]。因此，若將重力方向設為y向，顆粒團的運動方程可寫為

其中，m為顆粒團質量，m=ρVt；ρ為顆粒材料密度；v為顆粒團速度，下標x、y表示在x向和y向的分量；μg為氣體動力黏度；g為重力加速度；Cc為Cunningham修正系數，其表達式為[9,18]

式(9)中，Kn為Knudsen數，Kn=λg/Rh，λg為氣體分子平均自由程。

1.4 數值計算方法

依據式(7)、式(8)給出的顆粒團運動方程，采用四階變步長Runge-Kutta算法計算出經歷一個時間步長后顆粒團的速度。在一個時間步長的初速度和末速度已知的前提下，采用二階隱式Adams插值算法求解顆粒團的位移，即

其中，X和Y分別為顆粒團的x向和y向位移；?t為時間步長。?t通常要比聲波周期T=1/f和顆粒團弛豫時間τ小得多，其中弛豫時間τ可寫為

2 數值模擬結果

2.1 模型驗證

趙兵等[24]曾采用高速顯微攝像系統對水平駐波聲場中單個顆粒團的運動軌跡進行可視化實驗。圖2給出了趙兵等[24]的實驗及本文數值模擬所得到的顆粒團運動軌跡。圖2(b)給出的數值模擬結果中顆粒團的尺寸采用流體動力學直徑(2Rh)表示。數值模擬采用的參數條件與實驗一致，具體數值見表1。其中，T為氣體溫度，p為氣體靜壓，x0為顆粒初始位置。根據顆粒團的組成和結構參數(a、N、Df)，利用式(1)計算得到顆粒團的流體動力學半徑Rh=2.17μm，據此計算顆粒團的運動。結果顯示，數值模擬得到的顆粒團運動軌跡與實驗數據吻合良好，都反映出由于水平方向聲波對顆粒的黏性夾帶作用，顆粒團的運動軌跡呈現出“S”形；駐波聲場作用使得顆粒團在水平方向發生往復振動的同時，振動中心向波節點漂移，這種定向漂移效應是由駐波聲場中氣體介質的非對稱運動引起的[9,18,31]。這表明本文的數學模型和計算方法能夠合理預測顆粒團在駐波聲場作用下的運動特性。

表1 數值模擬參數Table 1 The numerical simulation parameters

圖2 實驗與數值模擬得到的顆粒團運動軌跡的對比Fig.2 Comparison of the motion trajectories of particle agglomerate obtained by experiment and numerical simulation

2.2 顆粒團的運動特性

由式(1)、式(2)可知，顆粒團的分形維數和流體動力學半徑由組成顆粒團的原生顆粒半徑、數目和排列情況所決定，而流體動力學半徑直接影響到顆粒團的動力學特性。鑒于駐波聲場中的顆粒團在聲波波動方向上存在邊振動、邊漂移的運動特性，同時由于慣性顆粒團的振動相位滯后于聲波引起的氣體介質振動相位，下文利用夾帶系數和相位滯后表示顆粒團的振動特性，利用漂移系數表示顆粒團的漂移特性。其中，夾帶系數定義為顆粒團速度振幅與顆粒團振動中心位置氣體介質速度振幅之比，相位滯后定義為顆粒團振動相位滯后于氣體介質的量，漂移系數定義為顆粒團振動中心的漂移速度與顆粒團振動中心位置氣體介質速度振幅之比。本部分數值模擬采用的氣相溫度T、氣相靜壓p、顆粒材料密度ρ、時間步長?t同表1一致，采用的聲場為L=150 dB、f=2000 Hz的駐波聲場，顆粒初始位置為x0=7λ/8。

2.2.1 原生顆粒半徑對顆粒團運動特性的影響

圖3 原生顆粒半徑對顆粒團運動特性的影響Fig.3 effect of the radius of primary particles on the motion characteristics of particle agglomerate

圖3給出了由半徑為a1和a2的兩個原生顆粒組成的顆粒團的夾帶系數、相位滯后和漂移系數，及與等體積球形顆粒的對比關系。等體積球形顆粒的半徑為圖中顆粒團及相應的等體積球形顆粒的曲線顏色相同。可見，當a1<0.1μm或a1>5μm時，顆粒團與等體積球形顆粒的夾帶系數、相位滯后和漂移系數的差異很小；當0.1μm 6a16 5μm時，顆粒團的夾帶系數低于等體積球形顆粒，相位滯后高于等體積球形顆粒；a2較小時，顆粒團的漂移系數低于等體積球形顆粒，a2較大時，顆粒團的漂移系數大于等體積球形顆粒。這是因為當a1<0.1μm或a1>5μm 時，兩原生顆粒的半徑相差較大，顆粒團的運動特性由粒徑較大的原生顆粒決定；同時，等體積球形顆粒的半徑也由較大的原生顆粒所主導，使得顆粒團的運動特性與等體積球形顆粒差異很小。然而，當兩原生顆粒的半徑相差不多時，顆粒團和等體積球形顆粒運動特性的差異變得明顯，此時，若將顆粒團當作球形顆粒進行處理，將帶來較大的計算誤差。

由圖3還可以看出，當原生顆粒2的半徑a2為0.25μm、0.5μm和1μm時，隨著原生顆粒1的半徑a1的增加，顆粒團的漂移系數先增加后減小，存在一個漂移系數峰值；當a2=2.5μm時，隨著a1的增加，顆粒團的漂移系數趨于減小，這是由顆粒團的慣性決定的。當a2較小(a2<1μm)時，a1由0.01μm增加至10μm的過程中，顆粒團的質量增加，顆粒團由零慣性顆粒向有限慣性顆粒轉變[32]；顆粒團處于零慣性區域時，其能夠被聲波充分夾帶，運動特性表現為隨聲波往復振動，此時夾帶系數接近1，相位滯后接近0，漂移系數接近0；隨著慣性增加，顆粒團處于有限慣性區域，夾帶系數單調減小，相位滯后單調增加，起初顆粒運動過程中經歷的流場非對稱性增強，顆粒團漂移系數增加，但是當顆粒團慣性增加到一定程度時，顆粒團將難以被聲波夾帶，特別是在極限情況下顆粒團既不發生振動又不發生漂移，漂移系數為0，因此存在一個漂移系數峰值。然而，當a2較大(如a2=2.5μm)時，顆粒團慣性始終保持在較高的水平，顆粒團漂移系數由0增大的階段將不再存在，因此顆粒團漂移系數呈現出單調減小的特征。

2.2.2 原生顆粒數目對顆粒團運動特性的影響

為了探究原生顆粒數目對駐波聲場中顆粒團運動特性的影響，在顆粒團分形維數Df=1.85、原生顆粒半徑a=0.25μm條件下，對原生顆粒數目N=2～100時顆粒團的運動特性進行數值模擬，得到了顆粒團及其等體積球形顆粒的夾帶系數、相位滯后和漂移系數隨原生顆粒數目的變化關系，如圖4所示。結果表明，隨著原生顆粒數目的增加，顆粒團的慣性趨于增大，導致其夾帶系數減小，相位滯后增加，漂移系數先增大后減小。此外，當原生顆粒數目較少時，顆粒團及相應的等體積球形顆粒的質量均較小，二者都接近于零慣性區域，對聲波有很強的跟隨性，顆粒團與等體積球形顆粒的運動特性差異較?。坏S著原生顆粒數目的增加，顆粒團及等體積球形顆粒慣性增加，由于顆粒團流體動力學半徑和等體積球形顆粒半徑的差異帶來的運動特性的差異變得顯著。需要指出的是，圖4(c)中，顆粒團和等體積球形顆粒的漂移系數曲線存在一個交點，交點對應的原生顆粒數目為72，這表明此時顆粒團和等體積球形顆粒漂移系數相等。究其原因是，漂移系數隨流體動力學半徑的增加先增加后減小[31]，由于非球形顆粒團的流體動力學半徑大于等體積球形顆粒的流體動力學半徑，等體積球形顆粒漂移系數的變化滯后于非球形顆粒團，進而引起顆粒團和等體積球形顆粒的漂移系數曲線相交。圖4中結果還表明，在聲凝并過程中，隨著聲波作用時間的延長，形成了更多顆粒團，并且顆粒團的粒度也更大，考慮到顆粒團和等體積球形顆粒在運動特性上的顯著差異，在對聲凝并進行建模時若將顆粒團視為球形，將帶來較大的誤差，因此有必要在聲凝并模型中充分考慮顆粒團的分形結構對其動力學行為的影響。

圖4 原生顆粒數目對顆粒團運動特性的影響Fig.4 effect of the number of primary particles on the motion characteristics of particle agglomerate

2.2.3 原生顆粒排列情況對顆粒團運動特性的影響

分形維數Df的大小直接反映了顆粒團中原生顆粒排列的致密程度。Df越小，原生顆粒排列越松散；Df越大，原生顆粒排列越致密。為了探討原生顆粒排列情況對顆粒團運動特性的影響，對由5個相同半徑的原生顆粒組成的顆粒團在駐波聲場中的運動特性進行數值預測。圖5給出了顆粒團中原生顆粒的排列示意圖。圖中，θ為旋轉角，即以直鏈形顆粒團為基準，右邊3個顆粒繞基準顆粒團的質心旋轉，與基準位置的夾角。對于圖5所示的顆粒團結構，在原生顆粒半徑a和旋轉角θ給定的條件下，根據幾何關系確定顆粒團的半徑R，聯立式(1)、式(4)，求解獲得顆粒團的分形維數Df和流體動力學半徑Rh。

圖5 顆粒團中原生顆粒排列情況示意圖Fig.5 Schematic diagram of the packing structure of primary particles in particle agglomerate

圖6給出了數值模擬得到的顆粒團的夾帶系數、相位滯后和漂移系數隨分形維數的變化關系。在θ由0?增加至120?的過程中，原生顆粒的排列變得更加致密，顆粒團的分形維數增加，由式(1)可見顆粒團的流體動力學半徑減小，顆粒團更容易被聲波夾帶，導致夾帶系數增加，同時相位滯后減小。由于等體積球形顆?？梢暈樵w粒的最致密排列形成的顆粒團(對應于Df=3)，等體積球形顆粒的運動特性不隨顆粒團分形維數的變化而變化，且其夾帶系數比實際顆粒團高，而相位滯后比實際顆粒團低，等體積球形顆粒與線形顆粒團振動特性差異最大。圖6的結果還表明，顆粒團與等體積球顆粒漂移系數的差異隨分形維數的增加而減小，這是因為分形維數越大，顆粒團的結構越致密，顆粒團流體動力學半徑與等體積球顆粒粒徑差異越小，導致漂移速度差異縮小。此外，在原生顆粒粒徑較小時(a6 0.5μm時)，漂移系數隨分形維數的增加而減?。辉谠w粒粒徑較大時(a>1.0μm時)，漂移系數隨分形維數的增加而增大。其原因是原生顆粒粒徑的增加導致顆粒團的流體動力學半徑增加，當顆粒團的流體動力半徑達到一定值時，顆粒團的定向漂移運動受到抑制。

圖6 原生顆粒排列情況對顆粒團運動特性的影響Fig.6 effect of the packing structure of primary particles on the motion characteristics of particle agglomerate

3 結論

基于分形理論的質量-半徑關系式和顆粒團無量綱流體動力學半徑分形維數關系式，建立外加聲場作用下球形原生顆粒組成的顆粒團的動力學模型，利用四階變步長Runge-Kutta算法和二階Adams插值算法對顆粒團運動方程進行求解，將模擬結果與實驗結果進行對比，驗證了數值模擬結果的正確性。在此基礎上，對顆粒團在駐波聲場中的夾帶系數、相位滯后和漂移系數隨原生顆粒半徑、數目和排列情況的變化特性進行數值預測，并與等體積球形顆粒進行比較。通過本文研究，得出以下結論：

(1)對于兩個原生顆粒組成的顆粒團，兩個原生顆粒的半徑相差較大，顆粒團與等體積球形顆粒的運動特性差異很?。划攦稍w粒的半徑相差不多時，顆粒團和等體積球形顆粒運動特性的差異變得明顯。

(2)分形維數一定時，隨著原生顆粒數目的增多，顆粒團的夾帶系數減小，相位滯后增加，漂移系數先增大后減小，顆粒團與等體積球形顆粒的動力學行為存在顯著差異。

(3)原生顆粒的排列情況決定了顆粒團的分形維數，原生顆粒的排列趨于致密時，分形維數增加，引起顆粒團的夾帶系數增加，相位滯后減小，漂移系數發生單調變化，顆粒團與等體積球形顆粒的運動特性的差異趨于減小。