祖曉萌， 朱曙光

(南京理工大學(xué)能源與動(dòng)力工程學(xué)院， 江蘇南京210018)

聲波能穿透不透光介質(zhì)，使得超聲測(cè)量技術(shù)能夠用于高含量氣固兩相系統(tǒng)的特性研究。理論模型將超聲衰減和聲速變化與兩相系統(tǒng)的特征參數(shù)相關(guān)聯(lián)，聲速測(cè)量更為簡(jiǎn)單和準(zhǔn)確，且對(duì)顆粒分布的不均勻性更為敏感[1]；但由于聲速的變化還依賴于兩相介質(zhì)的物性參數(shù)如黏度、密度、顆粒材料等[2]，所以理論模型更多集中于超聲衰減的研究。

在實(shí)際應(yīng)用中，即使在理想的測(cè)量條件下，測(cè)得的衰減值與理論預(yù)測(cè)值的一致性也僅能達(dá)到30%左右[3]。 高含量?jī)上嘞到y(tǒng)聲衰減理論模型的建立難點(diǎn)在于多次散射現(xiàn)象和顆粒間相互作用所導(dǎo)致的聲衰減與所求介質(zhì)參數(shù)的非線性對(duì)應(yīng)關(guān)系，且對(duì)于濃度越高、 顆粒越小時(shí)，多次散射和顆粒相互作用所發(fā)生的可能形式則越多，即從理論層面對(duì)兩相介質(zhì)間的分布和與聲波相互作用的隨機(jī)規(guī)律進(jìn)行分析相對(duì)困難[4]。 蒙特卡羅方法[5-7]因能夠?qū)?fù)雜抽象的物理過程簡(jiǎn)單具體化而得到快速發(fā)展， 并被成功地應(yīng)用于預(yù)測(cè)高含量液固兩相系統(tǒng)的聲衰減。本文中利用蒙特卡羅方法對(duì)超聲波透射、散射和衰減的發(fā)生概率和作用方式進(jìn)行建模，來模擬高含量氣固兩相系統(tǒng)聲衰減。

1 蒙特卡羅模型建立

1.1 模型的實(shí)現(xiàn)過程

通過將入射聲波按照“聲子”概念離散化，每個(gè)聲子的初始聲能量為1。聲子傳播會(huì)發(fā)生多次的透射和散射，聲子透射時(shí)認(rèn)為無能量衰減，當(dāng)聲子發(fā)生散射時(shí)，存在黏性衰減和熱衰減，并將兩者的作用效果用單次散射能量損失Eloss來表示。其模型示意圖如圖1所示,其中H為發(fā)射和接收間距離，D為上下邊界距離。

1.2 透射概率

兩相系統(tǒng)的模型是邊長(zhǎng)為H=0.02 m的正方體，通過橫截面的縱向顆粒數(shù)即可計(jì)算顆粒體積分?jǐn)?shù)，

圖1 蒙特卡羅模型示意圖Fig.1 Diagrammatic sketch of Monte Carlo model

(1)

式中：φi為不同粒徑顆粒體積分?jǐn)?shù)；Ri為顆粒半徑；rnum為橫截面的縱向顆粒數(shù)；Vs為顆粒所占總體積；Va為兩相系統(tǒng)的總體積。

將后續(xù)實(shí)驗(yàn)中不同粒徑顆粒進(jìn)行單層堆積時(shí)的體積分?jǐn)?shù)進(jìn)行測(cè)量，體積分?jǐn)?shù)的實(shí)測(cè)值見表1，具體計(jì)算如式(2)。

(2)

式中：Vr為顆粒的實(shí)際體積；Vb為顆粒表觀體積；m為總的顆粒質(zhì)量；ρr為顆粒密度；ρb為顆粒表觀密度。根據(jù)ρr和ρb便可得到顆粒堆積層顆粒體積分?jǐn)?shù)。

表觀密度的測(cè)量是將相同粒徑的玻璃微珠填滿量筒，并通過多次測(cè)量整體質(zhì)量取平均即可求出該粒徑顆粒的表觀密度[8]，顆粒材料的真實(shí)密度由廠家提供。

將式(1)進(jìn)行變形得到式(3)，將表1中的體積分?jǐn)?shù)帶入式(3)中，即可得到橫截面的縱向顆粒數(shù)，再根據(jù)式(4)可最終求出聲子的透射概率P，不同粒徑對(duì)應(yīng)的透射概率見表2。

表1 不同顆粒粒徑的體積分?jǐn)?shù)測(cè)量值

表2 不同粒徑顆粒的透射概率

(3)

(4)

式中：Sp為單層顆粒的投影面積；Stot為橫截面的面積；R為顆粒半徑。

εn為0到1的隨機(jī)數(shù)，由εn和P的關(guān)系來判斷聲子是否發(fā)生透射，即

εn

(5)

εn≥P散射。

(6)

1.3 聲子能量損失

對(duì)于高含量氣固兩相系統(tǒng)來說，多次散射過程除了會(huì)造成散射衰減，還會(huì)因多次散射而產(chǎn)生更多的熱衰減和黏性衰減。當(dāng)體積一定時(shí)，顆粒體積分?jǐn)?shù)和顆粒數(shù)量成正比。對(duì)低含量氣固兩相系統(tǒng)來說，聲波傳播過程中，每個(gè)顆粒只發(fā)生單次散射，并產(chǎn)生熱黏性衰減，不考慮顆粒間相互作用，所以總的聲衰減與體積分?jǐn)?shù)成正比。而當(dāng)體積分?jǐn)?shù)不斷增大，會(huì)產(chǎn)生熱、黏性相互作用和多次散射，將導(dǎo)致系統(tǒng)中單個(gè)顆粒多次與聲波相互作用產(chǎn)生更多的熱黏性損失，進(jìn)而使得總的聲衰減與體積分?jǐn)?shù)呈現(xiàn)非線性的對(duì)應(yīng)關(guān)系。同時(shí)對(duì)于高含量系統(tǒng)來說，聲衰減主要取決于黏性衰減和熱衰減，散射衰減則要小得多[9]。

1.3.1 熱衰減和黏性衰減

EACH模型[10-11]較為全面地考慮了熱傳導(dǎo)、黏性和散射對(duì)聲傳播特性的影響，當(dāng)聲波入射到兩相界面，發(fā)生散射的同時(shí)會(huì)產(chǎn)生熱衰減和黏性衰減，并認(rèn)為熱衰減和黏性衰減主要是由于兩相介質(zhì)的可壓縮性和密度差異所導(dǎo)致的。

基于上述理論，模型將由于顆粒自身及其相鄰顆粒間的熱、黏性相互作用產(chǎn)生的聲衰減過程與散射過程相結(jié)合，通過單次散射能量損失Eloss對(duì)其衰減強(qiáng)度進(jìn)行表征。高含量?jī)上嘞到y(tǒng)中，就單個(gè)顆粒來看，小顆粒與相鄰顆粒接觸面積更小，作用強(qiáng)度低，但是相鄰顆粒數(shù)量更多，所以本文中對(duì)不同粒徑顆粒發(fā)生散射進(jìn)行建模時(shí)，其單次散射能量損失Eloss取值相同且為常數(shù)。Eloss中包含了顆粒間相互作用引起的熱衰減和黏性衰減，因此在模型計(jì)算時(shí)，僅對(duì)單個(gè)顆粒進(jìn)行判斷即可，無需考慮其他顆粒對(duì)其的影響。

1.3.2 散射衰減

與熱衰減和黏性衰減的能量損失機(jī)制相比，散射衰減的機(jī)理有所不同。顆粒對(duì)聲波進(jìn)行散射時(shí)不會(huì)產(chǎn)生能量損失，它類似于光散射，只是使得部分聲能量沒有被換能器接收到而產(chǎn)生聲衰減[12]。當(dāng)Eloss取值為0時(shí)，即在聲子的整個(gè)傳播過程均沒有能量損失，只改變傳播路徑，通過程序計(jì)算即可得到散射衰減系數(shù)。

1.4 散射角度

平面波入射到球形顆粒表面，顆粒的散射聲壓分布[13]計(jì)算公式

(7)

將其計(jì)算結(jié)果進(jìn)行均一化后得到函數(shù)B(θ)。將散射聲壓分布角度分為360等份，ζn是服從[0,1]區(qū)間均勻分布的隨機(jī)數(shù)，當(dāng)滿足

(8)

散射角度θ取m+1進(jìn)行散射[14]。

1.5 聲子位置坐標(biāo)

當(dāng)聲子進(jìn)行透射或單次散射后，需要確定過程結(jié)束時(shí)的位置坐標(biāo)xn、yn、zn。單次的透射和散射，其聲子的行進(jìn)路線均被約束于邊長(zhǎng)為h1=H/rnum的正方體內(nèi)。利用角度θi的范圍對(duì)不同情況的聲子位置坐標(biāo)進(jìn)行推導(dǎo)，其示意圖如圖2—3所示。

圖2 單次散射后聲子位置示意圖Fig.2 Diagrammatic sketch of position of phonon圖3 單次透射后聲子位置示意圖Fig.3 Position of phonon after single transmission

1)聲子單次散射過程結(jié)束后位置確定，見式(9)—(12)。

當(dāng)0≤θi≤45 °或315 °<θi≤360 °，

(9)

當(dāng)45 °<θi≤135 °，

(10)

當(dāng)135 °<θi≤255 °，

(11)

當(dāng)225 °<θi≤315 °，

(12)

2)聲子單次透射過程結(jié)束后位置確定如式(13)—(17)所示。

當(dāng)0≤θi≤45 °或315 °≤θi≤360 °，

xn=xn-1+h1,yn=yn-1,zn=zn-1。

(13)

當(dāng)45 °<θi≤135 °，

xn=xn-1,yn=yn-1+h1,zn=zn-1+h1。

(14)

當(dāng)135 °<θi≤225 °，

xn=xn-1-h1,yn=yn-1,zn=zn-1。

(15)

當(dāng)225 °<θi≤315 °，

xn=xn-1,yn=yn-1-h1,zn=zn-1-h1。

(16)

1.6 聲衰減系數(shù)計(jì)算方法

聲子的初始能量M0=1， 變量E的初始值為0， 聲子發(fā)生單次散射時(shí)， 其出射能量與入射能量損失用Eloss表示， 當(dāng)聲子經(jīng)過n次散射后聲子能量為M0(1-Eloss)n。 如果聲子最終被接收， 則將此時(shí)的聲子能量計(jì)入E中， 每接收一個(gè)聲子，E在原有基礎(chǔ)上加上這個(gè)聲子能量。E和聲衰減系數(shù)α按下式計(jì)算。

E=E+M0(1-Eloss)n，

(17)

(18)

式中：E為接收聲子能量總和；Ntotal為發(fā)射的總聲子數(shù)；H為發(fā)射與接收間距離。

1.7 蒙特卡羅程序算法

對(duì)參數(shù)進(jìn)行初始配置，發(fā)射與接收間距離H=0.02 m，上下邊界距離D=0.02 m。發(fā)射的總的聲子數(shù)為Ntotal=100 000，為了對(duì)聲子在傳播過程中的運(yùn)行狀態(tài)進(jìn)行判定，引入了Nin吸收聲子數(shù)和Nout逃逸聲子數(shù)(即沒有被接收到的聲子數(shù))，利用MATLAB進(jìn)行程序編寫，程序框圖如圖4所示。

2 實(shí)驗(yàn)

2.1 裝置

圖5為實(shí)驗(yàn)裝置設(shè)計(jì)圖，工控機(jī)發(fā)出指令控制驅(qū)動(dòng)器驅(qū)動(dòng)發(fā)射換能器發(fā)出超聲波，同時(shí)發(fā)送指令給采集器準(zhǔn)備進(jìn)行信號(hào)采集，聲波穿過顆粒堆積層后被接收換能器接收，通過濾波放大后，經(jīng)采集器送給工控機(jī)進(jìn)行進(jìn)一步數(shù)據(jù)處理。

圖4 程序框圖Fig.4 Diagram of process

2.2 結(jié)果

采用不同粒徑的玻璃微珠緊密堆積來模擬高含量氣固兩相系統(tǒng)， 為了得到此實(shí)驗(yàn)系統(tǒng)的聲衰減系數(shù)， 由公式(19)可知， 在發(fā)射聲壓一定時(shí)， 接收聲壓隨著傳播距離呈指數(shù)形式衰減。 在實(shí)際的測(cè)量工作中是由換能器產(chǎn)生超聲波， 并將聲壓通過接收換能器轉(zhuǎn)換為電壓進(jìn)行聲傳播特性的分析。 實(shí)驗(yàn)中可通過改變堆積顆粒的堆積厚度即改變傳播距離來得到不同傳播距離下的接收電壓值， 并將式(19)中的聲壓用電壓表示， 如式(20)所示， 利用接收電壓的對(duì)數(shù)值與堆積厚度進(jìn)行線性擬合的斜率即為聲衰減系數(shù)。

Pi=Poe-αx，

(19)

Vi=Voe-αx，

(20)

式中：Pi為接收聲壓值；Po為發(fā)射聲壓值；Vi為接收聲壓；Vo為發(fā)射電壓；α為聲衰減系數(shù)；x為聲傳播距離。

利用接收電壓的對(duì)數(shù)值與堆積厚度進(jìn)行線性擬合的斜率即為聲衰減系數(shù)。40 kHz的驅(qū)動(dòng)頻率下不同粒徑玻璃微珠所對(duì)應(yīng)的線性擬合曲線見圖6，衰減系數(shù)值見表3。

圖5 實(shí)驗(yàn)系統(tǒng)設(shè)計(jì)圖Fig.5 Design picture of experimental equipment

3 討論

通過改變程序中單次散射能量損失Eloss的取值， 進(jìn)行程序計(jì)算， 并將計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)值進(jìn)行擬合。Eloss首次嘗試取值為0.01， 即發(fā)生散射時(shí)的聲能量損失為入射到顆粒表面總能量的1%。 運(yùn)行程序后的聲衰減計(jì)算結(jié)果明顯低于實(shí)驗(yàn)值， 當(dāng)進(jìn)一步增大Eloss值， 計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果隨粒徑變化的整體趨勢(shì)相同且逐漸逼近實(shí)驗(yàn)結(jié)果值， 不同的Eloss值所對(duì)應(yīng)的結(jié)果曲線與實(shí)驗(yàn)結(jié)果曲線對(duì)比如圖7所示， 當(dāng)Eloss取5%、 6%時(shí)與實(shí)驗(yàn)結(jié)果曲線擬合程度較好。

將Eloss在5%～6%之間取值來進(jìn)一步提高2條曲線的擬合程度，并利用相關(guān)系數(shù)R2對(duì)擬合程度進(jìn)行表征，數(shù)值越接近1證明擬合程度越高，Eloss取值精確到小數(shù)點(diǎn)后四位。當(dāng)Eloss在5%～6%之間取值時(shí)，所對(duì)應(yīng)R2值變化趨勢(shì)呈先增大后減小，最終得到Eloss=5.5%，R2=0.996 3為最佳結(jié)果，其相應(yīng)的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果曲線如圖8所示。

利用相同思路對(duì)60、 80 kHz的實(shí)驗(yàn)結(jié)果進(jìn)行擬合，結(jié)果分別為Eloss=6.04%，R2=0.995 3;Eloss=7.37%，R2=0.992 6，擬合曲線如圖9和圖10所示。當(dāng)頻率增大，由純空氣所導(dǎo)致的聲衰減變化值相較于熱衰減、黏性衰減和散射衰減來說近乎可以忽略，則因頻率增大所引起的聲衰減值的增大主要是由黏性衰減、熱衰減和散射衰減引起的，單次散射能量損失Eloss的增加表征了熱黏性衰減的增加，而不同頻率對(duì)應(yīng)的散射聲壓分布的變化也最終導(dǎo)致散射衰減的變化。

圖7 不同Eloss值對(duì)應(yīng)聲衰減規(guī)律曲線與實(shí)驗(yàn)結(jié)果對(duì)比Fig.7 Calculated attenuation curves correspouding to different Eloss compare to experimental results圖8 實(shí)驗(yàn)數(shù)據(jù)與模型計(jì)算結(jié)果的曲線最佳擬合Fig.8 The best fitted curve of model calculated compare to experimental curve

圖9 實(shí)驗(yàn)數(shù)據(jù)與模型計(jì)算結(jié)果的曲線最佳擬合Fig.9 The best fitted curve of model calculated compare to experimental curve圖10 實(shí)驗(yàn)數(shù)據(jù)與模型計(jì)算結(jié)果的曲線最佳擬合Fig.10 The best fitted curve of model calculated compare to experimental curve

聲衰減系數(shù)的計(jì)算公式中存在模型尺寸變量H，計(jì)算不同H所對(duì)應(yīng)的聲衰減系數(shù)，其結(jié)果曲線如下圖11所示，可以看出，聲衰減系數(shù)的計(jì)算結(jié)果不會(huì)隨模型尺寸H的改變而變化，即可認(rèn)為聲衰減僅與兩相系統(tǒng)內(nèi)部介質(zhì)參數(shù)有關(guān)。

圖8—11可知，本文中所建立的模型對(duì)于超聲波透射、散射及其能量損失的處理方法符合其在高含量氣固兩相系統(tǒng)中的傳播規(guī)律，能夠反映出高含量氣固兩相系統(tǒng)的衰減特性。

當(dāng)Eloss取值為0時(shí)，通過程序計(jì)算，可得到40 kHz驅(qū)動(dòng)頻率下的散射衰減系數(shù)，將散射衰減隨粒徑的變化曲線與總衰減隨粒徑的變化曲線進(jìn)行對(duì)比，如圖12所示。

不同的模型尺寸H所對(duì)應(yīng)的散射衰減系數(shù)曲線如圖13所示，隨著H的增大，散射衰減系數(shù)略有減小，整體趨勢(shì)相同，相較于總的衰減系數(shù)大小來說，其不同H所產(chǎn)生的散射系數(shù)的差異所引起的總衰減系數(shù)的變化可忽略不計(jì)。

由圖12—13可知，由多次散射過程所造成的散射衰減相對(duì)于總衰減來說可以忽略不計(jì)，故其聲能量的損失主要由熱、黏性損失及其相互作用所引起的。

圖11 不同H對(duì)應(yīng)的聲衰減計(jì)算結(jié)果曲線Fig.11 Calculated attenuation curves corresponding to different H圖12 散射衰減曲線與總衰減曲線對(duì)比Fig.12 Scattering sttenuation curve compare with total ultrasonic attenuation

通過將Eloss取值為0，得到40 kHz驅(qū)動(dòng)頻率下，不同粒徑的散射衰減隨體積分?jǐn)?shù)的變化規(guī)律曲線如圖14所示，可以看出，對(duì)不同粒徑顆粒來說，散射衰減值隨體積分?jǐn)?shù)的增大先增大后減小。其模型計(jì)算的規(guī)律曲線與理論計(jì)算結(jié)果趨勢(shì)一致[15]。在體積分?jǐn)?shù)較小時(shí)，散射顆粒數(shù)量較少，即大部分聲能量可以直接透射被接收器接收，故散射衰減系數(shù)較小。而當(dāng)體積分?jǐn)?shù)不斷增大，散射顆粒數(shù)量逐漸增多，散射概率大大增加，但是多次散射作用不明顯，散射衰減系數(shù)增大。當(dāng)濃度增大到一定程度時(shí)，散射概率大大增加的同時(shí)，由于顆粒間的多次散射，且顆粒間的間距越來越小，使得散射部分的聲能量反而更少地被散射到發(fā)射與接收間的柱形區(qū)域外，導(dǎo)致散射衰減系數(shù)結(jié)果減小。

圖13 不同H值對(duì)應(yīng)的散射衰減結(jié)果曲線Fig.13 Scattering attenuation curve corresponding to different H圖14 不同粒徑下散射衰減系數(shù)隨體積分?jǐn)?shù)的變化曲線Fig.14 Change law of scattaring attenuation with volume concentration under different particle sizes

4 結(jié)論

1)通過程序的計(jì)算結(jié)果曲線與40、 60、 80 kHz的3種驅(qū)動(dòng)頻率的實(shí)驗(yàn)測(cè)量結(jié)果曲線間進(jìn)行擬合，單次散射能量損失分別為Eloss=5.5%，Eloss=6.04%和Eloss=7.37%，且相關(guān)系數(shù)均可達(dá)到0.99以上，證明此模型的過程設(shè)計(jì)符合超聲波在高含量氣固兩相中傳播特性，能夠用于高含量氣固兩相系統(tǒng)的衰減特性研究。

2)將Eloss取0得到的散射衰減系數(shù)計(jì)算結(jié)果表明，頻率為40 kHz，體積分?jǐn)?shù)約為63%的玻璃微珠-空氣高含量氣固兩相系統(tǒng)的聲衰減主要由于熱、黏性及其相互作用所產(chǎn)生，而由多次散射過程所產(chǎn)生的散射衰減幾乎可以忽略不計(jì)。

3)通過程序計(jì)算得到40 kHz驅(qū)動(dòng)頻率下，不同粒徑下散射衰減系數(shù)隨體積分?jǐn)?shù)的變化規(guī)律曲線，變化趨勢(shì)呈現(xiàn)先增大后減小。