基于蒙特卡羅原理的混合顆粒三相體系聲衰減計算模型研究*

趙寧寧 肖新宇 凡鳳仙 蘇明旭

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

從單個固體和液滴顆粒的聲吸收和散射特性計算入手,基于概率統計的蒙特卡羅方法(MCM),將聲波以離散化的聲子加以處理,通過追蹤其運動歷程并進行事件統計,建立一種氣體介質中球形混合顆粒的聲衰減預測模型.對空氣中鋁粉顆粒和亞微米級水滴的聲衰減分別計算和驗證后,將模型推廣至含有混合顆粒的三相體系,對鋁粉和液滴構成的單、多分散混合顆粒體系進行數值研究.結果表明:兩類顆粒的聲吸收特性差異明顯,其散射聲壓均隨顆粒無因次尺寸kR 的增加呈現從后向散射占主要地位逐漸過渡到前向增強的趨勢.氣液固混合顆粒三相體系中,顆粒類型對于聲衰減影響明顯、且隨濃度的增加不同顆粒的衰減貢獻不再遵循隨混合比的線性遞變關系;對于多分散體系而言,聲衰減譜受平均粒徑影響較大,對于粒徑分布寬度參數則不敏感.模型可進一步結合數學反演形成混合顆粒體系測量的理論基礎.

1 引言

氣液固共存的顆粒三相體系廣泛存在于石油化工、煤炭液化、生物化工、環境工程等領域,如天然氣在輸送過程中管道內氣(天然氣)液(液態水)固(水合物顆粒)[1]、電力行業脫硫塔內氣(空氣)液(霧)固(煙塵)[2]、化工行業氣液固三相反應器內[3,4]、環境問題中氣(空氣)液(霧)固(霾)[5]等三相體系.氣液固三相顆粒體系中顆粒粒徑、相含率(各相占比)和濃度的準確監測對生產和研究均具有重要意義.圖像法[3]、層析成像法[4]、光散射法[5]等測量方法已得到學者關注,而超聲衰減法具備穿透性好、適應性強、非接觸式測量等優點,其中三相體系的超聲傳播規律、聲衰減預測及物理模型的建立尤為重要,需要開展專門研究.

超聲在顆粒體系中傳播時,會引起聲學基本物理參數的變化,如聲壓、聲速、聲衰減和聲阻抗等.針對兩相流顆粒體系中聲衰減特性的研究,Epstein 與Carhart[6]和Allegra 與Hawley[7]建立ECAH(Epstein-Carhart-Allegra-Hawley)模型,同時考慮液體介質的黏性、顆粒的彈性效應、熱傳導等作用,根據質量、動量和能量守恒定律,通過求解球坐標下的波動方程得到散射系數,奠定了超聲波衰減理論模型的基礎.之后,McClements[8-10]提出了“長波長”條件下聲衰減與聲速計算理論模型(McC模型),引入黏性厚度和熱厚度兩個參數,簡化聲衰減模型,研究了超聲吸收衰減效應和長波長區的聲衰減特性,并進行了實驗驗證.近年來,Parker 等[11]與Liu[12]研究了納米顆粒懸濁液的聲衰減,對顆粒體系進行測量和表征.董黎麗等[13]建立了乳濁液的聲衰減反演模型,研究了多分散脂肪乳濁液的粒徑分布測量問題.Wang 等[14,15]建立了耦合相修正模型,計算了空氣介質中固體顆粒的聲衰減.杜娜與蘇明旭[16]建立了水中氣泡的聲散射模型,研究了有黏條件下氣泡的聲學特性.侯森等[17]建立了聲衰減反演氣泡分布模型.陳時等[18]建立了含混合氣泡液體中的聲傳播模型.郭盼盼等[19]、李運思等[20]和Huang 等[21]建立了蒙特卡羅原理的聲衰減計算模型,應用于液固兩相、液液兩相單分散和多分散體系的聲衰減譜預測.

對兩相體系中顆粒的聲衰減特性研究已經趨于成熟,但是對于氣體介質中液、固混合顆粒構成的三相體系,其聲傳播規律和聲衰減計算仍屬難題.此外,實際的混合顆粒系粒徑分布并非均一,而是呈現多分散分布,ECAH 模型等均不能直接應用.因此,作者通過引入蒙特卡羅方法,利用其數理統計的特點,研究一種適用于三相、混合多分散顆粒體系的超聲波衰減模型.

2 模 型

蒙特卡羅方法是一種以概率統計理論為基礎的數值方法,可將其引入到顆粒兩相及多相流的模擬計算[19-23].類比于光散射計算中的“光子”概念,其核心思想是將入射聲波按照“聲子”概念抽象并離散化處理,即將換能器發出的連續超聲波離散成大量時序上相互獨立的聲子.

如圖1 所示,氣液固三相體系中混合了兩類球形顆粒,藍色為液滴顆粒,紅色為鋁粉顆粒,介質為空氣.結構參數L為發射器(T)和接收器(R)之間距離,d為接收器的直徑,H為上下邊界距離,ln為第n次(n=1,2,3···)隨機散射后聲子運動自由程.發射器發射聲波能量表示成大量非連續的聲子并對其進行追蹤.聲子行為包括透射、吸收、越界或者復散射后被接收器接收.通過統計接收器獲取聲子數,即可模擬在一定的顆粒粒徑、體積濃度和超聲頻率條件混合顆粒系中聲波行為并計算聲衰減系數.

圖1 聲子在混合顆粒體系中傳播過程Fig.1.Schematic of a phonon propagation through a mixed particle system.

鑒于本文涉及混合顆粒的氣液固三相體系,當聲子和顆粒發生碰撞時,須判斷顆粒類型,即鋁粉或液滴:

式中,ε1是在[0,1]區間內服從均勻分布的隨機數;混合比φ為液滴在整個混合顆粒體系中所占的數目比,例如φ取0 時顆粒全為鋁粉,1 時則全為液滴.

下一步,判斷聲子碰撞顆粒后可能發生事件,定義聲散射系數和消聲系數比值P:

其中消聲系數QTQa+Qs.Qa和Qs分別為聲吸收系數和聲散射系數,由Hay和Mercer[24]方法算得.接下來,根據設定條件判斷聲子的下一事件:

式中,ε2是[0,1]區間服從均勻分布的隨機數;n為散射次數;x為聲子沿著x方向運動的路程;z為聲子沿著z方向運動的路程.通常進入氣液固三相體系聲子可能被顆粒吸收、被接收器直接接收(透射)、越界逃逸或復散射.前三種情況下,聲子歷程終止,而對于復散射,則需進一步追蹤聲子.

在空氣中傳播的聲子遇到顆粒發生聲散射,需確定聲子散射方向及聲子兩次散射之間的自由程.采用(4)式計算歸一化散射聲壓f(θ):

式中,θ為散射角,取值范圍為0到360°.p(θ)是顆粒表面散射聲壓分布函數(顆粒為球形時與方位角無關),對于彈性固體顆粒,可由Faran 理論計算[25],對于液滴顆粒,采用液體球的散射聲壓計算方法[26].(5)式和(6)式分別給出其核心計算式:

式中,jn和nn分別是第一類球Bessel 函數和第二類球Bessel 函數,k為入射聲波波數;r為接收點距離,取顆粒半徑的100 倍;Pn(cosθ)是勒讓德多項式,散射系數An由Faran理論算出.

為確定聲子出射方向,將散射角θ從0到360°劃分為360 份,引入另一隨機數ε3與歸一化散射聲壓分布f(θ)比較,如果:

則聲子散射后的出射角就為θM,M取值范圍為1—360.之后確定隨機散射聲子的自由程l:

ε4是在[0,1]區間服從均勻分布的另一隨機數.

重復(2)式—(8)式的聲子運動過程,在確定經過n+1次散射后聲子仍在顆粒體系中后,可得其散射坐標:

式中,xn和zn是聲子第n次散射后的位置.至此,所有聲子都要經歷上述過程,根據(3)式條件判斷并統計其最終去向,模型中聲子運動全過程可由如下流程圖所示,可得聲衰減系數計算公式為

式中,Ndet是探測接收器接收聲子數目;Ntot是聲子樣本容量.蒙特卡羅模型的算法流程圖如圖2 所示.

圖2 蒙特卡羅模型的算法流程圖Fig.2.Flow chart of Monte Carlo model.

3 計算過程和結果驗證

3.1 單顆粒聲散射

計算對象鋁粉和液(水)滴顆粒的物性參數見表1,為比較不同類型顆粒的聲散射和聲吸收特性,定義系數Ps為液滴與鋁粉顆粒聲散射系數Qs之比,Pa為液滴與鋁粉顆粒聲吸收系數Qa之比.在超聲頻率為20 kHz、體積濃度為0.01%、顆粒半徑為0.1—1 μm 時,計算其隨顆粒半徑變化情況.由圖3可知,兩類顆粒的聲散射特性較接近,而吸收特性差異明顯并隨顆粒半徑變化,同粒徑下,鋁粉顆粒聲吸收更強,對應聲衰減也會更大.

圖3 液滴和鋁粉的聲散射及吸收系數比Fig.3.Ratio of scattering and absorption coefficients of droplet and aluminum particle.

表1 數值計算中介質和顆粒物性參數(25 ℃)Table 1.Parameters in the numerical simulation (25 ℃).

由于低頻條件下空氣中固體和液滴顆粒散射特性差異不大,圖4 僅給出液滴的歸一化散射聲壓((6)式計算).定義平面聲波傳播方向即圖中0°—90°和270°—360°為前向,反之即為后向.當無因次參量kR=0.5和kR=1 時,顆粒后向散射比較均勻且占主要地位;此后隨著無因次參量的增大,散射旁瓣增多,前向散射效應逐漸增強、向前集中,且前向主瓣散射角變小.如前所述,單顆粒的吸收和散射結果將直接耦合到模型中用于確定聲子的吸收概率和出射方向.

圖4 液滴顆粒散射聲壓Fig.4.Scattering pressure of droplet.

3.2 聲子統計

采用蒙特卡羅方法計算聲衰減系數時,其計算精度受限于聲子的樣本容量,樣本容量越高,精度越高,圖5(a)展示了液滴顆粒體系在聲波頻率為50 kHz,顆粒半徑為0.1 μm,體積濃度為0.01%時,分別采用1×104,1×105和1×106樣本容量,計算100次統計衰減系數平均值和相對偏差.可知采用聲子1×106樣本容量,數據相對偏差在—0.59%至0.57%,此時樣本容量足夠大,模擬結果受隨機因素影響的波動較小.

圖5 不同條件聲子統計結果 (a) 不同聲子數目時相對偏差;(b) 不同體積濃度時聲子去向Fig.5.Statistic of phonons under different conditions:(a)Relative deviations at different number of phonons;(b)phonon events at different particle volume concentrations.

為研究氣體介質中聲波傳播物理過程和聲子最終去向,對于鋁粉、液滴顆粒體系,設定聲波頻率為20 kHz,顆粒半徑為0.1 μm,體積濃度為0.01%—0.1%,聲子樣本容量為1×106.分別統計了聲子發生復散射、吸收、越界以及透射去向的數量.

統計結果如圖5(b)所示,聲子主要歷經吸收和直接透射過程.隨體積濃度增加,聲子在運動過程中與顆粒碰撞概率加大,碰撞后的散射聲子易偏離接收范圍,表現為透射聲子數減少,而越界的聲子數遞增.同時,經歷復散射過程的聲子數逐漸增加,但其數目總量不大(小于 4×104),小于越界逃逸和被吸收聲子數,因此并不占主導地位.對于兩種類型顆粒,被吸收聲子數差異明顯,這由不同類型顆粒自身吸收差異決定(如3.1 節),從而對聲衰減影響亦不同.

3.3 模型驗證

為驗證本文發展蒙特卡羅模型和聲衰減計算程序準確性,首先對單一類型顆粒兩相體系,計算了不同粒徑下聲衰減系數數值,將結果和經典ECAH模型以及McC 模型預測結果進行對比(對比模型結果均通過文獻公式編程運算而得).圖6(a)給出液滴、鋁粉顆粒在體積濃度Cv=0.01%,超聲頻率f=50 kHz 時采用蒙特卡羅方法、ECAH 模型、McC 模型三種結果對比.由圖6(a)可知,三種模型的計算結果較為吻合.在體積濃度和頻率一定條件下,受顆粒的黏性耗散和熱耗散影響,聲衰減系數隨著顆粒半徑增大呈現先增后減變化趨勢,該聲衰減極大值位于0.2—0.4 μm 之間.計算條件為顆粒半徑0.1 μm、體積濃度0.01%、頻率50 kHz 時,各物性參數分別單獨增加10%條件下,聲衰減系數對顆粒物性參數的敏感性如下圖6(b).由于顆粒相和連續相之間的高密度差,顆粒的黏性耗散和熱耗散效應對聲衰減結果的影響最大,其相關物性參數主要為顆粒密度(鋁粉19.12%、液滴19.54%)、比熱容(鋁粉4.35%、液滴17.30%)、導熱系數(鋁粉1.39%)、聲吸收系數(液滴0.42%).

圖6 模型驗證及物性敏感性 (a) 聲衰減系數隨顆粒半徑變化;(b) 聲衰減系數隨顆粒物性參數變化Fig.6.Model validation and the sensitivity of physical parameters:(a) Ultrasonic attenuation coefficient varies with the particle radius;(b) ultrasonic attenuation coefficient varies with the particle parameters.

將蒙特卡羅方法計算結果分別與鋁粉顆粒和液滴顆粒的實驗結果對比.可以看出,無論圖7(a)所示鋁粉顆粒在半徑R=2.8 μm,質量濃度Cm=0.032 kg/kg 條件下聲衰減系數[14],還是圖7(b)所示液滴顆粒在超聲頻率f=40,200 kHz,體積濃度Cv=0.0043%,0.0114%條件下聲衰減系數[27],同等條件下作者采用蒙特卡羅方法的數值計算結果與實驗數據整體吻合.

圖7 實驗結果對比 (a) 鋁粉顆粒;(b) 液滴顆粒 (R1和R2 分別為文獻中超聲和圖像法測得液滴半徑)Fig.7.Comparison with experimental results:(a) Aluminum particle;(b) droplet (R1 and R2 are droplet radius measured by ultrasonic and image method in the reference respectively).

4 混合顆粒系結果

4.1 單分散系

對單一類型固體顆粒、液滴顆粒的計算和驗證表明蒙特卡羅模型從建模思路到編程可靠,但在實際測量對象中,常包含多種類型顆粒,因此進一步將其拓展至鋁粉-液滴混合顆粒三相體系.

對于單分散混合顆粒三相體系,圖8 所示為頻率f=40 kHz,顆粒半徑R=0.1 μm,不同體積濃度條件下聲衰減系數隨顆粒數目混合比的變化曲線.由圖8 可知,隨著混合比φ增大即液滴所占比例增大(φ為1 時全為液滴),混合顆粒體系的聲衰減系數逐漸減小,說明液滴顆粒對聲衰減的貢獻較鋁粉顆粒小.在體積濃度較低時(如Cv=0.01%),聲衰減系數隨混合比近似呈現線性變化,而隨著體積濃度的進一步增加,尤其到0.05%,出現非線性變化趨勢,由于體積濃度增加,聲子碰撞概率增加,鋁粉顆粒的聲吸收作用進一步強化.值得注意的是,根據圖示聲衰減譜,在已知體積濃度和顆粒粒度前提下,可利用反演問題求解顆粒體系的混合顆粒數目比,從而獲知不同類型顆粒混合時的占比信息.

圖8 單分散三相混合體系聲衰減Fig.8.Ultrasonic attenuation coefficient in monodisperse mixed system.

4.2 多分散系

實際的顆粒系中粒徑往往并非均一(即非理想的單分散系),對于呈多分散分布混合顆粒三相體系,假設顆粒粒徑分布服從正態分布函數,即:

其中N為顆粒數目;R為特征尺寸參數(函數期望值);σ為分布寬度參數(函數標準差).

設定混合顆粒系的粒徑分布為雙峰模態,兩類顆粒占比均為0.5,鋁粉顆粒特征尺寸參數R1為0.10 μm,分布寬度參數σ1為0.015 μm;液滴顆粒特征尺寸參數R2為0.15 μm,分布寬度參數σ2分別為0.010,0.015和0.020 μm,另外設定一組R1為0.10 μm,R2為0.15 μm,σ1=σ2=0 的無分布寬度混合顆粒組對比,顆粒混合后的粒徑分布如圖9 所示.

圖9 混合顆粒體系粒徑分布Fig.9.Particle size distribution of mixed particle system.

圖10(a)為體積濃度Cv=0.02%,四種不同粒徑分布組合下聲衰減譜曲線,可知在頻率低于40 kHz 時,液滴粒徑分布參數的變化對聲衰減譜影響較小,只有在足夠頻率帶寬(如至100 kHz)的譜,能較好地同步分析平均粒徑和分布信息.圖10(b)為頻率f=40 kHz,聲衰減系數隨體積濃度增加曲線呈現非線性變化趨勢,結合圖5(b)中復散射聲子討論,說明顆粒體系中發生復散射的概率增大,進而引起衰減系數的增長趨勢變化.此外,與分布寬度為零的單分散混合顆粒的結果對比,正態分布的多分散體系由于顆粒尺寸落在衰減峰值區外,導致分布寬度參數越大,聲衰減數值越小.同樣,當兩種顆粒分布有重疊或不同比例時,超聲衰減譜同樣具備差異性,可以有效區分不同混合比時的聲衰減結果.

圖10 多分散三相混合體系聲衰減 (a) 隨頻率變化;(b) 隨濃度變化Fig.10.Ultrasonic attenuation coefficient in polydisperse mixed system:(a) Curves with frequency;(b) curves with volume concentration.

5 結論

通過蒙特卡羅方法建立了預測氣、液、固混合顆粒三相體系聲衰減的計算模型,研究混合顆粒體系中粒徑分布、頻率、體積濃度、混合顆粒數目比等多因素對聲衰減的影響,得出以下結論:

1) 發展了一種基于蒙特卡羅原理的聲衰減建模方法,較為準確地預測單一顆粒體系如液滴、鋁粉顆粒的聲衰減系數,其結果和ECAH 模型、McC模型較吻合,且與實驗數據整體吻合.

2) 該方法擴展到了混合顆粒三相體系中聲衰減的預測.計算結果表明,由于不同類型顆粒的聲吸收特性不同,在體積濃度較低條件下,隨著混合顆粒數目比的增加,聲衰減呈線性變化,而當體積濃度較高時,聲衰減呈現非線性變化趨勢.考慮顆粒本身物性影響,相同粒徑條件下鋁粉顆粒比液滴對聲衰減的貢獻更大.

3) 對于多分散混合顆粒體系,聲衰減結果對特征尺寸和分布寬度的敏感性不同.由于體積濃度增加時的復散射效應,超聲頻率和體積濃度改變對衰減結果有不同影響,需要考慮不同類型顆粒物性及顆粒粒徑分布的綜合作用.

由于本文僅針對了低濃度條件下的球形顆粒系進行建模和聲衰減預測,后續有望繼續進一步拓展至高濃度條件下的混合顆粒體系以及反演問題研究.此外,可以將模型拓展至非球形顆粒的研究,例如橢球形固體顆粒、液滴顆粒等.