中国科技核心期刊      中国指挥与控制学会会刊     军事装备类重点期刊
Weapon & Information System

Rigid Spherical Acoustic Scattering Simulation Based on LMS VirtualLab

  • ZHU Han
Expand
  • 713th Research Institute of China State Shipbuilding Corporation Limited, Zhengzhou 450000, China

Received date: 2022-06-21

  Revised date: 2022-07-19

  Online published: 2022-10-20

Abstract

Aiming at the calculation problem of underwater target sound scattering, this paper simulates the scattered sound field of a rigid sphere under the plane sound wave. Boundary element method is used to calculate the directivity function of rigid sphere under different ka values, and compare with the theoretical analytical solution. The results show that the spherical scattering sound field and its directivity corresponding to different ka values have obvious differences. As the ka value increases, the directivity becomes more obvious, and the forward scattering intensity increases. The solution calculation results are in good agreement. When the output sound pressure value at the calculation point is small, the default sound pressure output accuracy in the parameter settings of the Virtual.Lab Acoustics module is improved. Export calculation results that meet the needs. Boundary element method can be used to predict the acoustic scattering characteristics of actual underwater target.

Cite this article

ZHU Han . Rigid Spherical Acoustic Scattering Simulation Based on LMS VirtualLab[J]. Command Control and Simulation, 2022 , 44(5) : 102 -106 . DOI: 10.3969/j.issn.1673-3819.2022.05.018

仿真技术在现代工程中得到广泛应用,在声学仿真领域,球体的散射声场仿真具有严格的理论解析解,经常作为标准模型用于验证各种方法和理论[1-3],同时,由于球体本身可以作为某些特定问题的简化物理模型,常被作为简化模型进行实际问题的仿真计算[4-5]。因此,球体的声散射仿真研究一直是人们关注的问题。本文首先对刚性球体散射声场远场理论解析解的推导过程进行了说明,并计算了刚球远场散射声压,其次采用直接边界元方法,应用Virtual.Lab软件对水下刚性球体散射声场进行仿真计算,并将计算结果与理论解析解进行比较,得到了令人满意的结果。

1 理论解析解

1.1 波动方程及其解

图1所示,刚性不动球位于无限流体介质中,半径为a,表面光滑,在沿x轴方向传播的平面波P0ej(kx-ωt)作用下,其散射声场声压ps在球坐标系中满足波动方程:
1 r 2 r r 2 p s r+ 1 r 2 s i n θ θ s i n θ p s θ+ 1 r 2 s i n 2 θ α 2 p s φ 2+k2ps=0
式中,rθφ为球坐标系中的坐标变量,k=ω/c为波数,等于入射声波圆频率ω和球体周围介质中的声速c之比。考虑球体的对称性,且入射波和散射波关于x轴对称,因此,方程应与变量φ无关,式(1)简化为
1 r 2 r r 2 p s r+ 1 r 2 s i n θ θ s i n θ p s θ+k2ps=0
采用分离变量法求解上述方程可得
ps= m = 0 amPm(cosθ) h m ( 1 ) (kr)
式中,am为待定常数,由边界条件确定。
图1 平面波在球面上的散射

1.2 边界条件

刚性球体表面介质质点径向振速为零,即
ur|r=a= i ρ 0 ω p r|r=a=0
式中,ρ0是球体周围介质中的密度,p是介质中的总声压,即入射声压和散射声压之和,ur是介质质点振速的径向分量,即入射波引起的介质质点振速的径向分量和散射波引起的介质质点振速的径向分量之和。
用勒让德函数(Legendre Function)和球贝塞尔函数(Spherical Bessel Function)表示入射波,并利用关系式:
eikrcosθ= m = 0 (2m+1)imjm(kr)Pm(cosθ)
可求得常数am
am=[-im(2m+1)P0 rjm(kr)/ r h m ( 1 ) (kr)]|r=a
式中,jm(kr)为m阶球贝塞尔函数。

1.3 散射声场的远场解

将式(6)代入式(3)并加入时间因子e-iωt得到散射声表达式为
ps= m = 0 -im(2m+1)P0 j ' m ( k a ) h ' m ( 1 ) ( k a )·Pm(cosθ) h m ( 1 ) (kr)·e-iωt
式(7)是散射声场的一般表达式,但人们一般更关心它的远场特性,因此,可用球汉克尔函数(Spherical Hankel Function)在大宗量条件下的渐近展开式:
h m ( 1 ) ( k r ) k r 1 k r e i ( k r - m + 1 2 π )
代入式(7)得
ps=- P 0 k rei(kr-ωt) m = 0 im(2m+1) j ' m ( k a ) h ' m ( 1 ) ( k a )· e - i m + 1 2 π·Pm(cosθ) kr≫1
若记
D(θ)= 1 k a m = 0 im(2m+1) j ' m ( k a ) h ' m ( 1 ) ( k a ) e - i m + 1 2 πPm(cosθ)
则式(9)转化为
ps(r,θ)=- P 0 a rD(θ)ei(kr-ωt)
式(11)为刚性不动球散射声场的远场解,其中,rθ为球坐标的坐标变量,P0为入射声波声压幅值,a为刚球半径,D(θ)为指向性函数,jmm阶第一类球贝塞尔函数, h m ( 1 )m阶第一类球汉克尔函数,Pm为勒让德函数,其他符号定义同前所述。

1.4 解析解计算

上文中提到的指向性函数D(θ)是用于描述刚球散射声场的一个重要参量,因此,为计算刚球的指向性函数,计算模型选取半径为1 m的刚性不动球,入射声波选取单位波幅平面波,取计算阶数m=100,无因次系数ka={1,2,3,4,5,6,7,8},声速c=1 500 m/s,为满足远场条件,需要计算获得距刚球L2之外的声压,其中,L为刚球尺度,λ为声波波长,因此,计算r=50 m处的散射声压,并利用球面衰减规律归算到刚球等效声中心1 m处,如图2所示。
图2 不同ka值下的刚球指向性图

2 边界元方法

2.1 Helmholtz积分方程

刚球声散射计算属于外部问题[6],即所要计算的声场变量都在刚球封闭表面外部,满足Helmholtz方程:
2 p + k 2 p = 0 2 G + k 2 G = 0
式中,p为声压,G为三维自由空间格林函数:
G(X,Y)≡ e - i k r 4 π r
式中,rr(X,Y)是体积V中任意两点XY间的距离,i= - 1,k为波数。
图3所示,刚球表面为S,在体积为V表面为Σ半径无限大的区域内,除X点为格林函数G的奇异点外,其余各点的声压在这个空间内是平稳和非奇异的,因此,可环绕X点做半径ε→0的表面为σ的小球,在表面积为S+σ+Σ的空间V-Vσ内就不再有奇异性。
图3 计算区域定义
将式(12)代入格林公式,则有
S+σ p ( Y ) G ( X , Y ) n - G ( X , Y ) p ( Y ) ndS(Y)=0
经推导可知,式(13)左边积分在σ上为-p(X),在Σ上为零,因此,在S上的积分为p(X)。同理可得,点X在刚球表面S上时,积分值为 1 2p(X),即
S p ( Y ) G ( X , Y ) n - G ( X , Y ) p ( Y ) ndS(Y)= o ( X )     X V 1 2 p ( X ) X S
式中,V表示球体外部空间,S表示球体表面。

2.2 边界元方法

将式(15)中的球体表面S离散成四边形单元,则第i个单元的全局参量pi和其沿法线方向的导数 p i n与单元的节点上参量pij和其导数 p i j n的关系为
p i ( ξ , η ) = j = 1 4 N j ( ξ , η ) p i j p i ( ξ , η ) n = j = 1 4 N j ( ξ , η ) p i j n
式中,Nj(ξ,η)为单元形函数,当单元为四节点四边形时,有
Nj(ξ,η)= 1 4(1+ξjξ)(1+ηjη)(i=1,2,3,4)
式中,ξjηj为单元第j个顶点的局部坐标,其值为
jj)= ( - 1 , - 1 )     j = 1 ( 1 , - 1 ) j = 2 ( 1,1 ) j = 3 ( - 1,1 ) j = 4
将式(16)、(17)代入离散后的式(15),经整理得
[A]{p}=[B] p n
式中,[A]和[B]为系数矩阵,若p表示球面声压,则式(19)给出了声压及其法向导数的关系,通过求解该微分方程组即可得到刚球表面声压分布,进而得到空间V内任意一点的散射声压值。

2.3 Virtual.Lab Acoustics声学仿真计算

平面波作用下,刚球散射声场可以采用Virtual.Lab Acoustics声学模块进行仿真计算,本文建立刚球四边形边界元网格,如图4所示,网格划分尺寸通常取波长的1/6,即
l= λ 6= 2 π f 6 c
式中,l为网格最大尺寸,λ为声波波长,f为频率,c为水中声速。
分别计算ka={1,2,3,4,5,6,7,8}时的散射声场,图5给出了刚球周围2 m范围内近场的声压干涉图。
图5 刚球声散射近场声压干涉图
可以看出,边界元算法能够对刚性球表面声压和球体周围声场分布给出较好的计算结果,不论球体的声学“照亮区”还是“影区”,算法均能给出准确的计算结果,同时,随着ka值即计算频率的增大,刚性球近场声压干涉条纹随之加密,能够较好地模拟水下目标近场声压分布。
为了与上文中的理论解析解进行比较,本文需要计算刚球指向性函数,为此利用软件自带的Field Point Meshes功能下的Directivity Field Point Mesh功能,在距球心50 m处每隔1°设置一个计算场点,用于提取该处的声压。将提取到的声压根据球面衰减规律归算到刚球等效声中心1 m处,并与解析解进行比较,如图6所示。可以看出,数值算法的计算结果与理论解析解吻合较好,能够准确地给出刚性球不同ka值下的指向性,证明该方法能够准确模拟任意分置角下水下目标收发分置声散射强度。
图6 远场声压指向性数值解与解析解比较

3 结束语

本文利用LMS Virtual.Lab软件,采用直接边界元法计算水下刚性球体在平面声波作用下的指向性函数,并与理论解析解进行比较,得到如下结论:
1)无论解析解还是边界元方法数值解,计算时需满足远场条件,否则,无法得到稳定的、吻合度较高的计算结果;
2)采用Virtual.Lab Acoustics声学模块计算远场声压时,当计算点输出声压数值较小时,需要提高软件参数设置中默认声压输出精度,否则,导出的计算结果数据精度将无法满足需要;
3)边界元法能够准确计算水下目标的表面声压分布、近场声场分布以及收发分置声散射强度,可用于水下目标声散射特性预报。
[1]
徐忠昌, 张明敏. 刚性球的散射方向性图仿真研究[J]. 海军工程大学学报, 2012, 24(6): 65-69.

[2]
Waterman P C. New Formulation of Acoustic Scattering[J]. Journal of the Acoustical Society of America, 1969, 45(6): 1417-1429.

DOI

[3]
汤渭霖, 范军. 水中弹性球壳的共振声辐射理论[J]. 声学学报, 2000, 25(4): 308-312.

[4]
宋延勇, 苏明旭, 蔡小舒, 等. 刚性颗粒声散射计算的边界元法实现[J]. 过程工程学报, 2009, 9(增2): 107-111.

[5]
李运思, 苏明旭, 宋延勇, 等. 单个椭球颗粒声散射数值研究[J]. 南京大学学报(自然科学), 2015, 51(6): 1139-1147.

[6]
詹福良, 徐俊伟. Virtual.Lab Acoustics 声学仿真计算从入门到精通[M]. 西安: 西北工业大学出版社, 2013.

Outlines

/