基于相位共轭方法识别结构表面法向振速

刘 松 黎 胜 赵德有

1大连理工大学 工业装备结构分析国家重点实验室，辽宁 大连116024 2大连理工大学 运载工程与力学学部 船舶工程学院，辽宁 大连116024

基于相位共轭方法识别结构表面法向振速

刘 松1，2黎 胜1，2赵德有2

1大连理工大学 工业装备结构分析国家重点实验室，辽宁 大连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）可以看到，以上基于相位共轭方法通过计算反向传播的相位共轭声场来进行声源定位和识别具有实施简单，可直接使用平面阵列等来进行局部重建，不存在解不稳定性问题等优点。本文数值模拟计算均基于声压梯度测量使用偶极子源（PCD）识别结构表面声压。

2.3 结构法向振速识别

得到识别的结构表面声压分布后，利用以下两种方法识别结构的法向振速。

方法一：已知结构在声场中所处的位置及尺寸，以及离散的节点和单元数，基于声压梯度按照式（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）模态；板边界条件取为四边简支，并嵌入无限大刚性障板中；流体取为空气，声速为350 m／s，密度为1.21 kg／m3。

相位共轭测量阵列信息为：阵元距离板表面d＝0.1λ＝0.07 m，阵元个数221，阵元间距Δ＝0.05 m＝0.07λ，与板形状一致。

本节数值模拟计算将按照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所示。基于声压梯度测量并采用偶极子源得到板表面声压分布后，按照公式（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）梁模态；流体取为水，声速为1 500 m／s，密度为1 000 kg／m3。圆柱壳表面离散为554个节点，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