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

(上海理工大學(xué)能源與動(dòng)力工程學(xué)院,上海 200093)

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

1 引言

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

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

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

2 模 型

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

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

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

鑒于本文涉及混合顆粒的氣液固三相體系,當(dāng)聲子和顆粒發(fā)生碰撞時(shí),須判斷顆粒類(lèi)型,即鋁粉或液滴:

式中,ε1是在[0,1]區(qū)間內(nèi)服從均勻分布的隨機(jī)數(shù);混合比φ為液滴在整個(gè)混合顆粒體系中所占的數(shù)目比,例如φ取0 時(shí)顆粒全為鋁粉,1 時(shí)則全為液滴.

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

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

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

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

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

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

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

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

ε4是在[0,1]區(qū)間服從均勻分布的另一隨機(jī)數(shù).

重復(fù)(2)式—(8)式的聲子運(yùn)動(dòng)過(guò)程,在確定經(jīng)過(guò)n+1次散射后聲子仍在顆粒體系中后,可得其散射坐標(biāo):

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

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

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

3 計(jì)算過(guò)程和結(jié)果驗(yàn)證

3.1 單顆粒聲散射

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

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

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

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

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

3.2 聲子統(tǒng)計(jì)

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

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

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

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

3.3 模型驗(yàn)證

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

圖6 模型驗(yàn)證及物性敏感性 (a) 聲衰減系數(shù)隨顆粒半徑變化;(b) 聲衰減系數(shù)隨顆粒物性參數(shù)變化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.

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

圖7 實(shí)驗(yàn)結(jié)果對(duì)比 (a) 鋁粉顆粒;(b) 液滴顆粒 (R1和R2 分別為文獻(xiàn)中超聲和圖像法測(cè)得液滴半徑)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 混合顆粒系結(jié)果

4.1 單分散系

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

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

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

4.2 多分散系

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

其中N為顆粒數(shù)目;R為特征尺寸參數(shù)(函數(shù)期望值);σ為分布寬度參數(shù)(函數(shù)標(biāo)準(zhǔn)差).

設(shè)定混合顆粒系的粒徑分布為雙峰模態(tài),兩類(lèi)顆粒占比均為0.5,鋁粉顆粒特征尺寸參數(shù)R1為0.10 μm,分布寬度參數(shù)σ1為0.015 μm;液滴顆粒特征尺寸參數(shù)R2為0.15 μm,分布寬度參數(shù)σ2分別為0.010,0.015和0.020 μm,另外設(shè)定一組R1為0.10 μm,R2為0.15 μm,σ1=σ2=0 的無(wú)分布寬度混合顆粒組對(duì)比,顆?；旌虾蟮牧椒植既鐖D9 所示.

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

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

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

5 結(jié)論

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

1) 發(fā)展了一種基于蒙特卡羅原理的聲衰減建模方法,較為準(zhǔn)確地預(yù)測(cè)單一顆粒體系如液滴、鋁粉顆粒的聲衰減系數(shù),其結(jié)果和ECAH 模型、McC模型較吻合,且與實(shí)驗(yàn)數(shù)據(jù)整體吻合.

2) 該方法擴(kuò)展到了混合顆粒三相體系中聲衰減的預(yù)測(cè).計(jì)算結(jié)果表明,由于不同類(lèi)型顆粒的聲吸收特性不同,在體積濃度較低條件下,隨著混合顆粒數(shù)目比的增加,聲衰減呈線性變化,而當(dāng)體積濃度較高時(shí),聲衰減呈現(xiàn)非線性變化趨勢(shì).考慮顆粒本身物性影響,相同粒徑條件下鋁粉顆粒比液滴對(duì)聲衰減的貢獻(xiàn)更大.

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

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