基于相位共軛方法識(shí)別結(jié)構(gòu)表面法向振速

劉 松 黎 勝 趙德有

1大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連116024 2大連理工大學(xué) 運(yùn)載工程與力學(xué)學(xué)部 船舶工程學(xué)院，遼寧 大連116024

基于相位共軛方法識(shí)別結(jié)構(gòu)表面法向振速

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

1大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連116024 2大連理工大學(xué) 運(yùn)載工程與力學(xué)學(xué)部 船舶工程學(xué)院，遼寧 大連116024

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

相位共軛；法向振速；輻射聲場(chǎng)識(shí)別；聲壓梯度；阻抗

1 引言

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

2 理論

2.1 結(jié)構(gòu)聲輻射理論

結(jié)構(gòu)在簡(jiǎn)諧力作用下考慮流體加載效應(yīng)的有限元形式的運(yùn)動(dòng)方程為［9］：

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

｛fp｝和結(jié)構(gòu)表面聲壓｛p｝的關(guān)系式為：

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

對(duì)從聲源向外輻射正向傳播的聲場(chǎng)，根據(jù)Helmholtz－Kirchhoff積分公式，聲源表面S的聲壓p（Q）和聲壓梯度（或法向振速vn（Q））的關(guān)系為：

式中，［Z］為阻抗矩陣；｛vn｝為結(jié)構(gòu)表面法向速度向量。

由｛vn｝和結(jié)構(gòu)速度向量｛v｝之間的轉(zhuǎn)換關(guān)系式：

可得求解結(jié)構(gòu)－聲耦合問(wèn)題的方程為：

求出結(jié)構(gòu)速度向量｛v｝后，進(jìn)而可求出｛vn｝和｛p｝等。

對(duì)簡(jiǎn)諧激勵(lì)力作用下位于無(wú)限障板上的板結(jié)構(gòu)在板一側(cè)半無(wú)限域流體介質(zhì)和板表面Sp上產(chǎn)生的輻射聲壓p（P）也可由Rayleigh積分求得［10］：

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

2.2 相位共軛理論

對(duì)反向傳播的聲場(chǎng)，根據(jù)Helmholtz－Kirchhoff積分公式，在相位共軛陣列測(cè)量到聲壓和聲壓梯度后，其相位共軛聲場(chǎng)為：

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

實(shí)際的相位共軛陣列都是離散的，對(duì)包含N個(gè)陣元的離散的有限陣列，其相位共軛聲場(chǎng)為：

如基于測(cè)量聲壓使用單極子源來(lái)進(jìn)行時(shí)間反轉(zhuǎn)，其相位共軛聲場(chǎng)為：

如基于測(cè)量聲壓梯度使用偶極子源來(lái)進(jìn)行時(shí)間反轉(zhuǎn)，其相位共軛聲場(chǎng)為：

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

2.3 結(jié)構(gòu)法向振速識(shí)別

得到識(shí)別的結(jié)構(gòu)表面聲壓分布后，利用以下兩種方法識(shí)別結(jié)構(gòu)的法向振速。

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

式中，p1為結(jié)構(gòu)在聲場(chǎng)中所在位置識(shí)別的聲壓值，p2是與p1對(duì)應(yīng)的沿結(jié)構(gòu)表面法線方向增加距離為Δ1時(shí)識(shí)別的聲壓值。

方法二：除了考慮結(jié)構(gòu)在聲場(chǎng)中的位置和尺寸之外，還可以通過(guò)引入結(jié)構(gòu)表面阻抗關(guān)系對(duì)法向振速進(jìn)行識(shí)別。將識(shí)別的結(jié)構(gòu)表面聲壓代入公式（4），并對(duì)阻抗矩陣求逆，得到法向振速，如式（13）所示。

3 數(shù)值模擬

3.1 平板輻射聲場(chǎng)和法向振速識(shí)別

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

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

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

3.1.1 板表面聲壓的識(shí)別

采用Rayleigh積分計(jì)算和按照公式（11）基于聲壓梯度測(cè)量使用偶極子源（PCD）方法識(shí)別的平板表面聲壓分布，如圖1和圖2所示（聲壓?jiǎn)挝唬篜a；法向速度單位：m／s）。

圖1 Rayleigh積分計(jì)算的結(jié)構(gòu)表面聲壓Fig.1 The surface pressure calculated by Rayleigh integral

圖2 PCD識(shí)別的表面聲壓結(jié)果Fig.2 The surface pressure reconstructed by PCD

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

3.1.2 板表面法向振速的識(shí)別

采用有限元方法計(jì)算得到的板表面法向振速如圖3所示。基于聲壓梯度測(cè)量并采用偶極子源得到板表面聲壓分布后，按照公式（12）僅考慮平板在聲場(chǎng)中的位置和尺寸識(shí)別的法向振速如圖4所示，其中Δ1也取為0.001 m。按照公式（13）考慮平板表面阻抗關(guān)系識(shí)別的法向振速如圖5所示。

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

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

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

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

3.2 水下圓柱殼法向振速識(shí)別

采用有限元和邊界元耦合方法對(duì)浸入無(wú)限水深的有限長(zhǎng)圓柱殼的聲輻射進(jìn)行計(jì)算。圓柱殼的幾何尺寸及材料常數(shù)為：殼長(zhǎng)度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。圓柱殼的幾何中心取為原點(diǎn)。邊界條件僅以一端固定，另一端自由為例。在圓柱的（0.09，－0.052，0.036）處加一幅值為10 N，沿X軸正方向的簡(jiǎn)諧力。計(jì)算頻率f＝500 Hz，對(duì)應(yīng)于圓柱殼的（3，1）梁模態(tài)；流體取為水，聲速為1 500 m／s，密度為1 000 kg／m3。圓柱殼表面離散為554個(gè)節(jié)點(diǎn)，552個(gè)shell單元。

相位共軛陣列信息為：陣列與圓柱的徑向距離R1＝0.15 m＝0.05λ，高度為0.6 m，陣元個(gè)數(shù)為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

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

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

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

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

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

4 識(shí)別誤差分析

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

圖11 沿板長(zhǎng)邊方向y＝0.1 m時(shí)質(zhì)點(diǎn)振速比較曲線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

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

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

圖13 測(cè)量陣列處長(zhǎng)邊方向y＝0.1 m聲壓幅值Fig.13 The pressure amplitudes at y＝0.1 m

5 結(jié)論

本文基于相位共軛方法對(duì)平板和水下圓柱殼輻射聲場(chǎng)的識(shí)別進(jìn)行了數(shù)值仿真計(jì)算，得到結(jié)構(gòu)表面聲壓分布后，通過(guò)兩種方法識(shí)別結(jié)構(gòu)的法向振速，得到如下結(jié)論：

1）基于聲壓梯度計(jì)算結(jié)構(gòu)法向振速，能夠得到法向振速幅值的大致分布，該方法計(jì)算簡(jiǎn)單，對(duì)于無(wú)法得到聲源表面阻抗關(guān)系的情況非常實(shí)用，而引入結(jié)構(gòu)表面阻抗關(guān)系則能得到更加精確的結(jié)果。

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

［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，F(xiàn)INK 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，F(xiàn)INK 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，F(xiàn)aculty 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

國(guó)家自然科學(xué)基金 （10972046）

劉 松（1982－），男，博士研究生。研究方向：噪聲源識(shí)別和定位及聲場(chǎng)重構(gòu)。E-mail：Lius1982＠gmail.com

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

10.3969/j.issn.1673-3185.2011.04.004