有限浸深圆柱壳振动及远场声辐射的解析方法

郭文杰 ，李天匀1，2，3，朱翔1，张帅1

1 430074华中科技大学 船舶与海洋工程学院，湖北 武汉2 200240高新船舶与深海开发装备协同创新中心，上海3 430074船舶与海洋水动力湖北省重点实验室，湖北 武汉

摘 要：［目的］针对目前对于自由液面影响下圆柱壳—流场耦合系统振动及声辐射解析研究的匮乏，提出一Graf种有限浸没深度下有限长圆柱壳振动及远场声辐射的解析求解方法。［方法］采用镜像原理和 加法定理得到流体速度势的解析表达式，然后再结合能量泛函变分方法推导出计及自由液面影响的壳—液耦合振动方程，从而可以求解系统受迫振动响应。［结果］研究表明，相比于无限域，自由液面的存在会增大同阶次共振频率，但Nastran随着浸没深度的逐渐增加，均方振速很快趋于无限域工况。与 软件计算结果对比表明所提出的方法准确、可靠，且具有方法简便、计算量小的优点。利用求得的振动响应，通过傅里叶变换和稳相法可得到远场辐射声压，计算结果表明，自由液面会使得远场声压指向性和波动性出现类偶极子效应；但是不同于振动特性，远场声压并不会随浸没深度增大而很快趋于无限域工况。［结论］所提出的方法实现了外力激励下计及自由液面影响的水下圆柱壳远场声辐射快速预报，对于半空间结构声振问题的研究具有一定的指导意义。关键词：自由液面；镜像原理；Graf加法定理；远场辐射声压中图分类号：U661.44 文献标志码：A DOI：10.3969/j.issn.1673-3185.2017.04.010

Analytical research of vibration and far-field acoustic radiation of cylindrical shell immersed at finite depth

GUO Wenjie1，LI Tianyun1，2，3，ZHU Xiang1，ZHANG Shuai1 1 School of Naval Architecture and Ocean Engineering，Huazhong University of Science and Technology， Wuhan 430074，China 2 Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration，Shanghai 200240，China 3 Hubei Key Laboratory of Naval Architecture and Ocean Engineering Hydrodynamics，Wuhan 430074，China Abstract：Aiming at the current lack of analytical research concerning the cylindrical shell-flow field coupling vibration and sound radiation system under the influence of a free surface, this paper proposes an analytical method which solves the vibration response and far-field acoustic radiation of a finite cylindrical shell immersed at a finite depth. Based on the image method and Graf addition theorem, the analytical expression of the fluid velocity potential can be obtained, then combined with the energy functional of the variation method to deduce the shell-liquid coupling vibration equation, which can in turn solve the forced vibration response. The research shows that, compared with an infinite fluid, a free surface can increase at the same order of resonance frequency; but as the depth of immersion gradually increases, the mean square vibration velocity tends to become the same as that in an infinite fluid. Compared with numerical results from Nastran software, this shows that the present method is accurate and reliable, and has such advantages as a simple method and a small amount of calculation. The far-field radiated pressure can be obtained by the vibration response using the Fourier transformation and stationary phase method. The results indicate that the directivity and volatility of the far-field acoustic pressure of a cylindrical shell is similar to that of an acoustical dipole due to the free surface. However, the far-field acoustic pressure is very different from the vibration characteristics, and will not tend to an infinite fluid as the submerging depth increases. Compared with the numerical method, the method in this paper is simpler and has a higher computational efficiency. It enables the far-field acoustic radiation of an underwater cylindrical shell to be predicted quickly under the influence of external incentives and the free surface, providing guiding significance for acoustic research into the half space structure vibration problem. Key words：free surface；image method；Graf addition theorem；far-field radiated pressure

0引言

圆柱壳—流场耦合振动及声辐射特性的研究工作很多［1-4］，但是这些工作主要针对的是无限流域这类工况。在实际工程问题中，自由液面是常见的边界类型，因此需要对考虑自由液面影响的壳—液耦合振动及声辐射问题进行研究。自由液面对壳—液耦合振动影响的研究工作较少，Ergin等［5］基于实验和三维水弹性软件对有限浸没深度下圆柱壳振动特性进行分析，发现结构离自由液面越近，同阶次固有频率越大。Amabili 6 ［ ］对部分充液圆柱壳进行研究，提出了用扇形边界替代自由液面的近似方法。之后， Amabili［7 ］将该方法拓展到处理部分浸没问题。随后，Ergin等［8］利用边界积分法和镜像原理对部分充液（浸没）圆柱壳振动特性进行了研究，结果与实验数据符合良好。王鹏等［9-10］基于波传播方法和虚源法，分析了浅水域圆柱壳的自由振动特性，但是该方法假设振型不随浸没深度变化并且忽略了虚源声波对结构振速的影响。Guo等［11］建立了有限浸没深度下圆柱壳流固耦合物理模型，考虑了虚源速度势对结构振速的影响，并采用边界元方法研究了远场声辐射特性。王斌等［12］从圆柱壳表面的均方振速和辐射声功率的角度，对半浸状态和全浸状态下圆柱壳在无限长线激励作用下的声振特性进行比较分析，指出了二者之间的ANSYS差别与联系。刘佩等［13］采用有限元软件 对有限深度浸没圆柱壳进行仿真，得到了与文献［5］类似的结论，并指出自由液面对圆柱壳自由振动4的影响在浸没深度大于 倍半径时可以忽略不计。半空间声学问题的解析研究工作较多， Huang［14］研究了平面波入射下二维圆柱壳的散射声场，Hasheminejad等［15-16］开展了有限空间谐振动二维圆柱的声场研究工作。白振国等［17］采用镜像法建立了有限水深环境中二维圆柱壳的振动声辐射数学物理模型，初步计算了浅水对圆柱壳振动声辐射的影响规律及水深、潜深对声场分布和衰减特性的影响规律。Li等［18］基于镜像原理进行了自由液面下有限潜深无限长圆柱壳结构的声辐射性能研究，基于稳相法，最终得到了有限浸没深度下圆柱壳结构远场辐射声压的计算表达式。但是，这些工作的研究对象皆是无限长圆柱壳，当圆柱壳长度为有限时，半空间声学问题的求解难度较大，目前尚缺乏解析研究。圆柱壳振动时会向外辐射声波，由于自由液面的反射作用，反射声波会在结构表面发生刚性 散射，而散射声触及自由液面又会形成回波，继而在结构表面产生多次散射和在界面产生多次反射（互散射），并最终形成稳态声场。由于辐射声和散射声的相似性，可将其统一表示。基于镜像原理，反射声均可认为由虚源发出，故将所有声波分2为 类，即实源声和虚源声。若流体假设为不可压缩，虽然流体中不存在声波，但是也可通过镜像原理，将速度势设为实源速度势和虚源速度势来分析此类问题。本文将基于此类思路分析计及自由液面影响的有限长圆柱壳振动响应及声辐射问题，实际上考虑了互散射效应。采用镜像原理来处理自由液Graf面处声压为零的边界条件，再利用 加法定理2对实源和虚源这 种坐标系进行转换得到速度势在流场的分布，然后再结合能量泛函变分方法得到壳—液耦合振动方程，进而可以求解其振动特性。之后，利用傅里叶变换构建波数域壳体表面连续条件，再利用傅里叶逆变换和稳相法求得其远场声压。

1 理论分析

圆柱壳长度为L，厚度为h，中面半径为R，浸没深度为H，u v w分别表示轴向、周向和径向的中面位移，壳体材料的密度为 ρ ，弹性模量为E，泊松比为 μ ，流体密度为 ρ 。取圆柱壳左端面中f心为坐标原点o，对应直角坐标 (x y z) 。实际分析中选择柱坐标系 (x r φ) ，其中 x 表示轴向，r表1示径向，φ为周向角（与 y 轴的夹角），如图 所示。

1.1 振动分析

Wmn 为三向位移幅值；km =m π/L；ω 为角频率；t为时间；i为虚数单位。本文采用能量泛函变分的方法研究有限浸没深度圆柱壳的振动特性，故首先应得到各部分能量的表达式。壳体应变能（基于Love壳体理论［19］）可表示为

U = 1  εT σdV （2） 2 v式中： ε 为应变向量； σ 为应力向量； V为圆柱壳体积分域。根据位移函数的正交性，积分后，式（2）可写成+¥ +¥ U = 1 å å ξ K ξ T （3 ） 2 m n mn m -n m = 1n = -¥式中： ξ = [Umn Vmn Wmn ] ；刚度矩阵 K 为三m n mn Hermite阶 矩阵。壳体动能可表示为ρhω2 L π T = 0  (u2 + v2 + w 2)Rdφdx （4） 2 -π同理，根据位移函数正交性，积分后可表示为+¥ +¥

T = 1 å å ξ M ξ T （5 ） 2 m n mn m -n

m = 1n = -¥

式中，质量矩阵 M 为三阶对角矩阵。mn

为求解流体做功，首先需要得到速度势函数的解析表达式。本文基于势流理论，将流体视为不可压缩、无旋、无粘性的理想流体，因此速度势Laplace函数 ϕ ( r x φ t )满足柱坐标系下的 方程： 1 ¶ ¶ϕ ¶2 ϕ ¶ 2ϕ 6 (r ) + 1 + = 0 （） r ¶r ¶r r 2 ¶φ2 ¶x2

对于水下圆柱壳，满足无穷远处速度势为零的条件： ¶ϕ

（7）

=0 ¶r r =¥

由于自由液面的存在，可以借鉴镜像原理进行分析，认为速度势可由结构振动直接引起的实源速度势和自由液面反射的虚源速度势叠加组( , , )与实源坐标系关于自成。虚源坐标系 x′ r′ φ′ 2由液面对称，如图 所示。o′ z′ 设流域中任意一点为点P，其速度势函数可以表示为ϕ(r x φ t) = ϕr(r x φ t) + ϕi(r′ x′ φ′ t) （8）式中： ϕr ( r x φ t ) 为实源流体速度势；

ϕi(r′ x′ φ′ t)为虚源流体速度势。满足式（6）和式（7）的速度势函数有以下形式： ϕr ( r x φ t )= å+¥ å+¥ r ϕmn Kn (kmr)sin(km x)exp(inφ)sin(ωt) ( )= （9） m = 1n = -¥ ϕi r′ x′ φ′ t

å+¥ å+¥ i ϕmn Kn (kmr′)sin(km x′)exp(inφ′)sin(ωt) m = 1n = -¥式中，Kn（）为第2类修正贝赛尔函数。由于自由液面处速度势为零，故有ϕr(r x φ t) + ϕi(r′ x′ φ′ t) =0 （10）当 P点位于自由液面上时，满足如下位置关系：

r = r′ x = x′ φ = φ′ =π （11）将式（9）和式（11）代入式（10），正交化处理后得到

ϕr mn + (-1)nϕi m -n =0 （ 12 ）即可以得到速度势函数的解析表达式：

ϕ(r x φ t) = å+¥ å+¥ ϕmn[Knr (kmr)exp(inφ)

m = 1n = -¥ (-1)n K (kmr′)exp(-inφ′)]sin(km x)sin(ωt) （13） -n根据Graf加法定理［20］： K -n (kmr ′)exp(inφ′) =

å+¥ (-1)a Ka + n (2kmH )Ia (kmr)exp(iaφ) r＜2H a= -¥ å+¥ + n  (-1)a Ia + n (2kmH )Ka (kmr)exp(iaφ) r 2H a= -¥

（14） 1式中，（）为第Ia a阶第 类修正贝赛尔函数。r≈R <对于有限深度浸没，结构表面处半径2H，因此速度势解析表达式ϕ(r x φ t) = å+¥ å+¥ ϕmn[Kn (kmr)exp(inφ)

m = 1n =-¥ å+¥ a +n

(-1) Ka + n (2km H )Ia (kmr)exp(-iaφ)] a= -¥

（15） sin(km x)sin(ωt)因为系数a和n地位等价，交换顺序级数求和后上式改写为å+¥ å+¥ [ϕmn å+¥ + nϕ

ϕ(r x φ t) = Kn (kmr) - ( -1)a mn m= 1n = -¥ a= -¥ Ka (2km H )In (kmr)]exp(-inφ)sin(km x)sin(ωt) （16） +n

根据圆柱壳外壁面处速度连续条件，有¶ϕ - = - ¶w （17） ¶r ¶t r =R r =R将式（16）代 式（17）中入 ，正交化处理后可以得到速度势幅值向量 φ 与位移幅值向量 ζ 的n n

关系： （18） φ = Qζn n式中：Q 为速度势幅值向量和位移幅值向量的关= =系矩阵； φn {ϕm -N ϕm  ϕm }T ； ζn -N+ 1 N {Wm -N Wm W } T ，即可将速度势幅值-N +1 m N向量用位移幅值向量表示。由伯努利方程可以得到壁面处的流体动压力

¶ϕ

19 P =ρ （ ） f f ¶r r =R

流体做功为

L -ππ 20 W =- 1  P wRdφdx （ ） 0 f f 2点激励力做功为（21） W = F0 w(x φ ) e 0 0式中， F 为激励力幅值，激励力位置为 (x0 φ 0) 0处，激励力频率为f。由上述各能量分量可得到能量泛函表达式： （22） Õ= U- W - T - We f根据变分原理，满足¶Õ ，¶Õ ，¶Õ （23） = 0 = 0 = 0 ¶Umn ¶Vmn ¶Wmn由对幅值Umn ，Vmn 的偏导为0可以得到其与幅值Wmn的线性关系，简写为如下所示的形式： ¶Õ = a1Um + b1Vm + c1Wm =0 ¶Umn -n -n -n 24 （ ） ¶Õ = a 2Um + b 2Vm + c2Wm =0 ¶Vmn -n -n -n式中， a1 ，b1 ，c1 和 a ，b ，c2 都是关于角频率2 2 ω 、刚度矩阵 K 和质量矩阵 M 的系数，即mn mn

Umn ，Vmn 可由Wmn 进行代换。由于轴向函数在域内正交，当n的截断数为-N～N时，最终根据 ¶Õ ¶Wmn = 0可以得到轴向波数m取任意值时的方程： 25 T ζ = γn （ ） m n =式中： γ 表示激励力展开后的幅值向量， γ F0 n n 0){exp(iNθ 0)}T sin(km x 0) exp[i(N - 1)θ 0]  exp(-iNθ ； 2N+1 T 为 阶矩阵。m 求解自由振动时，γ 为零向量，因为 ζ 中元n n

素不全为 0 ，所以 T 必然不是满秩矩阵，即m det(T m) = 0。由此可以求解出轴向波数m取任意值时各阶角频率ω，从而可以得到固有频率值。求解受迫振动时，给定激励力频率f，可求出对应的 ζ = T m-1γ ，从而得到结构任意位置的振n n动响应。1.2 远场声辐射分析本文采用傅里叶变换方法，将轴向x变换到波数域，构建新的壁面连续条件（假设圆柱两端有两个半无限长声障柱），再进行逆变换即可求得声压表达式。定义傅里叶变换及逆变换形式如下： +π

f(k) = 1  f (x)exp(-ikx)dx 2π -¥ 26 （ ） +π

f(x) = -¥ f (k)exp(ikx)dk

式中，为波数，由此可将轴向坐标转化到波数域。k声压可划分为实源声压 P 和虚源声压 P ， r i根据分离变量法，实源声压可以表示为+¥

å 27

P (r x φ)= p (r x)exp(inφ) （ ） r n n= -¥

式中，pn为与周向角无关的声压物理量。Helmholtz波数域流体声压满足 方程： ¶2 p͂ (r k) ¶p͂ (r k) n + 1 n + (k2 - k - n2 2 ) p͂ (r k) =0 ¶r2 r ¶r r n f 2 （28）式中， k = ω c ，为压缩波数，其中 c 为流体声f f f速，角频率 ω = 2πf ，f 为谐振动频率。h h式（28由 ）可以得到实源傅氏声压如下形式的解： +¥

å (1) （29）

P r(r k φ) = An (k)Hn (krr)exp(inφ) n= -¥

(1) 1 Hankel式中： H () 为第 n 阶第 类 函数； n kr = kf 2 - k 2 ，为径向波数；A (k)为对应的波数域n幅值。同理，虚源傅氏声压可以表示为+¥

å (1) 30

P i(r′ k φ′) = Bn (k)Hn (krr ′ )exp(inφ ′) （ ） n= -¥

式中，Bn (k)为虚源声波数域幅值。考虑到空气中的波阻抗远小于水中的波阻抗，声波由水中射向空气，自由液面处可以看成是绝对软的边界，自由液面处的声压可以近似认为是零，则自由液面处某点（ r = r′ ，φ = π - φ′ ）的傅氏声压满足 P +P = 0 ，可以推导出r i （31） Bn (k) =- A (k) -n

根据柱贝塞尔函数的Graf加法定理可以实现坐标迁移：

H (1) (krr ′)exp(-inφ′) må= ¥ -¥ (-1) m + mH m(1) + -nn (2kr H )Jm (kr r)exp(imφ)= r＜2H m å+¥ =¥ (-1)m + nJ m+ n (2kr H )Hm (1) (krr)exp(imφ) r  2H

（32）式中，Jn () 为第n阶第1类贝塞尔函数。同理，将结构径向位移 w 变换到波数域，有w͂ (k φ) = n å= +¥ -¥ w n (k)exp(inφ) （ 33 ）

Wmnkm [1 - (-1)m exp(-ikL)] 。式中，w n (k )= 2π1 m åM = km2 -k 2 1

根据壁面处连续条件，变换到波数域，有

-¶P +¶ Pi （34） r = ρ ω2 w͂ (k φ) ¶r f r= R

将虚源傅氏声压迁移到实源坐标系下，因r≈R <2 H，有P i = - å+¥ å+¥ (-1)m +n A (k)Hm (1) (2kr H )Jm (krr)exp(imφ) n + n n= -¥ m = -¥

35 （ ）交换积分顺序后，有

P = i +¥ +¥ - å å (-1)m +n A (k)Hm (1) (2kr H )Jn (kr r)exp(inφ) m + n n= -¥ m = -¥

（36）将式（36）代入式（34）后，级数做有限截断，均-N从 取到 N ，可以得到波数域声压幅值向量An (k)与位移幅值向量w (k) 的关系： n A (k)= T w (k) （37） n ran n式中： A (k) ={A (k) A 1(k)  AN (k)}T ；T n -N -N+ ran

为波数域声压幅值到位移幅值的迁移矩阵； w (k) ={w͂ (k) w͂ 1(k)  w͂ (k)}T 。n -N -N+ N通过傅里叶逆变换，即可求出任意场点的声压。

+¥ （38） P(r x φ) = -¥ P͂(r k φ)exp(ikx)dk对于远场声压，可以采用稳相法［21］求解。限于篇幅，本文略去详细推导过程，直接给出球坐标

系下声压表达式： P ( R 0 θ φ)= -2i exp(ik f R 0) ´ r R0

åN An (k cos θ)exp[in(φ - π 2)] n= -N f

(0 2i exp(ikf R 0) P R  θφ )= ´ i R0

N N å å

Aa (k cos θ)(-1)n +a Jn (2k H sin θ) f + a f n =- N a= -N exp[in(φ - π 2)]

（39）

式中：A (k cos θ) 可由式（ 37 ）求得；R0 为场点到a f原点的距离；θ 为观测角（与轴向的夹角）。此处是将柱坐标系（ x r φ ）转换到球坐标系 (R  θ φ) 0下，有 r = R0 sin θ ，x = R0 cos θ 。2 数值分析2.1 振动分析=1.284 m，半径R=计算模型参数如下：壳长L 0.18 m，厚度 h=0.003 m，壳 =7 850 kg/m3，体密度 ρ E=206 GPa =0.3杨氏模量 ，泊松比 μ ，流体密度=1 025 kg/m3，流 =1 500 m/s。点谐激ρ 体声速 c f f励力位置坐标（R，L/2，0），激励力幅值 =1。F0 2.1.1 收敛性分析

H=0.2 2m，为了说明方法收敛性，本文取 和200 400 Hz分别计算在激励力频率为 和 时径向均方振速Vm 随截断数M，N的变化规律，为便于分1 s析收敛性，假设M=N。定义 V = |Vn |2ds ，其m 2s

中，s为表面积，Vn为径向速度。3从图 可以看出，随着截断项数的增加，径向M=N=10均方振速Vm 值很快趋于稳定，当 时已经10。足够收敛，因此本文算例截断项数取值均为

Fig.2图2 镜像原理示意图Schematic diagram of image method

Fig.1图1 模型及坐标系Model and coordinate system

Fig.3图3 径向均方振速收敛性分析The convergence curves of the quadratic velocity