水下穩(wěn)態(tài)渦流場聲傳播的數(shù)值模擬和實(shí)驗研究

喻敏，樊丁繁，張燁，劉航

(1.武漢理工大學(xué)船海與能源動力工程學(xué)院，湖北武漢 430074；2.中國船舶集團(tuán)有限公司系統(tǒng)工程研究院，北京 100094)

0 引 言

聲波探測作為一種直接的、非侵入式的和全局式的測量方法，可用于水下旋渦特性的測量[1-2]。聲波在穿過渦流場時，其傳播路徑發(fā)生變化[3-4]，由此產(chǎn)生聲波相位的變化，相位變化量和渦流場的形態(tài)特征密切相關(guān)。因此，發(fā)展基于渦流場聲傳播特性發(fā)展水下渦流場探測識別技術(shù)成為可能，在軍事和工程中具有應(yīng)用前景。

使用聲波探測方法進(jìn)行渦特性的計算首先在超聲領(lǐng)域得到了應(yīng)用。1992年，文獻(xiàn)[5-6]提出了在非均勻介質(zhì)(渦流場)中，引入光學(xué)中相位共軛鏡的方法對目標(biāo)的反射聲進(jìn)行聚焦，以達(dá)到定位的目的。1997年，Roux等[7]通過實(shí)驗從接收信號的幅值與相位差兩個方面研究聲信號穿過渦流場后產(chǎn)生的變化。從1998年開始，文獻(xiàn)[8-10]先后使用多種實(shí)驗裝置對水下渦流場進(jìn)行深入研究，指出渦流同時改變了聲波的傳播方向和傳播速度。2001 年，Pag‐neux等[11]基于射線聲學(xué)研究了二維空間中的理想渦對高頻聲波的散射。2005年，Rosny等[12]針對低馬赫數(shù)渦流場引起的聲信號相位跳變量較小、不易測量的問題，采用了TRM 技術(shù)來增大相位跳變量、降低干擾，取得了較好的效果。2023 年，陳雨晨等[13]針對介質(zhì)運(yùn)動對聲傳播的影響，建立一種利用高斯波束法求解亞音速運(yùn)動介質(zhì)中的聲傳播模型，認(rèn)為在海流會明顯影響聲的傳播。

上述文獻(xiàn)研究中，多采用理想渦流速分布公式，甚少通過現(xiàn)代計算流體力學(xué)(Computational Fluid Dynamics,CFD)方法獲取渦流場流速參數(shù)。文獻(xiàn)[13]雖然采用射線聲學(xué)方法對運(yùn)動介質(zhì)中的聲傳播進(jìn)行研究，但沒有涉及復(fù)雜流場中的聲傳播問題。本文針對聲波通過水下穩(wěn)態(tài)渦流場的數(shù)值模擬方法及實(shí)驗展開研究，首先根據(jù)穩(wěn)態(tài)渦流場中的射線微分方程數(shù)值模擬聲線傳播軌跡，然后開展水池超聲傳感器實(shí)驗，最后將數(shù)值模擬和實(shí)驗測量的通過渦流場的聲信號時延量進(jìn)行對比驗證，證明了本文通過射線聲學(xué)方法數(shù)值模擬水下穩(wěn)態(tài)渦流場聲傳播的有效性。

1 渦流場聲線軌跡數(shù)值模擬

渦流場聲傳播特性的改變，主要表現(xiàn)為通過渦流場的聲波振幅與相位的變化[14]。渦流場聲傳播效應(yīng)示意圖如圖1所示，順流區(qū)聲速加速，逆流區(qū)聲速減速，最終導(dǎo)致在接收端聲波的相位超前或滯后。

圖1 渦流場聲傳播效應(yīng)示意圖Fig.1 Schematic diagram of sound propagation effect in vortex field

射線聲學(xué)是幾何聲學(xué)的近似理論，適用于解決高頻情況下的聲傳播問題。基于射線理論可以有效模擬聲波通過渦流場的聲線軌跡及信號變化，具有直觀、計算效率高的特點(diǎn)。為模擬聲波通過渦流場的聲線軌跡，根據(jù)穩(wěn)定移動介質(zhì)中的程函方程[15]，推導(dǎo)出可進(jìn)行聲線軌跡模擬計算的穩(wěn)態(tài)渦流場射線聲學(xué)微分方程組：

式中：γ(r)為渦度；α表示聲線掠射角；u表示渦流場速度。

聲線軌跡數(shù)值仿真示意圖如圖2所示。在理論計算過程中，對于某一條聲線s，固定每步迭代中x方向的增量Δx，通過式(1)獲取y方向的增量Δy，進(jìn)而確定該聲線的軌跡。

圖2 聲線軌跡數(shù)值仿真示意圖Fig.2 Schematic diagram of sound ray numerical simulation

聲線通過渦流場后的數(shù)值模擬軌跡如圖3 所示。由圖3可知，通過渦流場的聲線掠射角發(fā)生偏移，越靠近渦核中心其掠射角變化越大，并在渦核后方形成聲線焦散區(qū)。仿真中，渦流場切向速度(簡稱：流速)r分布采用二維穩(wěn)態(tài)蘭金渦流(Rankine vortex)：

圖3 a=10 mm，Ma=4.6×10-2時的聲線傳播軌跡Fig.3 Sound ray propagation trajectory when a=10 mm,Ma=4.6×10-2

式中：r表示渦流場某一測點(diǎn)到渦核中心點(diǎn)的距離，a為渦流場流速最大處渦核特征半徑；Γ=2πacMa，定義為渦流環(huán)量；c=1 500 m·s-1為水中聲傳播速度；Ma為馬赫數(shù)，定義為渦流最大切向速度與聲速之比。

本文數(shù)值模擬結(jié)果與文獻(xiàn)[16]實(shí)驗結(jié)果進(jìn)行對比，在相同條件設(shè)置下，聲信號相位差計算結(jié)果如圖4所示。由圖4可見，相位差的數(shù)值模擬結(jié)果與文獻(xiàn)實(shí)驗結(jié)果整體上吻合較好。為了進(jìn)一步驗證數(shù)值模擬方法，本文設(shè)計并搭建了超聲測量渦流場實(shí)驗平臺，開展相關(guān)實(shí)驗驗證。

圖4 聲信號相位差本文數(shù)值模擬結(jié)果與文獻(xiàn)[16]實(shí)驗結(jié)果的對比Fig.4 Comparison between the numerical simulation result in this paper and the experimental result from reference[16] regarding the phase difference of acoustic signal

2 實(shí)驗裝置及渦流場參數(shù)

2.1 實(shí)驗儀器設(shè)備

實(shí)驗設(shè)備實(shí)物圖如圖5(a)所示。實(shí)驗通過螺旋槳產(chǎn)生穩(wěn)態(tài)渦流場，螺旋槳直徑d=60 mm，固定在尺寸為2 m×0.6 m×0.8 m 的亞克力長方體水池中，實(shí)驗布置如圖5(b)所示。

圖5 實(shí)驗設(shè)備實(shí)物圖及實(shí)驗布置示意圖Fig.5 Photo of the experimental equipment and schematic diagram of experiment layout

螺旋槳位于整渦管頂部，由電機(jī)帶動旋轉(zhuǎn)產(chǎn)生渦流。渦流流出整渦管后，使用超聲傳感器發(fā)射和接收通過渦流的聲信號。為保證兩個超聲傳感器中軸在同一條直線上，超聲傳感器剛性固定在同一移動滑臺上。水池裝有導(dǎo)軌和滑尺，便于精確標(biāo)定傳感器位置。針對螺旋槳高速旋轉(zhuǎn)產(chǎn)生的振動，在實(shí)驗中將造渦裝置和測量裝置分別固定在兩組滑臺上，在螺旋槳與其固定滑臺聯(lián)接處使用螺絲螺母固定并加墊片以減少振動。儀器設(shè)備型號及參數(shù)如表1所示。

表1 實(shí)驗儀器設(shè)備型號及參數(shù)Table 1 Model and parameters of experimental equipment

2.2 渦流場CFD仿真

由于實(shí)驗中的渦流場速度分布不易獲取，本文基于現(xiàn)代流體仿真軟件Fluent，根據(jù)實(shí)驗布置對渦流場分布進(jìn)行數(shù)值仿真。本文使用Design Modeler建模軟件建立三維模型。計算域的劃分與邊界條件的設(shè)置如圖6所示。采用ICEM前處理軟件劃分網(wǎng)格，使用六面體結(jié)構(gòu)化網(wǎng)格對整體區(qū)域進(jìn)行劃分，并在旋轉(zhuǎn)域表面進(jìn)行局部加密，以更好地模擬近壁流動。

圖6 數(shù)值模擬計算域Fig.6 Numerical simulation computing domain

螺旋槳所在的旋轉(zhuǎn)域采用穩(wěn)態(tài)(Moving Refer‐ence Frame,MRF)方法，并繞z軸以1 400 r·min-1的轉(zhuǎn)速旋轉(zhuǎn)。基于有限體積法求解RANS方程，湍流模型選用標(biāo)準(zhǔn)k-ε模型，調(diào)用標(biāo)準(zhǔn)壁面函數(shù)進(jìn)行計算。采用三維雙精度基于壓力隱式求解器，耦合方式采用SIMPLE 算法，離散格式采用二階迎風(fēng)格式。

2.3 渦流場參數(shù)

定義整渦管出口所在平面為z=0平面，分別在整渦管出口下端z=10、30和50 mm處取截面，得到渦流場速度分布矢量圖如圖7所示，其中z=10 mm截面上的流速分布云圖如圖8所示。由圖7、8可以看到渦流場在流出整渦管后向外消散，符合自由場渦流運(yùn)動規(guī)律。

圖7 渦流場速度分布矢量圖Fig.7 Velocity distribution vector diagram of vortex field

圖8 在z=10 mm處的平面渦流場速度云圖Fig.8 Velocity cloud diagram of the plane vortex field at z=10 mm

渦流流速最大位置對應(yīng)渦核半徑尺寸，由圖8可知，z=10 mm平面上的最大速度umax=1.258 m·s-1，判斷渦核半徑a=30.6 mm。根據(jù)z=10、30和50 mm處的流場速度數(shù)據(jù)可得到不同截面的渦流場主要參數(shù)，如表2所示，為基于射線聲學(xué)的渦流場聲傳播數(shù)值模擬方法提供輸入?yún)?shù)。

表2 最大切向速度和渦核半徑Table 2 The maximum tangential velocity and vortex core radius obtained by simulation

3 實(shí)驗研究

3.1 實(shí)驗工況設(shè)置

為避免發(fā)生聲衍射，聲波波長需遠(yuǎn)小于渦流場尺度。根據(jù)表2中的參數(shù)，選擇聲波的發(fā)射頻率為2 MHz，幅度為3 V，示波器采樣率為1.0 GHz。實(shí)驗工況通過改變兩個傳感器位置獲取不同截面的聲信號時延。在整渦管出口位置下端10、30和50 mm處分別設(shè)置3個測量平面，如圖9(a)所示；在每個平面上，依次間隔10 mm設(shè)置a～k共11組測點(diǎn)，測點(diǎn)布置示意圖如圖9(b)所示。

圖9 測量平面正視及俯視示意圖Fig.9 Schematic diagram of the front and top views of the measuring plane

按照表3所示工況在每個測量平面上依次對11個測點(diǎn)進(jìn)行測量、記錄、存儲數(shù)據(jù)。每處測點(diǎn)測量5次。其中，工況1是無渦情況，只需測一組工況，即以z=10 mm 平面上的測點(diǎn)a測量無渦工況下聲時延。

表3 實(shí)驗工況表Table 3 Table of experimental conditions

3.2 數(shù)據(jù)處理與分析

針對有渦情況，分別獲取有渦與無渦工況下的超聲傳感器接收端的時域信號。其中一組接收信號的對比如圖10 所示。采用互相關(guān)算法即可得到此組數(shù)據(jù)的時延量為2.74×10-8s。

圖10 有渦和無渦工況下接收信號對比Fig.10 Comparison of the receiving signals in the vortex and no vortex conditions

不同截面的時延計算結(jié)果如圖11 所示。圖11中實(shí)線部分為根據(jù)實(shí)測數(shù)據(jù)擬合得到的時延曲線，虛線部分為針對表2中不同截面處的渦流場，采用基于射線聲學(xué)的數(shù)值模擬方法計算得到的聲信號時延值。由圖11可知，時延曲線呈現(xiàn)“S”形變化趨勢，說明渦流場對聲傳播有明顯影響。隨著z值增大，時延最大變化量減少，即渦流場最大切向速度變小；時延最大變化量對應(yīng)的x值增大，即渦流半徑增大，說明渦流在向下運(yùn)動的過程中不斷擴(kuò)散。

圖11 實(shí)驗測量時延數(shù)據(jù)及擬合曲線Fig.11 Time delay data and their fitting curves obtained by experimental measurement

采用相位重構(gòu)算法[9]，根據(jù)實(shí)測的時延擬合曲線反演出渦流場流速分布，并與CFD 仿真流速數(shù)據(jù)以及理想Rankine渦流速數(shù)據(jù)進(jìn)行對比分析，結(jié)果如圖12所示。由圖12可知，CFD仿真得到的流場速度分布與實(shí)測數(shù)據(jù)反演流場速度分布更為接近，而根據(jù)式(2)計算的Rankine渦流速分布在最大流速兩側(cè)速度衰減較慢，這是因為實(shí)驗環(huán)境與理想流場假設(shè)存在差異的。而基于CFD 計算渦流場數(shù)據(jù)，條件設(shè)置與實(shí)驗條件一致，相較于式(2)更能反映出真實(shí)的流場分布。

圖12 在z=10 mm處的三組平面流速分布數(shù)據(jù)對比圖Fig.12 Comparison of three sets of flow velocity distribution data in the plane at z=10 mm

以數(shù)值模擬的時延為真值，采用均方根誤差公式，計算圖11 中實(shí)測時延擬合曲線和數(shù)值模擬曲線之間的誤差。由于數(shù)值模擬中所代入的流場分布是理想Rankine渦流速，在最大流速兩側(cè)與實(shí)際流場分布有一定差距，因此只計算最大流速中間對應(yīng)的時延曲線誤差，計算結(jié)果如表4 所示。由表4可知，數(shù)值模擬計算的時延與實(shí)驗測量值基本吻合，誤差均小于10%，說明本文提出數(shù)值模擬方法的有效性。

表4 時延數(shù)據(jù)與誤差分析Table 4 Time delay data and error analysis

4 誤差分析

本實(shí)驗主要針對聲時延量進(jìn)行測量，該時延量為小量，為了進(jìn)一步分析實(shí)驗結(jié)果的可信性，下面針對影響時延測量的誤差因素進(jìn)行討論。

4.1 傳感器指向性

實(shí)驗中所使用的超聲傳感器具有指向性，主瓣-3 dB處的開角Φ為4°±1°，可通過改變發(fā)射端與接收端軸向偏移量研究指向性對時延測量的影響。

在整渦管出口下方z=10 mm 平面按照圖13 所示進(jìn)行傳感器的布置，測點(diǎn)a處發(fā)射端與接收端中軸嚴(yán)格對齊，隨后依次偏移5 mm 的測點(diǎn)設(shè)為b～e測點(diǎn)，在各測點(diǎn)處測量有渦和無渦條件下的聲信號。由于超聲傳感器指向性的影響，不同測點(diǎn)處的聲信號幅值發(fā)生改變。改變接收信噪比(Signal to Noise Ratio,SNR)，每個測點(diǎn)的信噪比測量結(jié)果如表5所示。

表5 不同測點(diǎn)的信噪比測量結(jié)果Table 5 Measurement results of SNR at different positions

圖13 傳感器布放位置示意圖Fig.13 Schematic diagram of the sensor placement

信噪比的改變將直接影響時延估計的精度。根據(jù)本文實(shí)驗條件，計算互相關(guān)時延估計算法在不同信噪比條件下的時延誤差如圖14所示。由圖14可知，當(dāng)信噪比大于10 dB時，時延估計誤差均小于5%。由表5可知，實(shí)驗時發(fā)射端與接收端中軸偏移量小于10 mm就可以保證10 dB以上的信噪比。本實(shí)驗中兩個超聲傳感器被剛性固定在同一移動滑臺上，中軸偏移量遠(yuǎn)小于10 mm，信噪比較高，因此超聲傳感器指向性對實(shí)驗結(jié)果影響較小。

4.2 整渦管

在z=10 mm平面上，整渦管的存在可能會對聲波產(chǎn)生反射、透射和衍射等產(chǎn)生影響。整渦管出口干擾示意圖如圖15 所示。由于本實(shí)驗所使用的傳感器發(fā)射和接收均具有較窄的指向性，則整渦管的存在對聲信號所產(chǎn)生的影響較小，下面通過實(shí)測數(shù)據(jù)分析整渦管對本實(shí)驗聲傳播的影響。

圖15 整渦管出口干擾示意圖Fig.15 Schematic diagram of vortex tube outlet interference

在無渦條件下，在z=10 mm平面處，對存在整渦管和無整渦管兩種情況，分別測量聲波從發(fā)射端到接收端的傳播時間。兩種情況下傳播時間相對誤差為2.27%，即整渦管對z=10 mm平面上的聲時延測量影響較小，其他平面離整渦管距離更遠(yuǎn)，則影響將更小。

5 結(jié) 論

本文圍繞水下穩(wěn)態(tài)渦流場的聲傳播展開研究，給出了穩(wěn)態(tài)渦流場中的聲線傳播的微分方程，通過基于射線聲學(xué)的數(shù)值模擬方法計算聲波通過渦流場的時延，并與實(shí)驗測量結(jié)果進(jìn)行對比分析，得到如下結(jié)論：

(1) 基于射線聲學(xué)的數(shù)值模擬方法能有效計算聲波通過渦流場后的時延，數(shù)值模擬誤差小于10%。

(2) 通過現(xiàn)代CFD方法獲取水下穩(wěn)態(tài)渦流場的流速分布，較理想渦流場流速的理論公式更接近實(shí)際渦流場聲速的分布。

(3) 通過誤差分析可知，本實(shí)驗中由傳感器指向性和整渦管產(chǎn)生的誤差較小，驗證了本實(shí)驗測量結(jié)果的可靠性。