基于蒙特卡羅法模擬高含量氣固兩相流聲衰減

祖曉萌， 朱曙光

(南京理工大學(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ù)雜抽象的物理過(guò)程簡(jiǎn)單具體化而得到快速發(fā)展， 并被成功地應(yīng)用于預(yù)測(cè)高含量液固兩相系統(tǒng)的聲衰減。本文中利用蒙特卡羅方法對(duì)超聲波透射、散射和衰減的發(fā)生概率和作用方式進(jìn)行建模，來(lái)模擬高含量氣固兩相系統(tǒng)聲衰減。

1 蒙特卡羅模型建立

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

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

1.2 透射概率

兩相系統(tǒng)的模型是邊長(zhǎng)為H=0.02 m的正方體，通過(guò)橫截面的縱向顆粒數(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è)值見(jiàn)表1，具體計(jì)算如式(2)。

(2)

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

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

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

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

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

(3)

(4)

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

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

εn

(5)

εn≥P散射。

(6)

1.3 聲子能量損失

對(duì)于高含量氣固兩相系統(tǒng)來(lái)說(shuō)，多次散射過(guò)程除了會(huì)造成散射衰減，還會(huì)因多次散射而產(chǎn)生更多的熱衰減和黏性衰減。當(dāng)體積一定時(shí)，顆粒體積分?jǐn)?shù)和顆粒數(shù)量成正比。對(duì)低含量氣固兩相系統(tǒng)來(lái)說(shuō)，聲波傳播過(guò)程中，每個(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)來(lái)說(shuō)，聲衰減主要取決于黏性衰減和熱衰減，散射衰減則要小得多[9]。

1.3.1 熱衰減和黏性衰減

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

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

1.3.2 散射衰減

與熱衰減和黏性衰減的能量損失機(jī)制相比，散射衰減的機(jī)理有所不同。顆粒對(duì)聲波進(jìn)行散射時(shí)不會(huì)產(chǎn)生能量損失，它類似于光散射，只是使得部分聲能量沒(méi)有被換能器接收到而產(chǎn)生聲衰減[12]。當(dāng)Eloss取值為0時(shí)，即在聲子的整個(gè)傳播過(guò)程均沒(méi)有能量損失，只改變傳播路徑，通過(guò)程序計(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)行透射或單次散射后，需要確定過(guò)程結(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)聲子單次散射過(guò)程結(jié)束后位置確定，見(jiàn)式(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)聲子單次透射過(guò)程結(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)過(guò)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ì)聲子在傳播過(guò)程中的運(yùn)行狀態(tài)進(jìn)行判定，引入了Nin吸收聲子數(shù)和Nout逃逸聲子數(shù)(即沒(méi)有被接收到的聲子數(shù))，利用MATLAB進(jìn)行程序編寫(xiě)，程序框圖如圖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)采集，聲波穿過(guò)顆粒堆積層后被接收換能器接收，通過(guò)濾波放大后，經(jīng)采集器送給工控機(jī)進(jìn)行進(jìn)一步數(shù)據(jù)處理。

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

2.2 結(jié)果

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

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

3 討論

通過(guò)改變程序中單次散射能量損失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%之間取值來(lái)進(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)致的聲衰減變化值相較于熱衰減、黏性衰減和散射衰減來(lái)說(shuō)近乎可以忽略，則因頻率增大所引起的聲衰減值的增大主要是由黏性衰減、熱衰減和散射衰減引起的，單次散射能量損失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í)，通過(guò)程序計(jì)算，可得到40 kHz驅(qū)動(dòng)頻率下的散射衰減系數(shù)，將散射衰減隨粒徑的變化曲線與總衰減隨粒徑的變化曲線進(jìn)行對(duì)比，如圖12所示。

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

由圖12—13可知，由多次散射過(guò)程所造成的散射衰減相對(duì)于總衰減來(lái)說(shuō)可以忽略不計(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

通過(guò)將Eloss取值為0，得到40 kHz驅(qū)動(dòng)頻率下，不同粒徑的散射衰減隨體積分?jǐn)?shù)的變化規(guī)律曲線如圖14所示，可以看出，對(duì)不同粒徑顆粒來(lái)說(shuō)，散射衰減值隨體積分?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í)，由于顆粒間的多次散射，且顆粒間的間距越來(lái)越小，使得散射部分的聲能量反而更少地被散射到發(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)通過(guò)程序的計(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以上，證明此模型的過(guò)程設(shè)計(jì)符合超聲波在高含量氣固兩相中傳播特性，能夠用于高含量氣固兩相系統(tǒng)的衰減特性研究。

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

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