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

祖曉萌， 朱曙光

(南京理工大學能源與動力工程學院， 江蘇南京210018)

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

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

1 蒙特卡羅模型建立

1.1 模型的實現過程

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

1.2 透射概率

兩相系統的模型是邊長為H=0.02 m的正方體，通過橫截面的縱向顆粒數即可計算顆粒體積分數，

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

(1)

式中：φi為不同粒徑顆粒體積分數；Ri為顆粒半徑；rnum為橫截面的縱向顆粒數；Vs為顆粒所占總體積；Va為兩相系統的總體積。

將后續實驗中不同粒徑顆粒進行單層堆積時的體積分數進行測量，體積分數的實測值見表1，具體計算如式(2)。

(2)

式中：Vr為顆粒的實際體積；Vb為顆粒表觀體積；m為總的顆粒質量；ρr為顆粒密度；ρb為顆粒表觀密度。根據ρr和ρb便可得到顆粒堆積層顆粒體積分數。

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

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

表1 不同顆粒粒徑的體積分數測量值

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

(3)

(4)

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

εn為0到1的隨機數，由εn和P的關系來判斷聲子是否發生透射，即

εn

(5)

εn≥P散射。

(6)

1.3 聲子能量損失

對于高含量氣固兩相系統來說，多次散射過程除了會造成散射衰減，還會因多次散射而產生更多的熱衰減和黏性衰減。當體積一定時，顆粒體積分數和顆粒數量成正比。對低含量氣固兩相系統來說，聲波傳播過程中，每個顆粒只發生單次散射，并產生熱黏性衰減，不考慮顆粒間相互作用，所以總的聲衰減與體積分數成正比。而當體積分數不斷增大，會產生熱、黏性相互作用和多次散射，將導致系統中單個顆粒多次與聲波相互作用產生更多的熱黏性損失，進而使得總的聲衰減與體積分數呈現非線性的對應關系。同時對于高含量系統來說，聲衰減主要取決于黏性衰減和熱衰減，散射衰減則要小得多[9]。

1.3.1 熱衰減和黏性衰減

EACH模型[10-11]較為全面地考慮了熱傳導、黏性和散射對聲傳播特性的影響，當聲波入射到兩相界面，發生散射的同時會產生熱衰減和黏性衰減，并認為熱衰減和黏性衰減主要是由于兩相介質的可壓縮性和密度差異所導致的。

基于上述理論，模型將由于顆粒自身及其相鄰顆粒間的熱、黏性相互作用產生的聲衰減過程與散射過程相結合，通過單次散射能量損失Eloss對其衰減強度進行表征。高含量兩相系統中，就單個顆粒來看，小顆粒與相鄰顆粒接觸面積更小，作用強度低，但是相鄰顆粒數量更多，所以本文中對不同粒徑顆粒發生散射進行建模時，其單次散射能量損失Eloss取值相同且為常數。Eloss中包含了顆粒間相互作用引起的熱衰減和黏性衰減，因此在模型計算時，僅對單個顆粒進行判斷即可，無需考慮其他顆粒對其的影響。

1.3.2 散射衰減

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

1.4 散射角度

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

(7)

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

(8)

散射角度θ取m+1進行散射[14]。

1.5 聲子位置坐標

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

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

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

當0≤θi≤45 °或315 °<θi≤360 °，

(9)

當45 °<θi≤135 °，

(10)

當135 °<θi≤255 °，

(11)

當225 °<θi≤315 °，

(12)

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

當0≤θi≤45 °或315 °≤θi≤360 °，

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

(13)

當45 °<θi≤135 °，

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

(14)

當135 °<θi≤225 °，

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

(15)

當225 °<θi≤315 °，

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

(16)

1.6 聲衰減系數計算方法

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

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

(17)

(18)

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

1.7 蒙特卡羅程序算法

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

2 實驗

2.1 裝置

圖5為實驗裝置設計圖，工控機發出指令控制驅動器驅動發射換能器發出超聲波，同時發送指令給采集器準備進行信號采集，聲波穿過顆粒堆積層后被接收換能器接收，通過濾波放大后，經采集器送給工控機進行進一步數據處理。

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

2.2 結果

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

Pi=Poe-αx，

(19)

Vi=Voe-αx，

(20)

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

利用接收電壓的對數值與堆積厚度進行線性擬合的斜率即為聲衰減系數。40 kHz的驅動頻率下不同粒徑玻璃微珠所對應的線性擬合曲線見圖6，衰減系數值見表3。

圖5 實驗系統設計圖Fig.5 Design picture of experimental equipment

3 討論

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

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

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

圖7 不同Eloss值對應聲衰減規律曲線與實驗結果對比Fig.7 Calculated attenuation curves correspouding to different Eloss compare to experimental results圖8 實驗數據與模型計算結果的曲線最佳擬合Fig.8 The best fitted curve of model calculated compare to experimental curve

圖9 實驗數據與模型計算結果的曲線最佳擬合Fig.9 The best fitted curve of model calculated compare to experimental curve圖10 實驗數據與模型計算結果的曲線最佳擬合Fig.10 The best fitted curve of model calculated compare to experimental curve

聲衰減系數的計算公式中存在模型尺寸變量H，計算不同H所對應的聲衰減系數，其結果曲線如下圖11所示，可以看出，聲衰減系數的計算結果不會隨模型尺寸H的改變而變化，即可認為聲衰減僅與兩相系統內部介質參數有關。

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

當Eloss取值為0時，通過程序計算，可得到40 kHz驅動頻率下的散射衰減系數，將散射衰減隨粒徑的變化曲線與總衰減隨粒徑的變化曲線進行對比，如圖12所示。

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

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

圖11 不同H對應的聲衰減計算結果曲線Fig.11 Calculated attenuation curves corresponding to different H圖12 散射衰減曲線與總衰減曲線對比Fig.12 Scattering sttenuation curve compare with total ultrasonic attenuation

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

圖13 不同H值對應的散射衰減結果曲線Fig.13 Scattering attenuation curve corresponding to different H圖14 不同粒徑下散射衰減系數隨體積分數的變化曲線Fig.14 Change law of scattaring attenuation with volume concentration under different particle sizes

4 結論

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

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

3)通過程序計算得到40 kHz驅動頻率下，不同粒徑下散射衰減系數隨體積分數的變化規律曲線，變化趨勢呈現先增大后減小。