張世瑋，李玉宇，孟磊，寧翔，蘇明旭

（1 上海理工大學(xué)上海市動力工程多相流動與傳熱重點實驗室，上海 200093；2 大唐環(huán)境產(chǎn)業(yè)集團股份有限公司，北京 100097）

隨著能源、化工、環(huán)保、生物、制藥和材料制備等行業(yè)的發(fā)展，兩相流中顆粒測量問題受到了國內(nèi)外學(xué)者的廣泛研究[1-2]。其中高濃度兩相流的在線測量越來越成為關(guān)注的焦點。例如，水煤漿作為一種新型潔凈燃料，其粒徑和濃度直接影響其穩(wěn)定性和燃燒效率[3]。河流泥漿中的粒徑和濃度變化可以控制河流含沙量輸運，以及監(jiān)控洪水或者泥石流帶來的危害[4]。燃煤電廠石灰石漿液粒徑直接影響煙氣脫硫效率和副產(chǎn)品石膏的質(zhì)量，漿液與煙氣的反應(yīng)效率隨固體顆粒粒徑的減小而提高，但粒徑過小會造成磨機過磨，增加設(shè)備的運行耗能[5-6]，對顆粒粒徑參數(shù)進(jìn)行在線測量具有重要意義。

篩分法和圖像分析是測量顆粒粒徑的常規(guī)方法，但均需對漿液進(jìn)行取樣，不適用于在線測量[7-8]。光散射法因其精度高在顆粒測量領(lǐng)域廣泛應(yīng)用，但由于光的穿透能力較弱，測量時需對漿液進(jìn)行稀釋，因此針對高濃度漿液粒徑的在線測量并不適用[9-10]。超聲波具有穿透能力強、非接觸、裝置簡單等優(yōu)點，尤其適用于較高濃度、光學(xué)不透明兩相流中顆粒粒徑的在線測量[11-13]。Falola等[14]針對超聲衰減譜方法ECAH模型中物性參數(shù)未知易導(dǎo)致多解的問題，提出全局最優(yōu)求解算法，經(jīng)過硅標(biāo)準(zhǔn)顆粒驗證后，用于SrCl2·6H2O 冷卻結(jié)晶過程的在線測量。Kitao等[15-16]提出了一種基于聚焦超聲的動態(tài)聲散射方法（dynamic sound scattering，DSS），基本思想是通過連續(xù)采集后向聲散射信號進(jìn)行相關(guān)分析，得到衰減分布函數(shù)，獲得顆粒體系平均擴散系數(shù)，進(jìn)而推算出納米量級顆粒粒徑。Wrobel等[17]利用脈沖寬頻超聲譜法測量高濃度液固兩相流，成功應(yīng)用于高濃度液固懸濁液和非均勻復(fù)合材料的在線測量。Abad等[18]利用多頻超聲衰減和散射信號測量高濃度漿液濃度、速度及顆粒粒徑分布，應(yīng)用于礦物泥漿、泥沙和城市污水研究，充分說明了超聲衰減譜測量高濃度漿料中顆粒粒徑和濃度的可行性。不過，對于流動狀態(tài)的高濃度漿液在線測量方法，尤其是動態(tài)測量和衰減譜誤差處理技術(shù)有必要開展進(jìn)一步研究。

本文建立描述高濃度漿液的非線性聲學(xué)測量模型，以石灰石漿液為例研究超聲法表征高濃度漿液粒徑分布的在線測量方法。設(shè)計并搭建超聲在線測量系統(tǒng)，基于管段透射法獲取漿液流動過程中的超聲透射信號和衰減譜。為減小衰減譜誤差對反演結(jié)果的影響，發(fā)展一種基于最優(yōu)正則化的加權(quán)最小二乘（optimum regularization technique weighted total least squares，ORT-WTLS）算法用于計算漿液兩相流中的顆粒粒徑分布，并使用圖像法對其測量結(jié)果進(jìn)行驗證。

1 顆粒粒徑聲學(xué)測量模型

1.1 耦合相模型及其散射疊加

奠定超聲法顆粒粒徑測量基礎(chǔ)的是ECAH（Epstein-Carhart-Allegra-Hawley）模型，該模型通過微體積元中質(zhì)量、動量和能量守恒定律出發(fā)，精確描述了球形顆粒兩相流中的聲波動行為，聲衰減和聲速的方程最終歸結(jié)為一個6階線性方程組的求解問題，最終得出的復(fù)波數(shù)方程如式(1)[19-20]。

式中，An為散射系數(shù)；φ為體積濃度；kc為連續(xù)介質(zhì)波數(shù)；R為顆粒半徑；k為兩相介質(zhì)的復(fù)波數(shù)，如式(2)。

式中，αs(ω)和cs(ω)分別為聲衰減系數(shù)和聲速；ω為角頻率；i為虛數(shù)單位。

由于ECAH模型并未考慮復(fù)散射和顆粒間相互作用對聲衰減的影響，因此該模型無法準(zhǔn)確預(yù)測高濃度顆粒兩相介質(zhì)中的聲波動行為。高濃度顆粒兩相流中需考慮顆粒間相互作用對超聲衰減的影響，其衰減特性隨濃度呈非線性相關(guān)[21]。Harker 和Temple（H & T）提出有別于ECAH 散射模型的觀點，根據(jù)水動力學(xué)角度出發(fā)提出耦合相模型的概念，獲得了描述連續(xù)相和顆粒相之間相互作用的黏性耗散方程，以及各相獨立的質(zhì)量和動量守恒方程。式(3)給出了復(fù)波數(shù)方程的表達(dá)式[22-23]。

式中，β和ρ分別為等溫壓縮系數(shù)和密度；顆粒相參數(shù)用上撇號表示；S是顆粒相濃度和粒徑、連續(xù)相黏度和密度的函數(shù)，如式(4)。

式中，δ為黏性集膚深度，如式(5)。

式中，η為連續(xù)相黏度。

式(3)建立了兩相流中顆粒相體積濃度、粒徑及物性參數(shù)與超聲衰減之間的定量關(guān)系。耦合相模型求解兩相空間的平均狀態(tài)方程，宏觀處理了液固兩相之間的界面，但并未對顆粒散射現(xiàn)象進(jìn)行描述，因此還需考慮顆粒散射對聲衰減的影響。Riebel 采用光學(xué)類比的方法給出BLBL（Bouguer-Lambert-Beer-Law）散射模型，其衰減系數(shù)的表達(dá)式描述為式(6)[24]。

式中，σ=ωR/c為顆粒量綱為1參量；c為連續(xù)相壓縮波聲速。對式(3)、式(6)給出的H & T和BLBL模型衰減系數(shù)求和得到描述高濃度漿液兩相流中的聲場波動行為。圖1為石灰石漿液體積分?jǐn)?shù)10%、顆粒直徑1～100μm、超聲頻率1～10MHz 時的模型矩陣。

圖1 石灰石漿液粒徑測量聲學(xué)模型矩陣

1.2 獨立模式反演算法

顆粒粒徑分布測量原理是基于不同粒徑產(chǎn)生的超聲譜特性差異，通過尋找理論譜和實驗譜之間的最佳擬合值確定。借助波的散射理論聯(lián)系單個顆粒的散射效應(yīng)和整個離散顆粒群的宏觀效果，得到懸濁液中的復(fù)波數(shù)，進(jìn)一步求解得到的超聲衰減系數(shù)表達(dá)如式(7)[25]。

式中，qi為顆粒半徑在第i個尺寸區(qū)間[Ri,Ri+1]的體積分?jǐn)?shù)；N為顆粒粒徑分散等級。

為獲得準(zhǔn)確的粒徑分布，將模型改寫為標(biāo)準(zhǔn)數(shù)學(xué)求解問題，改變系數(shù)矩陣和超聲譜向量，將式(7)離散并在多頻條件下轉(zhuǎn)化為線性方程如式(8)[26]。

式中，A為系數(shù)矩陣；F為待求的離散化顆粒尺寸數(shù)目分布；G為實驗超聲譜組成的向量。

該矩陣方程非適定且高度病態(tài)，常用矩陣求逆無法得到準(zhǔn)確結(jié)果，通常由最優(yōu)正則化（optimum regularization technique，ORT）算法結(jié)合非負(fù)最小二乘求解[27-28]，通過引入正則化因子γ和光順矩陣H，對病態(tài)方程進(jìn)行改善，提高解的穩(wěn)定性。表達(dá)如式(9)。

式中，AT為A的轉(zhuǎn)置矩陣；γ為正則化因子，表示矩陣的權(quán)重。

1.3 加權(quán)最小二乘算法求解粒徑分布

為減小顆粒兩相流中超聲衰減譜的測量誤差對粒徑分布求解精度的影響，發(fā)展一種基于最優(yōu)正則化的加權(quán)最小二乘算法求解顆粒系粒徑分布。引入權(quán)重矩陣W，粒徑分布求解公式表達(dá)如式(10)。

其中，超聲衰減譜權(quán)重矩陣W的計算表達(dá)如式(11)。

式中，ej為超聲衰減譜與其擬合譜的殘差絕對值；μ為ej的標(biāo)準(zhǔn)差。根據(jù)式(11)將超聲衰減譜的殘差分為三類：①殘差不是粗大誤差，殘差ej小于k0μ，權(quán)重因子為1；②殘差可能是粗大誤差，此時殘差ej在k0μ和k1μ之間，權(quán)重因子為k0/(|ej/μ|)；③殘差是粗大誤差，此時殘差ej大于k1μ，權(quán)重因子為0。通常取k0為1.5，k1為2.5[29]。

圖2為ORT-WTLS算法流程。由權(quán)重系數(shù)dw構(gòu)建權(quán)重矩陣W，并根據(jù)權(quán)重系數(shù)對實驗超聲衰減譜進(jìn)行修正，直至權(quán)重系數(shù)dw中無粗大誤差或殘差的標(biāo)準(zhǔn)差小于0.01結(jié)束迭代。ORT-WTLS算法令超聲衰減譜中具有粗大誤差的權(quán)重系數(shù)置零，能夠有效降低實驗衰減譜測量誤差對反演結(jié)果的影響。

圖2 ORT-WTLS算法流程

2 數(shù)值模擬

2.1 ORT-WTLS算法驗證

為驗證ORT-WTLS 算法的性能，以體積分?jǐn)?shù)1%的水-玻璃微珠懸濁液為研究對象，采用H & T和BLBL疊加模型分別計算顆粒直徑10～80μm懸濁液中的理論超聲衰減譜。采用ORT-WTLS 算法計算并統(tǒng)計不同算例粒徑分布體積中位徑DV50（小于該直徑的顆粒體積各占顆?？傮w積的50%），如圖3 所示。計算結(jié)果與設(shè)定值的相對誤差均小于5%，驗證了ORT-WTLS算法的準(zhǔn)確性。

圖3 ORT-WTLS計算結(jié)果與設(shè)定直徑對比

2.2 算法抗噪性分析

為模擬實驗超聲衰減譜引入的測量誤差，在理論譜中加入10dB 隨機噪聲，使用ORT-WTLS 算法計算顆粒系粒徑分布。圖4 為10dB 含噪衰減譜和由兩種算法得到的顆粒系粒徑分布。含噪衰減系數(shù)隨超聲頻率變化趨勢與理論衰減譜一致，但個別頻率下衰減系數(shù)的殘差為粗大誤差，致使圖4(b)中Twomey[30]算法得到的DV50與預(yù)設(shè)直徑差異較大，顆粒直徑10～40μm 時的相對誤差分別為55.3%、16.4%、23.8%、14.6%。而對超聲譜做加權(quán)處理的ORT-WTLS算法得到的特征參數(shù)DV50與設(shè)定直徑基本一致，不同粒徑下的相對誤差均小于5%。

圖4 不同粒徑含噪衰減譜及反演結(jié)果對比

分析ORT-WTLS 算法的魯棒性，統(tǒng)計連續(xù)20 組含噪衰減譜的反演結(jié)果。圖5 為10dB 含噪超聲衰減譜連續(xù)20 次計算粒徑分布DV50的統(tǒng)計結(jié)果。ORT-WTLS算法得到的特征參數(shù)DV50與設(shè)定直徑更接近，粒徑10～40μm時的相對誤差均小于5%，且標(biāo)準(zhǔn)差均小于1μm。顆粒直徑10μm 時，ORTWTLS算法得到的特征參數(shù)DV50與設(shè)定值基本一致，但Twomey 算法的計算結(jié)果偏差較大。當(dāng)顆粒直徑大于20μm 時，兩種算法得到的DV50與設(shè)定值基本一致，但ORT-WTLS 算法標(biāo)準(zhǔn)差明顯小于Twomey算法。因此ORT-WTLS 算法對超聲衰減譜的測量誤差有較好的魯棒性。

圖5 含噪衰減譜兩種算法計算結(jié)果對比

3 實驗測量及結(jié)果分析

3.1 在線測量系統(tǒng)及樣品配置

圖6為基于管段式透射法原理的漿液粒徑在線測量系統(tǒng)實物，主要包含超聲測量傳感器、超聲波發(fā)射接收儀（內(nèi)含數(shù)據(jù)采集卡，最高采樣率100MHz）、漿液循環(huán)裝置和上位機測量系統(tǒng)。傳感器測量區(qū)寬度18mm，在厚度7.5mm 的石英玻璃兩端分別安裝超聲發(fā)射、接收換能器（Olympus V310，中心頻率5MHz）。使用汕頭超聲CTS-04UT發(fā)射接收儀激勵換能器產(chǎn)生超聲波，另一換能器接收通過漿液后的超聲透射信號，波形數(shù)據(jù)的采樣率為100MHz。使用大唐環(huán)境產(chǎn)業(yè)集團股份有限公司某電廠提供的粉末配置不同質(zhì)量濃度的石灰石漿液。在7kg 純凈水中分別加入0.78kg、1.75kg、3.00kg、4.67kg 顆粒粉末配置質(zhì)量分?jǐn)?shù)（Cm）為10%～40%的石灰石漿液。分別將背景介質(zhì)水和被測漿液裝入漏斗中，打開球閥使水或漿液以重力作用流經(jīng)測量管段（流速3～4m/s）。通過網(wǎng)線將波形數(shù)據(jù)傳輸至上位機并計算實驗超聲衰減譜，結(jié)合ORT-WTLS 算法反演漿液體系的粒徑分布，實現(xiàn)漿液粒徑分布的動態(tài)測量。

圖6 漿液粒徑超聲在線測量系統(tǒng)

3.2 石灰石漿液粒徑在線測量

圖7(a)為背景介質(zhì)水在流動過程中不同時間的超聲透射信號，其一致性表明超聲測量系統(tǒng)具有較好的穩(wěn)定性。圖7(b)為質(zhì)量分?jǐn)?shù)10%漿液的超聲透射信號，由于固體顆粒的沉降現(xiàn)象，漿液流經(jīng)測量管段時的濃度隨時間的增大而減小，在t=5s時濃度為最低值，超聲透射信號幅值最大。漿液中顆粒相濃度越大，沉降速度越慢，濃度較小時信號幅值隨時間的增大現(xiàn)象明顯。圖8 為質(zhì)量分?jǐn)?shù)10%～40%漿液在t=3s時的超聲透射信號，信號幅值隨質(zhì)量濃度的增加明顯減小。液固兩相體系中顆粒相濃度增加到一定程度后，如漿液質(zhì)量分?jǐn)?shù)大于10%時，超聲透射信號幅值隨濃度呈現(xiàn)非線性變化。

圖7 不同測量時間的超聲透射信號

圖8 質(zhì)量分?jǐn)?shù)10%～40%漿液超聲透射信號

對背景介質(zhì)水和不同濃度漿液超聲透射信號作頻譜分析（fast Fourier transform，F(xiàn)FT），則超聲衰減系數(shù)的表達(dá)如式(12)[31]。

式中，Aw,i為水的透射信號在頻率fi的幅值；As,i為漿液透射信號在頻率fi的幅值；l2為測量區(qū)寬度。

圖9(a)為質(zhì)量分?jǐn)?shù)10%漿液在不同時間的超聲衰減譜，由于t=5s時超聲透射信號的幅值較大，衰減譜與其他時間相比在數(shù)值上存在較明顯差異，但t為1～5s時的歸一化衰減譜基本吻合。

圖9 不同測量時間和濃度的超聲衰減譜

圖9(b)為質(zhì)量分?jǐn)?shù)10%～40%漿液實驗衰減譜，同一濃度時衰減系數(shù)隨超聲頻率的增加而增大，漿液質(zhì)量分?jǐn)?shù)從10%增加至40%，頻率5MHz時的衰減系數(shù)由24Np/m 增加至139Np/m。質(zhì)量分?jǐn)?shù)30%和40%時，漿液對超聲高頻信號衰減嚴(yán)重，最大頻率由質(zhì)量分?jǐn)?shù)10%的9.4MHz 降至7.8MHz，但質(zhì)量分?jǐn)?shù)30% 和40% 時的有效頻率帶寬為6.7MHz，可用于漿液顆粒粒徑的反演計算。漿液質(zhì)量分?jǐn)?shù)40%時，超聲透射信號由背景介質(zhì)的7V降低至0.33V，但衰減系數(shù)隨超聲頻率的增加而增大，與模型預(yù)測的理論衰減譜一致。同時，以超聲頻率5MHz 為例，計算不同濃度漿液在時間t=1～5s時衰減系數(shù)平均值，當(dāng)質(zhì)量分?jǐn)?shù)為10%～40%時，漿液衰減系數(shù)并未隨濃度線性遞增，此時更適合采用H & T和BLBL非線性模型。

圖10為ORT-WTLS算法得到的質(zhì)量分?jǐn)?shù)10%～40%漿液在不同時間的粒徑分布。漿液質(zhì)量分?jǐn)?shù)10%時，粒徑分布雖稍微偏差，但其峰值對應(yīng)的顆粒直徑差異不大。分析質(zhì)量分?jǐn)?shù)20%時不同時刻的粒徑分布，5 個測量時間下共得到3 個不同的粒徑分布，其峰值對應(yīng)的顆粒直徑與質(zhì)量分?jǐn)?shù)10%時基本一致。漿液質(zhì)量分?jǐn)?shù)30%時，除t=5s 時粒徑分布峰值對應(yīng)的顆粒直徑稍偏小外，其余時間下的粒徑分布基本吻合。質(zhì)量分?jǐn)?shù)40%時，不同時間的粒徑分布雖稍有偏差，但粒徑分布峰值對應(yīng)顆粒直徑均在16μm左右。

圖10 不同濃度漿液在不同時刻的粒徑分布

3.3 算法實驗驗證和測量誤差分析

將上述實驗所用石灰石漿液進(jìn)行稀釋，使用丹東百特儀器有限公司BT-1600光學(xué)顯微鏡觀察漿液顆粒的大小及形態(tài)，圖11 為光學(xué)顯微鏡圖片。使用ImageJ 圖像處理軟件綜合分析多張圖片，石灰石漿液粒徑分布特征參數(shù)見表1。

表1 圖像法分析漿液粒徑分布特征參數(shù)

圖11 石灰石漿液顯微鏡圖像

使用第二節(jié)中的兩種反演算法對質(zhì)量分?jǐn)?shù)10%和30%漿液在時間t=1～5s 時的衰減譜進(jìn)行計算，ORT-WTLS 算法得到的體積中位經(jīng)與圖像法的相對誤差分別為3.81%和6.35%，而Twomey 算法的相對誤差為-4.87%和10.12%，表明ORTWTLS算法能夠有效提高漿液顆粒粒徑的在線測量精度。采用圖像法對超聲在線測量結(jié)果的準(zhǔn)確性進(jìn)行驗證。由于石灰石漿液測量過程中的固體顆粒沉降現(xiàn)象，不同時間的粒徑分布及特征參數(shù)稍有差異，計算時間t=1～5s 的平均值與圖像法測量結(jié)果進(jìn)行對比。表2 為粒徑分布特征參數(shù)DV50、DV90和分布寬度DV90-DV10的誤差統(tǒng)計結(jié)果。

表2 粒徑分布特征參數(shù)誤差統(tǒng)計結(jié)果

漿液質(zhì)量分?jǐn)?shù)10%時的測量誤差小于5%，質(zhì)量分?jǐn)?shù)20%～40%時相對誤差稍有增加，但均小于7%。漿液粒徑分布寬度的測量誤差相對較大，質(zhì)量分?jǐn)?shù)10%時誤差達(dá)到15.25%。后續(xù)研究考慮將多種群智能優(yōu)化算法應(yīng)用于多相流動過程中顆粒粒徑分布雙參數(shù)反演，以期得到更為準(zhǔn)確的顆粒粒徑及分布寬度。

4 結(jié)論

針對高濃度漿液兩相流超聲衰減譜的誤差處理及粒徑分布在線測量展開研究，建立顆粒測量的非線性聲學(xué)理論模型并研制漿液顆粒粒徑超聲在線測量系統(tǒng)，發(fā)展了一種基于加權(quán)最小二乘的ORTWTLS算法求解顆粒兩相流中的粒徑分布，主要結(jié)論如下。

（1）以不同粒徑理論衰減譜和10dB 含噪衰減譜為例分析了ORT-WTLS算法的準(zhǔn)確性和抗噪性，數(shù)值模擬算例表明ORT-WTLS 算法計算結(jié)果與設(shè)定值相對誤差小于5%。

（2）研制漿液細(xì)度超聲在線測量系統(tǒng)，采集質(zhì)量分?jǐn)?shù)10%～40%漿液在流動過程中的超聲透射信號并計算衰減譜，頻率5MHz處時衰減系數(shù)并未隨濃度線性遞增，表明高濃度漿液粒徑測量采用H & T和BLBL模型是合適的。

（3）漿液質(zhì)量分?jǐn)?shù)40%時，仍可得到有效帶寬6.7MHz 的超聲衰減譜，有效衰減系數(shù)的最大值可達(dá)240Np/m，表明超聲法測量高濃度漿液兩相流顆粒粒徑分布的可行性和準(zhǔn)確性。

（4） 通過ORT-WTLS 算法得到的質(zhì)量分?jǐn)?shù)10%～40%漿液粒徑分布與圖像法測量結(jié)果吻合，粒徑測量誤差小于7%。結(jié)合管段透射法超聲在線測量系統(tǒng)與ORT-WTLS 算法可實現(xiàn)高濃度漿液兩相流粒徑在線測量。

符號說明

A—— 模型系數(shù)矩陣

An—— 模型散射系數(shù)

Aw,i—— 背景介質(zhì)透射信號在頻率fi時的信號幅值，V

As,i—— 被測漿液透射信號在頻率fi時的信號幅值，V

c—— 連續(xù)相壓縮波聲速，m/s

cs—— 超聲波在兩相/多相介質(zhì)中的傳播相速度，m/s

fi—— 頻率，MHz

F—— 待求的離散化顆粒尺寸數(shù)目分布

G—— 實驗超聲衰減系數(shù)組成的向量，Np/m

H—— 光順矩陣

I—— 單位矩陣

k—— 兩相介質(zhì)的復(fù)波數(shù)

kc—— 連續(xù)介質(zhì)波數(shù)

l2—— 測量區(qū)寬度，m

N—— 顆粒粒徑分散等級

qi—— 顆粒半徑在第i個尺寸區(qū)間[Ri,Ri+1]的體積分?jǐn)?shù)

R—— 顆粒半徑，m

W—— 超聲衰減譜誤差權(quán)重矩陣

αs—— 超聲衰減系數(shù)，Np/m

φ—— 體積濃度

ρ—— 密度，kg/m3

σ—— 顆粒量綱為1參數(shù)

γ—— 正則化因子