劉 松 黎 勝 趙德有

1大連理工大學 工業(yè)裝備結構分析國家重點實驗室，遼寧 大連116024 2大連理工大學 運載工程與力學學部 船舶工程學院，遼寧 大連116024

基于相位共軛方法識別結構表面法向振速

劉 松1，2黎 勝1，2趙德有2

1大連理工大學 工業(yè)裝備結構分析國家重點實驗室，遼寧 大連116024 2大連理工大學 運載工程與力學學部 船舶工程學院，遼寧 大連116024

基于相位共軛方法對平板和水下圓柱殼輻射聲場的識別進行了數值仿真計算，得到結構表面聲壓分布后，通過兩種方法識別結構的法向振速：一種方法給出結構在聲場中的位置和尺寸，基于聲壓梯度計算法向振速；另一種方法根據結構表面阻抗關系計算法向振速。數值計算結果表明：基于聲壓梯度計算法向振速時，能得到法向振速幅值的大致分布，該方法計算簡單，適用于無法得到聲源表面阻抗關系的情況；而引入結構表面阻抗關系則能得到更加準確的識別結果。

相位共軛；法向振速；輻射聲場識別；聲壓梯度；阻抗

1 引言

噪聲源識別和定位是噪聲控制的關鍵問題，具有重要的應用價值和理論意義。目前聲場識別應用最廣泛的是近場聲全息方法 （Near－field A-coustic Holography，NAH）［1］。時域中的時間反轉（Time Reversal，TR）等價于頻域中的相位共軛（Phase Conjugation，PC）［2］。本文采用相位共軛（時間反轉）方法識別平板和水下圓柱殼的輻射聲場以及表面振動速度。時間反轉方法可實現聲波的反向傳播和自適應聚焦，可用于聲源定位和識別，在測得噪聲源正向傳播的聲場后，可以基于特定的時間反轉或相位共軛算法實現噪聲源的識別和聲場重構。時間反轉法用于聲源識別和定位的能力受到聲波衍射極限分辨率的限制，理論上其分辨率最高為二分之一波長大小，因此并不能提供出足夠的噪聲源位置和噪聲輻射特性信息，一些學者提出了解決時間反轉法中衍射極限的方法，如由Fink等提出的通過在近場使用亞波長散射體將倏逝波轉換為傳播波；Lerosey等［3］利用聲匯（acoustic sink）概念使用時間反轉鏡（單陣元）對超聲波實現了1／14波長的分辨率；Conti等［4］將近場傳聲器陣列和4個波長外的時間反轉鏡相結合實現了1／20波長的分辨率。Rosny等［5-8］采用兩個無限大平面陣列利用解析方法對點聲源的研究表明，基于聲壓梯度測量和偶極子源的近場時間反轉能突破衍射極限分辨率。本文基于相位共軛方法對平板和水下圓柱殼輻射聲場的識別進行了數值仿真計算，得到結構表面聲壓分布后，通過兩種方法識別結構的法向振速并給出相關誤差分析。

2 理論

2.1 結構聲輻射理論

結構在簡諧力作用下考慮流體加載效應的有限元形式的運動方程為［9］：

式中，［Zs］＝（－ω2［M］＋iω［C］＋［K］）／iω為結構阻抗矩陣，［M］、［C］和［K］分別為結構質量矩陣、阻尼矩陣和剛度矩陣，ω為激勵圓頻率；｛v｝為結構速度向量；｛fe｝為外激勵力向量；｛fp｝為結構表面聲壓所引起的流體對結構的作用力向量。

｛fp｝和結構表面聲壓｛p｝的關系式為：

式中，［G］cos為方向余弦轉換矩陣；［A］＝∫S［N］T［N］dS，［N］為形狀函數矩陣。

對從聲源向外輻射正向傳播的聲場，根據Helmholtz－Kirchhoff積分公式，聲源表面S的聲壓p（Q）和聲壓梯度（或法向振速vn（Q））的關系為：

式中，［Z］為阻抗矩陣；｛vn｝為結構表面法向速度向量。

由｛vn｝和結構速度向量｛v｝之間的轉換關系式：

可得求解結構－聲耦合問題的方程為：

求出結構速度向量｛v｝后，進而可求出｛vn｝和｛p｝等。

對簡諧激勵力作用下位于無限障板上的板結構在板一側半無限域流體介質和板表面Sp上產生的輻射聲壓p（P）也可由Rayleigh積分求得［10］：

同樣對板表面Rayleigh積分方程 （P∈Sp）進行離散，可得邊界元求解方程（4）。

2.2 相位共軛理論

對反向傳播的聲場，根據Helmholtz－Kirchhoff積分公式，在相位共軛陣列測量到聲壓和聲壓梯度后，其相位共軛聲場為：

式中，S′為閉合陣列表面。

實際的相位共軛陣列都是離散的，對包含N個陣元的離散的有限陣列，其相位共軛聲場為：

如基于測量聲壓使用單極子源來進行時間反轉，其相位共軛聲場為：

如基于測量聲壓梯度使用偶極子源來進行時間反轉，其相位共軛聲場為：

從公式（9）～（11）可以看到，以上基于相位共軛方法通過計算反向傳播的相位共軛聲場來進行聲源定位和識別具有實施簡單，可直接使用平面陣列等來進行局部重建，不存在解不穩(wěn)定性問題等優(yōu)點。本文數值模擬計算均基于聲壓梯度測量使用偶極子源（PCD）識別結構表面聲壓。

2.3 結構法向振速識別

得到識別的結構表面聲壓分布后，利用以下兩種方法識別結構的法向振速。

方法一：已知結構在聲場中所處的位置及尺寸，以及離散的節(jié)點和單元數，基于聲壓梯度按照式（12）識別結構法向振速。

式中，p1為結構在聲場中所在位置識別的聲壓值，p2是與p1對應的沿結構表面法線方向增加距離為Δ1時識別的聲壓值。

方法二：除了考慮結構在聲場中的位置和尺寸之外，還可以通過引入結構表面阻抗關系對法向振速進行識別。將識別的結構表面聲壓代入公式（4），并對阻抗矩陣求逆，得到法向振速，如式（13）所示。

3 數值模擬

3.1 平板輻射聲場和法向振速識別

計算采用的板結構參數為：板長0.8 m，板寬0.6 m，厚度0.005 m；材料參數為：E＝210 GPa，ρ＝7 833 kg／m3，v＝0.3；板幾何中心取為坐標原點，在板的 （－0.2，－0.1，0）處加單位幅值的簡諧激勵力，計算頻率f＝500 Hz，對應板的（5，1）模態(tài)；板邊界條件取為四邊簡支，并嵌入無限大剛性障板中；流體取為空氣，聲速為350 m／s，密度為1.21 kg／m3。

相位共軛測量陣列信息為：陣元距離板表面d＝0.1λ＝0.07 m，陣元個數221，陣元間距Δ＝0.05 m＝0.07λ，與板形狀一致。

本節(jié)數值模擬計算將按照Rayleigh積分法得到的平板在相位共軛陣列處的輻射聲壓值的共軛作為測量值，代入式（11）得到識別的平板表面聲壓幅值分布。其中相位共軛陣列上的聲壓可以用沿板法線方向相近位置兩次測量的聲壓值的算術平均，即代替；聲壓梯度可以取為?p?n＝（p2′－p1′）／Δ1，Δ1為相近兩次測量的距離，算例中取為0.001 m。

3.1.1 板表面聲壓的識別

采用Rayleigh積分計算和按照公式（11）基于聲壓梯度測量使用偶極子源（PCD）方法識別的平板表面聲壓分布，如圖1和圖2所示（聲壓單位：Pa；法向速度單位：m／s）。

圖1 Rayleigh積分計算的結構表面聲壓Fig.1 The surface pressure calculated by Rayleigh integral

圖2 PCD識別的表面聲壓結果Fig.2 The surface pressure reconstructed by PCD

從圖2的結果可以看出，基于聲壓梯度測量并采用偶極子源得到的相位共軛聲場幅值與通過Rayleigh積分得到的結果吻合較好，能夠突破聲波的衍射極限，準確給出板表面聲壓的分布。

3.1.2 板表面法向振速的識別

采用有限元方法計算得到的板表面法向振速如圖3所示?；诼晧禾荻葴y量并采用偶極子源得到板表面聲壓分布后，按照公式（12）僅考慮平板在聲場中的位置和尺寸識別的法向振速如圖4所示，其中Δ1也取為0.001 m。按照公式（13）考慮平板表面阻抗關系識別的法向振速如圖5所示。

從圖4和圖5的結果可以看出，采用公式（12）基于聲壓梯度計算平板表面法向振速，僅能夠得到板表面法向振速的大致分布，但是該方法計算簡單，對未知聲源表面阻抗關系的情況非常實用。按照公式（13）考慮了板平面的阻抗關系后，則能夠比較準確地識別平板法向振速幅值分布。

圖3 有限元法計算的板表面法向振速Fig.3 The normal velocity calculated by FEM

圖4 按照公式（12）識別的板表面法向振速Fig.4 The normal velocity of the plate identified by formula（12）

圖5 按照公式（13）識別的板表面法向振速Fig.5 The normal velocity of the plate identified by formula（13）

3.2 水下圓柱殼法向振速識別

采用有限元和邊界元耦合方法對浸入無限水深的有限長圓柱殼的聲輻射進行計算。圓柱殼的幾何尺寸及材料常數為：殼長度L＝0.5 m，半徑R＝0.104 5 m，殼的厚度h＝0.003 m，彈性模量E＝206 GPa，密度ρs＝7 850 kg／m3，泊松比γ＝0.3，水中聲波傳播速度c＝1 500 m／s，密度ρ0＝1 000 kg／m3。圓柱殼的幾何中心取為原點。邊界條件僅以一端固定，另一端自由為例。在圓柱的（0.09，－0.052，0.036）處加一幅值為10 N，沿X軸正方向的簡諧力。計算頻率f＝500 Hz，對應于圓柱殼的（3，1）梁模態(tài)；流體取為水，聲速為1 500 m／s，密度為1 000 kg／m3。圓柱殼表面離散為554個節(jié)點，552個shell單元。

相位共軛陣列信息為：陣列與圓柱的徑向距離R1＝0.15 m＝0.05λ，高度為0.6 m，陣元個數為60。陣元間距為：軸線方向Δ＝0.3 m＝0.1λ，周向Δ＝0.08 m＝0.03λ。算例中陣列處的聲壓梯度同樣取為：＝（p2′－p1′）／Δ1，其中Δ1取為沿圓柱殼徑向0.000 1 m。

采用耦合方法得到圓柱殼的表面聲壓和法向速度幅值分布，如圖6和圖7所示。

圖6 圓柱殼表面聲壓分布（單位：Pa）Fig.6 The surface pressure of the cylinder

圖7 圓柱殼表面法向速度分布（單位：m／s）Fig.7 The surface normal velocity of the cylinder

與平板法向振速識別方法一樣，基于聲壓梯度測量并采用偶極子源得到的圓柱殼表面聲壓分布，如圖8所示。得到識別的圓柱殼表面聲壓分布后，按照公式（12）僅考慮圓柱殼在聲場中的位置和尺寸識別的法向振速，如圖9所示，其中也取為0.000 1 m。按照公式（13）考慮圓柱殼表面阻抗關系識別的法向振速，如圖10所示。

圖8 基于PCD得到的表面聲壓幅值分布Fig.8 The surface pressure reconstructed by PCD

圖9 按照公式（12）識別的圓柱殼表面法向振速Fig.9 The normal velocity of the cylinder shell identified by formula（12）

圖10 按照公式（13）識別的圓柱殼表面法向振速Fig.10 The normal velocity of the cylinder shell identified by formula（13）

從圖9和圖10的結果可以看出，采用公式（12）基于聲壓梯度計算圓柱殼表面法向振速，能夠得到表面法向振速的大致分布，而按照公式（13）考慮了圓柱殼表面的阻抗關系后，則能夠比較準確地識別圓柱殼法向振速幅值分布。

4 識別誤差分析

為說明本文方法的有效性，對結構表面法向振速的歸一化幅值進行誤差分析。圖11是基于聲壓梯度測量采用偶極子源方法考慮板阻抗關系時，沿板長邊方向y＝0.1 m時歸一化的法向振速幅值的識別值，與有限元法計算得到結果的比較曲線。圖12是圓柱殼表面上z＝0處，采用60陣元基于PCD方法得到周向歸一化的法向振速識別值和采用耦合方法計算值的比較曲線。

圖11 沿板長邊方向y＝0.1 m時質點振速比較曲線Fig.11 The normalized particle velocity amplitudes at y＝0.1 m

圖12 徑向z＝0法向振速幅值分布比較曲線Fig.12 The normalized pressure amplitude at z＝0 circumferential ring

通過圖11和12可以看到，基于聲壓梯度測量采用偶極子源并考慮結構表面阻抗關系時，識別的歸一化表面法向振速幅值與有限元或耦合方法計算得到的結果吻合較好，并且平板的識別分辨率為0.15 m／0.7 m＝0.2λ，圓柱殼的識別分辨率為0.08 m／3 m＝0.03λ，均突破了聲波的衍射極限。

基于聲壓梯度測量計算出板表面法向振速后，重新計算了測量陣列處的聲壓級。圖13給出了沿測量陣列處長邊方向y＝0.1 m時，聲壓級識別值與測量值的比較曲線，可以看到重新計算得到的測量陣列處的聲壓級值與測量值基本一致，識別誤差較小。

圖13 測量陣列處長邊方向y＝0.1 m聲壓幅值Fig.13 The pressure amplitudes at y＝0.1 m

5 結論

本文基于相位共軛方法對平板和水下圓柱殼輻射聲場的識別進行了數值仿真計算，得到結構表面聲壓分布后，通過兩種方法識別結構的法向振速，得到如下結論：

1）基于聲壓梯度計算結構法向振速，能夠得到法向振速幅值的大致分布，該方法計算簡單，對于無法得到聲源表面阻抗關系的情況非常實用，而引入結構表面阻抗關系則能得到更加精確的結果。

2）通過誤差分析可知，基于聲壓梯度測量采用偶極子源并考慮結構表面阻抗關系時，識別的歸一化表面法向振速幅值與有限元或耦合方法計算得到的結果吻合較好，突破了聲波的衍射極限。通過識別的板表面法向振速，重新計算得到的測量陣列處的聲壓級值也與測量值基本一致。

［1］MAYNARD J D，WILLIAMS E G.Nearfield acoustic holography （NAH）：I.Theory of generalized holography and the development of NAH［J］.Journal of the Acoustical Society of America，1985，78（4）：1395－1413.

［2］JACKSON D R，DOWLING D R.Phase conjugation in underwater acoustics［J］.Journal of the Acoustical Society of America，1991，89（1）：171－181.

［3］LEROSEY G，ROSNY J D，TOURIN A，et al.Focusing beyond the diffraction limit with far－field time reversal［J］.Science，2007，315：1120-1122.

［4］CONTI S G，ROUX P，KUPERMAN W A.Near－field time－reversal amplification［J］.Journal of the Acoustical Society of America，2007，121（6）：3602-3606.

［5］ROSNY J D，FINK M.Overcoming the diffraction limit in wave physical using a time－reversal mirror and a novel acoustic sink［J］.Physical Review Letters，2002，89（12）：124301.

［6］ROSNY J D，FINK M.Focusing properties of near－field time reversal［J］.Physical Review A，2007，76：065801.

［7］CARMINATI R，SAENZ J J，GREFFET J J，et al.Reciprocity，unitarity and time reversal symmetry of the S matrix of fields containing evanescent components［J］.Physical review A，2000，62（1）：012712（7）.

［8］FINK M，CASSEREAU D，DERODE A，et al.Time－reversed acoustics［J］.Reports on Progress in Physics，2000，63（12）：1993－1995.

［9］EVERSTINE G C，HENDERSON F M.Coupled finite element／boundary element approach for fluid－structure interaction［J］.Journal of the Acoustical Society of America，1990，87（5）：1938－1947.

［10］FAHY F.Sound and structural vibration：radiation，transmission and response［M］.London：Academic Press，1985.

Identification of the Surface Normal Velocity of a Structure Based on Phase Conjugation Method

Liu Song1，2Li Sheng1，2Zhao De－you2
1 State Key Laboratory of Structural Analysis for Industrial Equipment，Dalian University of Technology，Dalian 116024，China 2 School of Naval Architecture，Faculty of Vehicle Engineering and Mechanics，Dalian University of Technology，Dalian 116024，China

s：The identification of the underwater radiated sound field of a plate and cylindrical shell were studied numerically by using the phase conjugation method.The normal velocity distribution of the structure was reconstructed by two different methods after obtaining the surface pressure distribution.One method was based on the pressure gradient and only required the dimension and the location of the structure in the sound field，the other was based on the surface impendence relationship of the structure vibration and the surface pressure distribution.The results show that the pressure gradient method can get the approximate distribution of the surface normal velocity and is simple to implement even if the surface impendence relationship is unknown.The method based on the surface impendence relationship achieves more accurate identification results of the surface normal velocity.

phase conjugation；normal velocity；identification；pressure gradient；impendence

TB532，U661.44

：A

：1673－3185（2011）04－19－06

2010-08-06

國家自然科學基金 （10972046）

劉 松（1982－），男，博士研究生。研究方向：噪聲源識別和定位及聲場重構。E-mail：Lius1982＠gmail.com

黎 勝（1973－），男，副教授，博士生導師。研究方向：船舶與海洋結構物振動噪聲分析及控制。E-mail：Shengli＠dlut.edu.cn

10.3969/j.issn.1673-3185.2011.04.004