多聲道超聲波流量計內部流-聲耦合特性仿真分析

馬雪林，徐 雅，謝代梁，徐志鵬，曹松曉，胡鶴鳴

（1.中國計量大學 浙江省流量計量技術重點實驗室，浙江 杭州 310018; 2.中國計量科學研究院，北京 100029）

0 引 言

多聲道超聲波流量計是通過測量超聲波在流體中的順、逆流時間差，計算聲道上的平均流速，再對不同聲道高度的平均流速進行積分，從而求得流量測量的非接觸式流量計[1]。與傳統(tǒng)流量計比較，它具有測量精度高、無阻流部件、量程比寬，更能適應被測流體溫度、壓力、密度等參數(shù)的變化等優(yōu)點，因此近年來在各個領域得到了廣泛的應用[2]。目前，多聲道超聲波流量計主要應用于國內外大型水電站大口徑輸水管道的流量計量，在在線檢測水輪機效率和狀態(tài)方面起到了重要作用[3]。

為了有效的降低超聲波流量計的測量誤差，提高測量精度，很多學者對其進行了研究。天津大學的SUN Y等[4]將計算流體力學、超聲波學和射線聲學理論相結合，提出了一種新的超聲波流量計流-聲分析的混合方案。范思航等[5]在多普勒效應的理論基礎上，分析了聲波在管道中的傳播過程等，進而計算得到相應的技術參數(shù)，設計出一款適用于多普勒流量計的超聲換能器。IOOSS B[6]在射線追蹤法的理論基礎上，推導了在二維模型下聲波在充分發(fā)展的流場中的傳播微分方程，利用Runge-Kutta進一步求解，得出了雷諾數(shù)在2.5×104與7.5×107之間的變化，以及聲波傳播的非線性軌跡所導致的測量誤差。趙文等[7]為探究超聲流量計換能器大小及安裝結構所引起的擾流效應，對換能器造成的擾流影響進行了數(shù)值仿真，得到了小口徑超聲波流量計換能器擾流造成的系統(tǒng)誤差。BEZDěK M等[8]將有限元法和Helmholtz積分射線追蹤法結合，使用FEM計算聲波的發(fā)射和接收，解決RTM不能很好計算管壁交界處波形轉換的問題。

實際上，由于是流場和聲場耦合作用共同形成了超聲波流量計的流量測量誤差。但是根據(jù)以上研究可知，近年來針對多聲道超聲波流量計的研究大多集中在其內部流場及其影響，針對流-聲耦合特性的研究較少。因此本文針對多聲道超聲波流量計流量測量的流場和聲場耦合問題，采用多物理場數(shù)值仿真軟件COMSOL對這一問題展開研究，探究了多聲道測量管道內部的流場分布以及不同聲道下聲波在管道流動中的耦合動態(tài)傳播過程，并通過互相關法計算誤差，探討了聲波傳播造成的擾動對流量測量的影響。

1 計算模型

本文流場和聲場耦合方法的數(shù)值模擬計算通過有限元軟件COMSOL Multiphysics完成，主要分為流場計算和聲場計算兩個部分。在研究流聲耦合動力學行為時，我們一般視流體為穩(wěn)定流動，聲波傳播對流體參數(shù)形成一定擾動影響造成系統(tǒng)誤差[9]。為了獲得聲波在穩(wěn)態(tài)背景流場下的瞬態(tài)傳播過程，超聲波流量計模型中的穩(wěn)態(tài)背景流場采用不可壓縮流動k-ε湍流模型模擬獲得，而瞬態(tài)高頻聲場則通過超聲波接口下時域求解器聯(lián)合流場求解模塊模擬實現(xiàn)。

本文所采用的仿真模型是基于一臺DN200的5聲道超聲流量計建立，即為交叉布置的2聲道和3聲道，采用圓柱形換能器，如圖1所示。該流量計實際口徑為200 mm，聲道角度為45°，換能器直徑是30 mm。

圖1 仿真模型（單位：m）

2 數(shù)值模擬計算

為了避免管道入口流體的不穩(wěn)定對測量精度的影響，則需要保證流體經(jīng)過超聲波流量計測量管道時已經(jīng)充分發(fā)展，所以在進行穩(wěn)態(tài)流場的模擬計算時，將設置流體已經(jīng)充分發(fā)展為入口邊界條件。

2.1 聲場的瞬態(tài)計算

超聲信號在流體的傳播中必然滿足基本的物理定律，即連續(xù)性方程、歐拉方程、物態(tài)方程：

其中ρ、p和u分別為流場的密度、壓強和速度參數(shù)。

超聲波流量計的工作原理主要是將超聲波作為信息的載體，測量其在介質中的順、逆流時間差，是小振幅的聲波，所以，在描述聲波在流體中的傳播過程時，對應的壓強、速度、密度等物性參數(shù)可以做以下線性化處理：

其中帶有下標“0”的表示穩(wěn)態(tài)背景流場的參數(shù)，帶有小下標“1”表示聲場的擾動量。

將上述公式結合處理后保留一階小擾動量，超聲傳播方程則表示為：

其中ρ0、p0和u0分別表示穩(wěn)態(tài)背景流場下的密度、壓強和速度參數(shù)，相應的聲波傳播對以上參數(shù)的影響分別表示為ρ1、p1和u1的小擾動量。

對基于穩(wěn)態(tài)背景流場下的超聲傳播方程進行間斷有限元法的空間離散和Runge-Kutta進一步的時間推進計算，結合相應的邊界條件即可獲得聲波在穩(wěn)態(tài)背景流場下的瞬態(tài)傳播過程。

2.2 聲學邊界條件

在聲傳播過程中，在信號管的換能器（信號）端指定超聲詢問信號。傳播時間法要求從兩端發(fā)送信號，因此在實際設備中，這兩端既是發(fā)射器也是接收器。在模擬計算時，只需將源邊界和接收器邊界互換，從而實現(xiàn)一對換能器交替發(fā)射和接收對方的超聲波信號[10]。詢問信號是采用高斯脈沖調制的正弦波，通過在源邊界指定法向速度vn（t）來模擬實現(xiàn)：

式中：f0——振動頻率；

A——信號幅值；

ω0——振動角頻率；

T0——振動周期。

取值如下：A=0.1 mm，ω0=2πf0，f0=1 MHz，T0=1/f0。

為了防止聲波在流量計管段中無限傳播，在出入口經(jīng)反射回到測量段，對測量精度產(chǎn)生影響，則在超聲波流量計模型的前后段分別設置100 mm的吸收層，使聲波能夠在傳播到出入口時能夠被有效的衰減，過濾及吸收。在吸收層的外邊界加上聲學阻抗：

式中：n——曲面法向量；

Z——介質水的特征阻抗。

2.3 網(wǎng)格劃分和時間步長的選取

在流場和聲場計算時采用兩種不同的網(wǎng)格劃分方式，如圖2所示分別為流場和聲場的網(wǎng)格。流場網(wǎng)格在換能器安裝方式所形成的凹凸部位進行加密處理，在采用流場求解模塊獲得超聲波傳播的穩(wěn)態(tài)背景流場后，由于流場網(wǎng)格和聲場網(wǎng)格的不一致性，通過內置的拉伸算子將流場信息插值到聲場網(wǎng)格上。聲場計算過程中，使用間斷伽遼金法（GDM）進行接口建模時，為了避免產(chǎn)生小尺寸的網(wǎng)格，取得較好的仿真結果，聲學網(wǎng)格采用自由四面體網(wǎng)格進行劃分，將網(wǎng)格尺寸控制在λ/1.5～λ/2，同時采用單元質量優(yōu)化功能，時間步長為Δt=T0/12，為了確保聲波在管道內可以完成從換能器發(fā)射端到達接收端的傳播過程，計算周期設置為220。

圖2 流場和聲場網(wǎng)格

本文的所有模擬計算均在實驗室的工作站上進行，其配置如下：2個 Intel Xeon E5-4 667 v4 @2.20 GHz CPU和128 GB（16×8 GB）內存。

3 計算結果及分析

3.1 穩(wěn)態(tài)背景流場

圖3所示為計算得到的v=3 m/s不同聲道截面流速分布云圖。從圖中看出換能器安裝效應對超聲波流量計的穩(wěn)態(tài)背景流場有一定影響。流體在換能器處產(chǎn)生滯流，形成低速區(qū)，這一點還可以從不同聲道的軸向速度分布圖得出（如圖4所示）。而且因為A1、A2、B2聲道的換能器分別有不同程度的嵌入，在該聲道換能器凸起或凹陷附近又形成一個低速的旋渦，上游換能器形成的低速區(qū)正好在聲波的傳播路徑內，聲波從發(fā)射到管道內部要經(jīng)過一個較大的流速變化，對聲波傳播產(chǎn)生了一定影響，對下游換能器來說，旋渦則出現(xiàn)在其后側，對聲波傳播的影響較小。

圖3 不同聲道的流速分布

圖4 測量段各聲道的軸向速度分布

圖4所示為背景流速為3 m/s時測量段各聲道的軸向速度分布情況，為了便于觀察分析，將對不同聲道下的測量結果進行如下的歸一化處理：流體流速/入口流速。

從圖中看出在換能器附近的區(qū)域產(chǎn)生了流速分布不均勻情況，這表明在上下游換能器的安裝結構內部產(chǎn)生了漩渦，且由于換能器安裝方式及聲道高度不同，各聲道軸向速度的分布及最大值受到影響較大，聲道高度越大則軸向速度最大值越小，分布范圍越窄。

通常，這是因為換能器上游流經(jīng)管道壁面的流體在突然經(jīng)過換能器安裝所引起的凸起或凹陷部分時，受到剪切力的影響而在換能器附近形成一定的速度梯度。當流體流速較大或者換能器嵌入較深的情況下，這個速度梯度將會更大，流體自身相互作用而產(chǎn)生流體動力聲源，作用于換能器安裝結構內部的滯留流體發(fā)生振動，形成噪聲，從而對超聲波流量計的測量精度產(chǎn)生一定的影響。

3.2 瞬態(tài)聲場

圖5展示了背景流速為2 m/s時，在不同時刻下B2聲道截面聲波在管道內部的傳播過程。

圖5 B2聲道截面聲波的瞬態(tài)傳播過程

以B2聲道截面為例，從圖中看出聲波在流體中以球面波的形式傳播，經(jīng)過一段時間，從測量管段的限制中脫離進入傳感器主軸管段內部。同時相鄰聲道換能器發(fā)射的聲波有少量以同樣的形式在B2截面內傳播。隨著時間的推進，聲波不斷前移，則會與相鄰聲道發(fā)射的少量聲波發(fā)生相交，以相互作用的形式繼續(xù)前移，如圖中所示，約在0.136 ms分別到達基表的下壁面和上壁面，兩道聲波與管壁相互作用產(chǎn)生第一次反射，與原始發(fā)射聲波在主軸管段內產(chǎn)生復雜的疊加并繼續(xù)傳播。約在0.165 ms時刻，聲波傳遞至接收換能器并與經(jīng)壁面反射的聲波進行復雜疊加。傳播至傳感器出入口的聲波及反射波均被兩側的吸收層吸收，對測量管段內的聲波沒有產(chǎn)生影響。

通過B2截面不同時刻聲波在管道內部的傳播過程可以看出，一次聲波傳輸過程中，聲波在管道內部壁面、換能器端面處發(fā)生多次復雜的反射及疊加，并與相鄰聲道發(fā)射的少量聲波相互作用影響，這將對接收換能器檢測到聲波信號產(chǎn)生一定的影響，從而形成系統(tǒng)誤差。

因此，在實際的超聲波流量計設計中，信號采集電路中添加了窗口延遲功能來避免上述提到的反射聲波先于原始發(fā)射聲波到達接收換能器從而影響到傳播時間差的確認[7]。

圖6展示了背景流速為2 m/s時，0.136 ms時不同聲道截面超聲波在管道內部的傳播過程。從圖中可以看出當聲波在某一時刻突破測量管段的限制以球面波的形式不斷前移進入主軸管段，將會與其他聲道發(fā)射的少量超聲波相遇、疊加繼續(xù)隨著時間推移前進。越靠近管道軸心的聲道上的超聲波因為傳播路徑長的緣故而到接收換能器的用時越長。

圖6 同一時刻下不同聲道聲波傳輸過程

圖7展示了背景流速為3 m/s時，聲波傳播至管道中心處和聲波到達接收端換能器時B2聲道接收端換能器端面的聲壓分布情況。從圖中可以看出，不同階段下接收換能器端面的聲壓分布并不均勻且相對復雜，端面周邊的聲壓與中心位置的聲壓相比較大。這主要是因為聲波在傳播階段進行多次復雜的反射及疊加而形成的。

圖7 聲道B2中段和末段換能器聲壓分布

圖8展示了背景流速為3 m/s順流情況下，0.165 ms時，不同聲道接收端換能器端面的聲壓分布對比情況。

從圖8中可以看出在順流傳播末段時，近似以管道軸心對稱分布的聲道，其接收換能器的聲壓分布也是基本對稱的，換能器靠近軸心的邊緣聲壓分布更復雜，且離軸心越遠的聲道聲波更早到達接收換能器。

圖8 傳播末段聲壓分布

3.3 數(shù)據(jù)分析

因為換能器端面的聲壓分布并不均勻，將通過平均聲壓采用互相關法計算傳播時間差。

如圖9所示是B2聲道上的一對接收換能器記錄的信號順、逆流傳播時的平均聲壓，采集過程通過切換換能器的源和接收器來完成。從圖中可以看出，平均聲壓變化的曲線，在逆流情況下，因為流體運動對聲波傳播產(chǎn)生一定的阻礙作用，聲波的反射及復雜疊加的影響，使換能器接收到聲波的波形比順流時較晚。

圖9 聲壓信號

由互相關原理可知，當互相關函數(shù)取得最大值時，其時間位移就是順逆流情況下，接收換能器接收聲波的波形時間差。假設在m0點處取得函數(shù)最大互相關，可以求得聲波傳播時間差Δt[11-12]：

在本文中測量介質為水，則超聲波在流體中的傳播速度c0=1 481 m/s，根據(jù)超聲波流量計的測量原理，聲速遠大于流體流速，故可以將測量原理數(shù)學表達式簡化為：

其中Li為聲道i的理論長度。將Δt代入上式中，即可解得聲道i聲道方向的面平均速度。

時差法超聲波流量計聲道上的面平均流速與管道內的體平均流速之間存在如下關系：

其中K為流速修正系數(shù)。

換算體積流量為：

式中：D——管道直徑；

Vi——聲道i上的線平均流速；

wi——聲道i的權重系數(shù)。

將K值設為1，可求得超聲流量計測得的體積流量Q。系統(tǒng)計算誤差可以用以下公式表示[13]：

式中：Qt——理論流量；

Vin——仿真中設置的入口速度。

由表1中的仿真結果可得出如下規(guī)律：

表1 不同流速下的測量誤差

1）斜插式管段多聲道超聲波流量計，入口速度為1～3 m/s時，系統(tǒng)計算誤差ε為負；

2）流量計管道直徑不變的情況下，系統(tǒng)計算誤差絕對值隨著入口速度增大（即Re增大）而增大。

4 結束語

通過有限元軟件COMSOL Multiphysics對一臺DN200的5聲道超聲波流量計進行了流-聲耦合的多物理場數(shù)值模擬研究，獲得了不同聲道不同換能器安裝條件下的流場分布特性、其耦合特性下聲波傳播的動力學過程及換能器接收端面的聲壓分布情況，分析了不同流速下測量誤差的變化趨勢。主要得出結論：

1）由于換能器安裝效應造成管道發(fā)生凸起或凹陷影響背景流場分布不均勻，在換能器端面附近形成低速區(qū)的旋渦，導致各聲道的線平均流速減小，積分體積流量比理論流量值小，則系統(tǒng)計算誤差為負。

2）超聲波在流量計管道內部傳播過程中，發(fā)生多次復雜的反射與疊加，導致接收換能器端面接收實際聲壓信號時產(chǎn)生誤判，從而使聲波傳播時間差有一定的誤差。

3）在該模型的換能器安裝條件下，入口速度為1～3 m/s時，系統(tǒng)計算誤差為-1.85%～-0.98%，系統(tǒng)計算誤差隨著入口速度增大（即Re增大）而增大。