楊德森，李中政，2，方爾正

(1.哈爾濱工程大學(xué)水聲工程學(xué)院，黑龍江哈爾濱150001;2.海軍92198部隊(duì)，遼寧興城125109)

?

三列不同相位原波形成的參量陣聲場研究

楊德森1，李中政1，2，方爾正1

(1.哈爾濱工程大學(xué)水聲工程學(xué)院，黑龍江哈爾濱150001;2.海軍92198部隊(duì)，遼寧興城125109)

摘要:參量陣能生成低頻、高指向性聲束，在海底地層剖面測量、掩埋目標(biāo)探測與識(shí)別、水下通信等水聲工程領(lǐng)域有著廣泛的應(yīng)用。理論研究表明，當(dāng)三列原波相互作用時(shí)，能形成具有兩個(gè)差頻分量的參量陣輻射系統(tǒng)，但尚未有文獻(xiàn)對(duì)受原波相位影響的多頻參量陣聲場進(jìn)行系統(tǒng)分析?；诖耍疚木C合計(jì)算精度和效率的需求，借鑒算子分離思想，利用時(shí)頻域相結(jié)合的方法求解KZK方程，實(shí)現(xiàn)對(duì)三列不同相位原波相互作用形成的參量陣聲場的描述。計(jì)算結(jié)果表明:參量陣差頻信號(hào)的聲壓幅值隨著原波相位的改變而改變，并且原波相位均為零時(shí)，兩個(gè)差頻分量的聲壓幅值最大。最后進(jìn)行了試驗(yàn)研究，實(shí)驗(yàn)結(jié)果與仿真結(jié)果吻合較好。

關(guān)鍵詞:參量陣;三波作用;原波相位;算子分離法;時(shí)頻域相結(jié)合

1962年，Westervelt［1］提出了參量陣的最初模型:同向傳播的兩列高頻原波利用媒質(zhì)的非線性效應(yīng)獲得差頻波的聲發(fā)射裝置。1965年，Berktay［2］引入調(diào)制“包絡(luò)”的概念，給出了聲學(xué)參量陣更精確和完整的理論解釋，使得聲源輻射信號(hào)不再局限于雙頻基波的情況。Lucilla Di Marcoberardin［3］，John A.Birken［4］對(duì)多列聲波輻射形成系列差頻分量的參量陣聲場進(jìn)行了理論及實(shí)驗(yàn)研究，為海底地質(zhì)特性、管道或電纜的鋪設(shè)提供技術(shù)指導(dǎo)。2002年，王潤田等［5］用有多個(gè)差頻檔位的“堤防隱患監(jiān)測聲吶”對(duì)堤防的損毀程度進(jìn)行探測和評(píng)估。顯然，多頻原波相互作用形成的參量陣在水聲工程領(lǐng)域有非常廣泛的應(yīng)用前景。

在眾多非線性聲場計(jì)算的理論模型中，KZK方程［6－8］在計(jì)算參量陣近場聲場時(shí)有著非常明顯的優(yōu)勢。Masahiko Akiyama等［9－11］利用KZK方程對(duì)不同相位雙頻原波形成的參量陣聲場特性進(jìn)行計(jì)算，并進(jìn)行了試驗(yàn)研究。Fenlon等［12］基于Burgers方程推導(dǎo)了多頻有限振幅聲波的Bessel－Fubini解，呂君等［13］根據(jù)Fenlon理論研究了二階近似條件下的雙頻聲源相互作用時(shí)的遠(yuǎn)場指向性。為研究有限振幅波在介質(zhì)中傳播的頻域能量變化問題，楊德森等［14］基于頻譜分解方法推導(dǎo)了單頻及多頻大振幅的傳播規(guī)律。

到目前為止，并沒有相關(guān)文獻(xiàn)對(duì)三列聲波相互作用形成的參量陣聲場特性進(jìn)行研究。本文借鑒算子分離思想［8］，綜合考慮計(jì)算效率和精度，將KZK方程分解為線性(包含衍射、吸收效應(yīng))和非線性(包含非線性效應(yīng))兩部分，線性部分采用二階龍哥庫塔有限差分法(DIRK2)和Crank－Nicolson有限差分法(CNFD)相結(jié)合［7，18］的方法;非線性部分則采用守恒型迎風(fēng)格式［15－17］積分求解。最后將兩部分進(jìn)行整合，對(duì)不同相位三列原波相互作用形成的參量陣聲場進(jìn)行分析，旨在為參量陣的進(jìn)一步工程應(yīng)用提供相應(yīng)的理論指導(dǎo)。

1　KZK方程的相關(guān)理論

綜合考慮有限振幅聲波傳播過程中的吸收、衍射和非線性效應(yīng)，建立適用于參量陣聲場計(jì)算的KZK方程［6－7］:

式中:gn、hn是第n次諧波聲壓在空間位置(rj，zm)處的諧波系數(shù)。

1.1二階龍格庫塔法計(jì)算KZK方程的線性部分

利用算子分離思想對(duì)式(1)的吸收、衍射部分進(jìn)行頻域計(jì)算，得到如下表達(dá)式:

用IBFD法進(jìn)行有限差分離散，化為矩陣有

式中:AJ×J=三對(duì)角矩陣，

將式(3)進(jìn)行改寫:

所以，利用IBFD法對(duì)KZK方程的線性部分進(jìn)行離散得到的局部截?cái)嗾`差為O(dz)。為了提高計(jì)算精度，與CNFD法的局部截?cái)嗾`差相匹配，根據(jù)二階龍格庫塔思想［17］，將既定的軸向步長dz分為b1dz、b2dz(b1+b2=1)兩部分，分別計(jì)算zm，zm+b1Δz、處的斜率s1、s2:

1.2守恒迎風(fēng)格式計(jì)算KZK方程的非線性部分

采用頻域方法對(duì)KZK方程進(jìn)行數(shù)值計(jì)算，非線性部分形象、直觀地體現(xiàn)了各階諧波的耦合，但計(jì)算時(shí)間正比于N2(N為最大截?cái)嘀C波階次，衍射、吸收項(xiàng)的計(jì)算時(shí)間正比于N)，故計(jì)算量非常大。參量陣的形成過程屬于弱非線性聲學(xué)范疇，其非線性系數(shù)較小，可以將頻域信號(hào)變換為時(shí)域波形，換取較大計(jì)算步長達(dá)到提高計(jì)算效率的目的。

1.1節(jié)的求解得到zm+1處聲壓幅值p(rj，zm+1)是考慮聲波非線性傳播衍射、吸收效應(yīng)，接下來將其作為節(jié)點(diǎn)(rj，zm)處的初始條件對(duì)KZK方程的非線性部分(無粘滯Burgers方程)求解［15－17］:

方程(7)兩邊同時(shí)對(duì)時(shí)間進(jìn)行積分，將聲壓轉(zhuǎn)換到時(shí)域，有如下表達(dá)式:

對(duì)無粘滯Burgers方程進(jìn)行守恒型迎風(fēng)格式表示:

由式(10)、(11)就可以計(jì)算空間層zm→zm+1的時(shí)域波形了，然后將其通過FFT變換到頻域，以其作為初始條件代入KZK方程的線性部分，求解軸向位置zm+1→zm+2的聲場傳播。依次類推，即可獲知整個(gè)求解區(qū)域的聲場空間分布特性。

2　仿真結(jié)果分析

按照上述KZK方程求解參量陣聲場的方案，設(shè)定其計(jì)算區(qū)域的網(wǎng)格劃分情況如下，軸向計(jì)算范圍為zmax=3.5r0，徑向計(jì)算范圍是rmax=41a，考慮到計(jì)算的方便，在頻域計(jì)算過程中，對(duì)諧波階數(shù)進(jìn)行截?cái)啵x定N=128。

圖1　參量陣聲場空間的網(wǎng)格劃分示意圖Fig.1 Mesh grids of the parametric array’s acoustic field

另外，為了盡量降低邊界的反射，本文借鑒求解電磁場的PML吸收邊界條件，將徑向網(wǎng)格的外圍邊界進(jìn)行處理，設(shè)定1倍半徑長度的邊界層作為PML邊界［15］(如圖1陰影部分所示)。

選定如下仿真參數(shù):水的密度ρ=1 000 kg/m3、聲速c=1 482 m/s，水的非線性系數(shù)為β=3.5，圓形活塞聲源的表面峰值聲壓為P0=60 kPa，并且隨著相位的改變而變化。半徑為a=10 cm，水中聲波衰減系數(shù)與頻率的關(guān)系為:α/f2=2.17×10－13dB/m;聲源輻射的聲波:

2.1同相位原波形成的參量陣聲場空間分布

本節(jié)設(shè)定三列原波相位滿足θf1=θf2=θf3=0，對(duì)參量陣各個(gè)差頻分量的聲場特性進(jìn)行研究。

由圖2可知:

1)差頻分量的聲能量累積過程大致相同，都存在一個(gè)拐點(diǎn)。原波在近場相互作用時(shí)，非線性效應(yīng)比吸收效應(yīng)明顯，差頻波的聲能量累積速率較快;到達(dá)拐點(diǎn)之后，非線性效應(yīng)的聲能量累積速率小于吸收效應(yīng)引起的聲衰減速率，聲吸收效應(yīng)在參量陣的軸向聲場傳播過程中占據(jù)主導(dǎo)地位，參量陣差頻波聲壓幅值越來越小;

2)三列原波相互作用形成兩個(gè)具有高指向性的參量陣差頻波束，并且兩個(gè)差頻分量的波束寬度基本一致，但差頻分量fd2=8 kHz的軸向聲壓幅值比fd1=8 kHz的聲壓幅值大。

圖2　參量陣的聲壓幅值分布特性Fig.2 The amplitude distribution of the parametric array’s pressure

2.2差頻分量fd1的聲壓幅值分布特性

對(duì)不同相位原波相互作用形成的參量陣差頻分量fd1=4 kHz的聲場進(jìn)行分析，給出如下組圖3，由圖3可知:

1)當(dāng)θf2=θf3=0時(shí)，隨著原波f1=41 kHz相位θf1的增加(如圖3(a)的上圖所示)，同一軸向位置處差頻分量fd1=4 kHz的參量陣聲壓幅值越來越低，當(dāng)θf1=θf3=0時(shí)，隨著原波f2=45 kHz相位θf2的增加(如圖3(a)的下圖所示)，同一軸向位置處差頻分量fd1=4 kHz的參量陣聲壓幅值越來越低，當(dāng)θf1=θf3=0時(shí)，隨著原波f3=49 kHz相位θf3的增加(如圖3(b)的上圖所示)，同一軸向位置處差頻分量fd1=4 kHz的參量陣聲壓幅值越來越低。這表明三列原波相互作用時(shí)，原波疊加形成的調(diào)制信號(hào)包絡(luò)隨著相位的改變而改變，進(jìn)而影響到差頻信號(hào)fd1的軸向聲壓幅值分布;

2)相較原波f1=41 kHz和f3=49 kHz相位改變引發(fā)的差頻分量fd1=41 kHz軸向聲壓幅值變化，原波f2=45 kHz相位改變對(duì)參量陣差頻分量fd1=41 kHz的影響更為明顯，這是因?yàn)椴铑l分量fd1的形成過程為fd1=f3－f2，fd1=f2－f1，顯然，其非線性作用機(jī)制較為復(fù)雜。

圖3　原波相位改變時(shí)，參量陣聲壓幅值變化規(guī)律Fig.3 The parametric array’s pressure amplitude characteristics due to the primary wave’s phase changed

為了對(duì)相同相位的原波相互作用形成的參量陣聲場進(jìn)行分析，給出仿真組圖4。

圖4　參量陣差頻分量fd1的聲壓幅值特性Fig.4　The pressure amplitude characteristic of the difference frequency component fd1

由圖4可知:

1)三列原波的相位滿足θf1=θf2=θf3=0時(shí)，軸向聲壓幅值最大，并且原波f2=45 kHz的相位θf2改變時(shí)，同一軸向位置處參量陣差頻分量fd1的聲壓幅值最低;

2)原波f1=41 kHz和f3=49 kHz具有相同的相位時(shí)，形成的參量陣差頻分量fd1聲壓幅值相等。說明非載波原波的相位變化對(duì)參量陣聲場空間分布具有相同的影響。

2.3差頻分量fd2的聲壓幅值分布特性

本節(jié)對(duì)不同相位原波相互作用形成的參量陣差頻分量fd2的聲壓幅值分布特性進(jìn)行研究，給出如下仿真結(jié)果。

圖5　原波相位改變形成的參量陣聲壓幅值變化特性Fig.5 The parametric array’s pressure amplitude characteristics due to the phase changed of the primary wave

相對(duì)fd1=41 kHz來講，參量陣差頻分量fd2=8 kHz的聲壓幅值分布受原波相位的影響較小，這是因?yàn)椴铑l分量fd2的形成過程是fd2=f3－f1，其非線性作用機(jī)制更為簡單。

另外，fd2的聲壓幅值有小范圍的波動(dòng)，并且同一軸向位置的聲壓幅值隨著相位的增加而降低，這是由于三列原波疊加形成的調(diào)制信號(hào)峰值降低導(dǎo)致參量陣聲源輻射系統(tǒng)的聲源級(jí)降低。如果忽略聲源級(jí)帶來的影響，那么就可認(rèn)為兩列不同相位原波形成的參量陣差頻波聲壓幅值分布受相位的影響較小，跟文獻(xiàn)［9－12］的研究結(jié)果相吻合。

同樣，為了對(duì)相同相位的不同原波相互作用形成的參量陣聲場進(jìn)行分析，給出仿真組圖6。

圖6　差頻分量fd2的聲壓幅值特性Fig.6　The pressure amplitude characteristic of the difference frequency component fd2

圖6表明，原波相位相同時(shí)，同一徑向位置的差頻分量幅值最大。隨著原波相位的增加，參量陣差頻分量fd2的聲壓幅值越來越低，并且原波f1，f3的相位相同時(shí)，聲壓幅值的分布相同，亦即非中心頻率分量的原波相位改變具有相同的功效。

3　試驗(yàn)研究

為了對(duì)不同相位原波形成的參量陣聲場進(jìn)行研究，本文在哈爾濱工程大學(xué)水聲工程學(xué)院的消聲水池對(duì)軸向聲壓幅值進(jìn)行測量的實(shí)驗(yàn)。

圖7　參量陣軸向聲壓幅值仿真值與實(shí)驗(yàn)值的比較=60°Fig.7 Comparison of the experimental value of the parametric array with the simulation values when=60°

選定原波f1、f2、f3的相位分別為=90°、==0和=60°、==0°兩組工況，對(duì)仿真結(jié)果和實(shí)驗(yàn)結(jié)果進(jìn)行比對(duì)。

對(duì)圖7進(jìn)行分析可知，三列不同相位的原波相互作用可以生成具有兩個(gè)高指向性、差頻分量的參量陣輻射系統(tǒng)。另外，由于KZK方程在數(shù)值計(jì)算參量陣聲場空間分布時(shí)，僅考慮傳播媒介的粘滯性吸收，導(dǎo)致遠(yuǎn)場聲場由于吸收效應(yīng)引起的聲衰減小于實(shí)驗(yàn)測量值。

圖8給出參量陣差頻分量fd1的聲壓幅值隨軸向距離的變化圖示，其中圖8(a)表示f1=41 kHz分別為30°、90°，θf1=θf3=0時(shí)的聲壓幅值軸向分布圖示，圖8(b)表示f1=45 kHz，θf2分別為0°、60°，θf1=θf3=0時(shí)的聲壓幅值軸向分布圖示。

圖8　參量陣軸向聲壓幅值仿真值與實(shí)驗(yàn)值比較Fig.8　Comparison of the experimental value of the parametric array with the simulation values of axial sound pressure amplitude

圖9　參量陣軸向聲壓幅值仿真值與實(shí)驗(yàn)值比較Fig.9 Comparison of the parametric array’s experimental value with the simulation values

圖9表明在誤差允許范圍內(nèi)，參量陣差頻分量fd1隨著軸向距離的增加逐漸降低，跟仿真結(jié)果吻合程度較高。

接下來選定原波頻率為f3=49 kHz、其相位θf3=0°、30°、60°，θf1=θf2=0時(shí)的參量差頻分量fd2的軸向聲壓幅值分布。

對(duì)實(shí)驗(yàn)結(jié)果分析可知，差頻分量fd2隨軸向距離的變化規(guī)律跟KZK方程仿真計(jì)算結(jié)果基本一致，并且在近場區(qū)域其聲壓幅值變化規(guī)律更為接近，亦即KZK方程描述參量陣的近場聲場優(yōu)勢較為明顯。

4　結(jié)論

本文借鑒KZK方程的時(shí)域算子分離思想，將時(shí)域有限差分與頻域有限差分相結(jié)合，把求解參量陣非線性聲場的KZK方程分為兩部分:采用守恒迎風(fēng)格式對(duì)KZK方程的非線性部分進(jìn)行計(jì)算，其中由于得到的線性部分采用二階龍格庫塔法對(duì)KZK方程的衍射、吸收效應(yīng)進(jìn)行計(jì)算，最后進(jìn)行整合，得到三列不同相位原波相互作用形成的參量陣聲場分布特性，通過研究得出如下結(jié)論:

1)DIRK2是具有二階計(jì)算精度的有限差分法，該算法可提高計(jì)算精度，并且能跟遠(yuǎn)場計(jì)算的Crank－Nicolson有限差分法計(jì)算精度相匹配;

2)相對(duì)頻域有限差分計(jì)算KZK方程，采用守恒迎風(fēng)格式變換到時(shí)域求解參量陣聲場的非線性部分可有效提高計(jì)算效率，降低計(jì)算量;

3)當(dāng)三列原波相互作用時(shí)，能獲得兩個(gè)差頻分量的高指向性聲束，并且三列原波的相位均為0時(shí)，其軸向聲壓幅值取最大值;

4)三列不同相位的原波相互作用時(shí)，頻率較低的差頻分量fd1由于其非線性作用機(jī)制較為復(fù)雜，其聲壓幅值隨相位改變有明顯的振蕩。與此同時(shí)，頻率較高的差頻分量fd2的聲壓幅值隨相位的改變僅有小范圍的波動(dòng)。

參考文獻(xiàn):

［1］WESTERVELT P J.Parametric acoustic array［J］.The Journal of the Acoustical Society of America，1963，35(4):535－537.

［2］BERKTAY H O.Possible exploitation of non－linear acoustics in underwater transmitting applications［J］.Journal of Sound and Vibration，1965，2(4):435－461.

［3］DI MARCOBERARDINO L，MARCHAL J，CERVENKA P.Nonlinear multi－frequency transmitter for sea－floor characterization［C］//Proceedings of the Acoustics Conference.Nantes，F(xiàn)rance，2012:2800－2805.

［4］BIRKEN J A.Empirical results from frequency－scanning nonlinear sonar in deep water［J］.The Journal of the A－coustical Society of America，1974，56(S1):S41.

［5］王潤田.海底聲學(xué)探測與底質(zhì)識(shí)別技術(shù)的新進(jìn)展［J］.聲學(xué)技術(shù)，2002，21(1/2):96－98.WANG Runtian.Progress in detecting the geological formations and sediment properties by sound［J］.Technical A－coustics，2002，21(1/2):96－98.

［6］KAMAKURA T，HAMADA N，AOKI K，et al.Nonlinearly generated spectral components in the near－field of a directive sound source［J］.The Journal of the Acoustical Society of America，1989，85(6):2331－2337.

［7］NAZE TJ?TTA J，TJ?TTA S，VEFRING E H.Propagation and interaction of two collinear finite amplitude sound beams ［J］.The Journal of the Acoustical Society of America，1990，88(6):2859－2870.

［8］LEE Y S，HAMILTON M F.Time－domain modeling of pulsed finite－amplitude sound beams［J］.The Journal of the Acoustical Society of America，1995，97(2):906－917.

［9］KAMAKURA T，SAKAI S，NOMURA H，et al.Parametric audible sounds fields by phase－cancellation excitation of primary waves［J］.The Journal of the Acoustical Society of A－merica，2008，123(5):3694.

［10］KAMAKURA T，NOMURA H，SAKAI S.Spatial phaseinversion technique for parametric source with suppressed carrier［J］.The Journal of the Acoustical Society of America，2009，125(4):2717.

［11］KAMAKURA T，NOMURA H，AKIYAMA M，et al.Parametric sound fields formed by phase－inversion excitation of primary waves［J］.Acta Acustica United with Acustica，2011，97(2):209－218.

［12］FENLON F H.An extension of the Bessel－Fubini series for a mutihle－frequency CW acoustic source of finite amplitude ［J］.The Journal of the Acoustical Society of America，1972，51(1B):284－289.

［13］呂君，趙正予，周晨，等.有限振幅聲波間的非線性相互作用對(duì)聲源遠(yuǎn)場指向性的影響［J］.物理學(xué)報(bào)，2011，60(8):084301.LV Jun，ZHAO Zhengyu，ZHOU Chen，et al.Effect of finite－amplitude acoustic wave nonlinear interaction on farfield directivity of sound source［J］.Acta Physica Sinica，2011，60(8):084301.

［14］楊德森，董雷，時(shí)潔，等.大振幅波非線性傳播的頻率特性［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2010，31(7):928－931.YANG Desen，DONG Lei，SHI Jie，et al.Frequency characteristics of nonlinear propagation of large amplitude waves［J］.Journal of Harbin Engineering University，2010，31(7):928－931.

［15］SONESON J E.A parametric study of error in the parabolic approximation of focused axisymmetric ultrasound beams ［J］.The Journal of the Acoustical Society of America，2012，131(6):EL481－EL486.

［16］SONESON J E.A user－friendly software package for HIFU simulation［C］//Proceedings of the 8th International Symposium on Therapeutic Ultrasound.Minneapolis，Minnesota:AIP，2009，1113:165－169.

［17］PINTON G F.Numerical methods for nonlinear wave propagation in ultrasound［D］.Durham:Department of Biomedical Engineering，Duke University，2007:37－46.

［18］LI Zhongzheng，F(xiàn)ANG Erzheng，SHI Shengguo.Study on the calculation method of the KZK equation for parametric array［C］//Processings of the 12th International Conference on Signal Processing(ICSP).Hangzhou，China:IEEE，2014:139－143.

Study of the parametric array acoustic field yielded by three collimated primary waves with different phases

YANG Desen1，LI Zhongzheng1，2，F(xiàn)ANG Erzheng1
(1.College of Underwater Acoustics Engineering，Harbin Engineering University，Harbin 150001，China;2.Unit 92198 of Naval，Xingcheng 125109，China)

Abstract:A parametric array can generate an acoustic beam of low frequency and high directivity.This has been widely used in acoustics engineering including seafloor profile measurement，buried underwater target detection and recognition，and underwater communications.Theoretical studies have shown that the interaction of three collimated primary waves can produce a parametric array radiation system with two difference frequencies.However，no systematic analysis has been done of a multi－frequency parametric array acoustic field affected by the primary wave phase.To ensure the precision and calculation efficiency of the computation，this study drew on ideas from operator separation，combining the time domain with the frequency domain to solve the KZK equation for the parametric array，thereby achieving a description of a parameter array acoustic field formed by the interaction of three collimated primary waves with different phases.The simulation results showed that the sound pressure amplitude of the parametric array changed with a change in the primary wave phases.The sound pressure amplitude of the parametric array difference frequency component was largest when the primary wave phase was zero.Finally，the experimental study and the simulation produced identical results.

Keywords:parametric array;three collimated primary waves;primary wave phase;operator splitting scheme;time－frequency domain

通信作者:楊德森，E－mail:dshyang@ hrbeu.edu.cn.

作者簡介:楊德森(1957－)，男，中國工程院院士，教授，博士生導(dǎo)師.

基金項(xiàng)目:水聲技術(shù)國防科技重點(diǎn)實(shí)驗(yàn)室基金資助項(xiàng)目(9140c2001 0613c20078).

收稿日期:2014－10－28.網(wǎng)絡(luò)出版時(shí)間:2015－12－21.

中圖分類號(hào):O427.9

文獻(xiàn)標(biāo)志碼:A

文章編號(hào):1006－7043(2016)01－0007－07

doi:10.11990/jheu.201410082

網(wǎng)絡(luò)出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20151221.1555.022.html