本发明涉及的是一种海洋信道下弹性结构辐射声场预报以及水声目标探测和信道海洋环境参数反演等领域下高效准确的研究方法。
背景技术:
弹性结构在海洋信道下的声辐射预报方法的研究,对开展结构振动辐射噪声实时预报和有效控制具有重要的理论研究价值,对海洋信道中辐射噪声的实时监测预报具有举足轻重的作用,是以后我国水声技术领域长期关注的热点和难点问题之一。
然而,目前弹性结构在海洋信道下的振动与声辐射研究的流体域大多考虑为无界或半空间流体域,对于海洋信道下弹性结构多物理场耦合声辐射问题的研究尚不多见。因为信道下结构声场问题涉及流固耦合、声壳耦合甚至声与复杂海底边界的声边界耦合等复杂多物理场耦合环境,在这些耦合条件下数学理论推导难以求解,数值模型无法建立、结构表面振动信息难以获取。传统数值法(边界元法bem、有限元法fem、统计能量法sea等)将受大网格计算量和复杂海洋环境和多物理场耦合等因素的严重限制,无法开展相关研究工作;解析解法只能分析部分简单结构在简单信道下的声场问题,对复杂海洋信道下任意弹性结构的声辐射研究也无能为力,也有部分学者提出直接忽略结构与流体、结构与环境的耦合作用,把结构视为点声源,但这样直接忽略了结构的近场声辐射特性。针对海洋信道环境的研究,通常采用海洋声传播理论(抛物线方程法,简正波法,波数积分法等)进行信道环境分析,抛物线方程法(pe)因为在进行远距离声场计算结果准确,计算速度快和适应性强等特点,受到了各国研究者的青睐,但研究对象通常为海洋信道下点源模型,对海洋信道下弹性结构声辐射研究很少涉及。以上共同导致了无法从理论解和数值法角度有效地开展信道下弹性结构声辐射研究,但它对海洋中水下结构声辐射、预报和识别有极为重要的意义,急需探索一种新的研究方法来解决海洋信道中结构声辐射问题。
技术实现要素:
本发明的目的在于提供一种对弹性结构和海洋信道适应能力强,结果准确效率高,使用简单易于推广的基于耦合fem-pe的海洋信道下弹性结构声辐射预报方法。
本发明的目的是这样实现的:
步骤一,建立海洋信道下弹性结构声辐射局域多物理场数值模型,采用有限元数值法分别对流固耦合、声边界耦合和无穷远处smerfield辐射条件进行声学定义,计算局部海洋信道流体域下的声场信息;
步骤二,进行三维信道空间下二维声场处理,通过远场场点与结构中心在竖直方向形成二维截面,并在距离结构中心一定距离处的深度方向上选取一维截线,提取该截线上的复声压信息作为声场外推抛物线方程法的初始场信息;
步骤三,对有限元计算提取的声场信息进行与声场外推抛物线方程法初始场的匹配,使有限元计算结果完整地匹配到抛物线方程法空间初始场的每个格点;
步骤四,设置与有限元相同的信道环境参数和场点位置,对信道下任意场点声压进行声场外推抛物线方程法的快速计算,得到海洋信道下结构声辐射特性,用于进行信道下弹性结构辐射声场快速预报。
本发明还可以包括:
1、所述在距离结构中心一定距离处的深度方向上选取一维截线中,一维截线水平距结构中心的距离为计算频率对应的波长即l≈λ。
2、采用分段三次艾尔米特插值技术对有限元法复声压与抛物线方程法声场进行空间格点的匹配,建立有限元法与抛物线方程法耦合条件,作为抛物线方程法有限差分法计算的初始场值。
本发明的理论模型:
如图1所示,本发明主要理论部分由两大部分组成:多物理场fem数值计算和pe声场计算,通过复声压线性插值作为fem与pe的耦合条件即pe声场计算初始条件。以下将说明两大部分的计算物理理论。
多物理场局域fem数值计算
海洋信道下弹性结构声辐射问题涉及流固耦合、声边界耦合和无限远边界等声学耦合边界处理,采用有限元建立多物理场下声学控制方程和边界连续条件如下:
流固耦合方程
在结构表面与外部流体接触的耦合面上,满足的边界条件为结构表面法向的振动速度与外部流体介质的振动速度相同,可写出结构与流体的耦合方程为
其中,刚度矩阵kij和阻尼矩阵cij、质量矩阵mij均为n×n阶矩阵。下标a、s和c为声学系统矩阵、力学结构矩阵和耦合矩阵;定义耦合矩阵kc、mc为
海洋信道的海面边界通常为dirichlet边界,满足的边界条件为界面声压为零
pa(x,y,z)|z=0=0(2)
其中,为了区分海水和海底,下标a和b分别标示海水流体层和海底层。
对于液态海底,满足的边界条件为声压p(x,y,z)连续,法向振速v(x,y,z)连续
pa(x,y,z)=pb(x,y,z)(3)
van(x,y,z)=vbn(x,y,z)(4)
对于各向同性弹性海底,满足法向上满足位移连续和应力连续,切向的应力为零
其中,u为弹性体中的水平位移、w为弹性体中的垂直位移,ρ为介质密度,σ为应力,λ、μ为拉米常数,定义δ为
海洋信道四周边界为无限远边界,有限元数值法采用pml(perfectlymatchedlayer,简称pml)技术进行模拟,pml通过对控制方程增加吸收系数转换为吸收层的控制方程,为了简化方程描述,令x轴为x1轴、y轴为x2轴,利用分离变量可写出频域下的pml方程
其中,σi为吸收系数,vi,pi为匹配层域的速度和声压幅值。
采用pml对边界进行处理后,使在边界上满足smerfield远场熄灭条件
p(x,y,z)|r=∞=0(9)
在边界上通过声吸收使声压为零,达到边界没有反射声以模拟信道四周的无限大空间。
通过流固耦合方程、声边界耦合方程和pml技术建立海洋信道下结构声辐射数值模型,计算获取结构辐射声场信息,便可提取从海面到底截线声压数据,为了使得fem计算结果与pe设置条件在空间格点完美匹配,需要通过分段三次艾尔米特插值法对包含幅值和相位信息声场数据进行处理,然后作为fem-pe耦合条件即作为pe法声场快速计算的初始条件。
海洋信道下pe法声场外推理论
推导抛物线方程的出发点为轴对称坐标系(r,z),简谐点源在水平变化声道中的声场由波动方程
其中,p为流体介质声压,ρ(r,z)为流体密度,c(r,z)为流体声速,zs为声源在z轴方向的位置。
在轴对称坐标系下,当流体介质密度为常数时,时间关系为e(-jωt)(ω为角频率)的简谐点源的亥姆霍兹方程为
其中,k=ω/c为波数
导出标准的抛物型波动方程有几种方法,本发明中严格按照tappert的方法,假设方程(11)的解形式为
p(r,z)=ψ(r,z)v(r)(12)
将(12)式代入(11)式,分离变量可得:
其中,k0参考波数,c0参考声速,且k(r,z)=k0n(r,z),n(r,z)=c0/c(r,z)。
在满足(k0r)>>1的前提条件下,式(13)的解可以用以下渐近式形式表示
接着做远场条件假设((k0r)>>1),可得到简化的椭圆形波动方程
为了求解以上方程,首先定义以下两个算符
根据定义的两个算符,可以把椭圆型波动方程写成如下形式
对式(18)进行因式分解,将方程分解为输入波和输出波两个分量,即
(p+ik0-ik0q)(p+ik0+ik0q)ψ-ik0(pq-qp)ψ=0(19)
当介质与距离无关即n≡n(z),算符p和q可以相互交换位置,因此式(19)最后一项为零,只考虑输出波分量,就得到
为了后续求解方便,这里令
从而式(17)给出的平方根算符可以写为
随着研究的不断深入,对根式(22)采取不同的近似处理,可以得到不同精度要求和不同传播角度的抛物方程,常用的处理方式是一种基于帕德级数展开的甚大角度范围的广义pe,这种根式的近似处理方法使抛物方程方法能够处理的传播角度几乎达到了90°,对算符q进行帕德级数展开,可得到
其中,
将式(23)带入方程(20),便可得到以水平算符的帕德级数展开式为基础的大角度抛物线方程
本发明采用claerbout得出的大角度pe,需用有限差分法(finitedifferencemethod,简称fdm)进行求解,数值求解过程中需要对流体域进行离散,海洋信道下离散各方向满足的步距为
式中:
可知,抛物线有限差分法是一个步进求解的过程,由前一个场信息便可求解下一个场的信息,采用该形式便可求解在初始距离r0处沿深度分布的声场。在声场求解过程中需要对海洋环境的初始条件和边界条件做出详细规定,通常,海面自由表面采用声场软边界处理,因而在表面满足ψ(r,0)=0,海底采用半无限均匀液态/弹性空间处理,把海底延续部分满足的辐射条件通过用厚度为几个波长的人工吸收层来限定。关于建立pe初始场数据的问题,可指定为试验测得或数值计算获取的具有一定指向性声源所产生的沿深度方向的声场(包括幅度和相位),本发明通过有限元数值法获取多物理场局域环境下复杂结构体声源在深度方向的声场,作为pe的计算的初始场条件。
根据有限元计算在深度方向提取声场结果p,可建立深度坐标zf与声压pf的关系,在pe计算域深度方向上水平距结构中心ro(ro=λ)处每个深度格点(r0,z)声场可通过以下分段三次艾尔米特插值法获取
其中,zp为最小区间[zk,zk+1]上一点,p'f,k为有限元提取结果函数在节点zk处的导数值,k=0,1,2,3···n,n为fem提取结果离散个数。
通过式(27)便可获取pe在深度方向的每个格点的声场值
本发明具有的突出优势为:
(a)该发明对弹性结构形状和海洋信道环境类型适应能力强,海洋信道下弹性结构多物理场局域模型采用有限元法建立,因此可对任意复杂结构在不同海洋信道下的辐射声场进行有限元数值计算,对结构类型具有很强的适应性。然后把局域下计算的声场信息作为pe法的声源输入条件,进行海洋信道下远场声场计算。通过推导标准化pe法可进行具有液态海底(含声吸收)的海洋信道下结构声辐射预报,若海洋海底为质地较硬的弹性层,则采用弹性pe对海底声场进行理论建立,求解在具有弹性层海底的海洋信道下弹性结构辐射声场,且对于浅海楔形海底或深海信道下结构声辐射预报问题也同样使用,所以该方法可对不同海洋信道下任意弹性结构辐射声场进行快速预报。
(b)该发明在进行计算时具有计算用时少、效率高的特点,通过对海洋信道下抛物线方程中平方根算符进行pade近似处理,加快了pe法声场计算的收敛和计算速度。在相同计算硬件条件下,可以快速的计算传统数值法(有限元法、边界元法、统计能量法等)无法模拟或因计算量大而无法计算的复杂声辐射问题,且计算结果准确高效。
(c)该发明操作过程简单使用方便,实施工作量小,易于在理论研究和实际工程中推广。其关键一步为只需通过数值法截线提取或通过试验法的垂直线列阵获取结构深度方向上线声压分布信息,进行pe法初始场完美匹配后,便可进行海洋信道下弹性结构的外场辐射声场快速预报。
与以往预报方法相比,该发明对弹性结构和海洋信道适应能力强,且计算结果准确效率高,使用简单易于推广。有效地解决目前在研究海洋信道下弹性结构声辐射预报时所遇到的计算量大、物理场多和信道环境复杂等诸多瓶颈性问题。
附图说明
图1a为本发明采用的海洋信道下弹性结构辐射声场预报模型。
图1b为本发明进行浅海信道下结构声场预报的原理示意图。
图2为浅海信道下点源声传播验证模型。
图3a为浅海信道下深度方向上20hz有限元与波数积分法计算结果对比图。
图4a为本发明方法在20z下计算结果与有限元计算结果对比图。
图4b为本发明方法在30hz下计算结果与有限元计算结果对比图。
图4c本发明方法在40hz下计算结果与有限元计算结果对比图。
图4d本发明方法在50hz下计算结果与有限元计算结果对比图。
图5a为采用本发明进行海洋浅海信道下弹性圆柱壳声辐射预报原理图
图5b为fem建立海洋浅海信道下弹性圆柱壳声辐射局域数值模型。
图6a为fem在频率为100hz下截面辐射声场空间分布。
图6b为fem在频率为200hz下截面辐射声场空间分布。
图6c为fem在频率为400hz下截面辐射声场空间分布。
图6d为fem在频率为800hz下截面辐射声场空间分布。
图7a为本发明计算在频率为100hz下辐射声压级随距离的变化曲线。
图7b为本发明计算在频率为200hz下辐射声压级随距离的变化曲线。
图7c为本发明计算在频率为400hz下辐射声压级随距离的变化曲线。
图7d为本发明计算在频率为800hz下辐射声压级随距离的变化曲线。
图8a为本发明计算在频率为100hz下辐射声场伪彩图。
图8b为本发明计算在频率为200hz下辐射声场伪彩图。
图8c为本发明计算在频率为400hz下辐射声场伪彩图。
图8d为本发明计算在频率为800hz下辐射声场伪彩图。
图9的表1为pekeris浅海信道环境和弹性结构参数。
图10的表2为本发明对不同计算频率的计算时间测试。
图11的表3为本发明对不同计算范围的计算时间测试。
具体实施方式
下面举例对本发明做更详细的描述。
1.多物理场局域环境下初始场获取
海洋信道下弹性结构声辐射问题涉及流固耦合、声边界耦合和无限远边界等声学耦合边界处理,采用有限元建立多物理场下声学控制方程和边界连续条件如下:
流固耦合方程
在结构表面与外部流体接触的耦合面上,满足的边界条件为结构表面法向的振动速度与外部流体介质的振动速度相同,可写出结构与流体的耦合方程为
其中,刚度矩阵kij和阻尼矩阵cij、质量矩阵mij均为n×n阶矩阵。下标a、s和c为声学系统矩阵、力学结构矩阵和耦合矩阵;定义耦合矩阵kc、mc为
海洋信道的海面边界通常为dirichlet边界,满足的边界条件为界面声压为零
pa(x,y,z)|z=0=0(2)
其中,为了区分海水和海底,下标a和b分别标示海水流体层和海底层。
对于液态海底,满足的边界条件为声压p(x,y,z)连续,法向振速v(x,y,z)连续
pa(x,y,z)=pb(x,y,z)(3)
van(x,y,z)=vbn(x,y,z)(4)
对于各向同性弹性海底,满足法向上满足位移连续和应力连续,切向的应力为零
其中,u为弹性体中的水平位移、w为弹性体中的垂直位移,ρ为介质密度,σ为应力,λ、μ为拉米常数,定义δ为
海洋信道四周边界为无限远边界,有限元数值法采用pml(perfectlymatchedlayer,简称pml)技术进行模拟,pml通过对控制方程增加吸收系数转换为吸收层的控制方程,为了简化方程描述,令x轴为x1轴、y轴为x2轴,利用分离变量可写出频域下的pml方程
其中,σi为吸收系数,vi,pi为匹配层域的速度和声压幅值。
采用pml对边界进行处理后,使在边界上满足smerfield远场熄灭条件
p(x,y,z)|r=∞=0(9)
在边界上通过声吸收使声压为零,达到边界没有反射声以模拟信道四周的无限大空间。
通过流-固耦合方程、声-边界耦合方程和pml技术建立海洋信道下结构声辐射数值模型,计算获取结构辐射声场信息。为了使得fem计算结果与pe设置条件在空间格点完美匹配,需要通过分段三次艾尔米特插值法对包含幅值和相位信息声场数据进行处理,然后作为fem-pe耦合条件即作为pe法声场快速计算的初始条件。
2.海洋信道下pe法声场外推理论
推导抛物线方程的出发点为轴对称坐标系(r,z),简谐点源在水平变化声道中的声场由波动方程
其中,p为流体介质声压,ρ(r,z)为密度,c(r,z)为声速,zs为声源在z轴方向的位置。
在轴对称坐标系下,当流体介质密度为常数时,时间关系为e(-jωt)(ω为角频率)的简谐点源的亥姆霍兹方程为
其中,k=ω/c为波数。
导出标准的抛物型波动方程有几种方法,本发明中严格按照tappert的方法,假设方程(11)的解形式为
p(r,z)=ψ(r,z)v(r)(12)
将(12)式代入(11)式,分离变量可得:
其中,k0参考波数,c0参考声速。且k(r,z)=k0n(r,z),n(r,z)=c0/c(r,z)。
在满足(k0r)>>1的前提条件下,式(13)的解可以用以下渐近式形式表示
接着做远场条件假设((k0r)>>1),可得到简化的椭圆形波动方程
为了求解以上方程,首先定义以下两个算符
根据定义的两个算符,可以把椭圆型波动方程写成如下形式
对式(18)进行因式分解,将方程分解为输入波和输出波两个分量,即
(p+ik0-ik0q)(p+ik0+ik0q)ψ-ik0(pq-qp)ψ=0(19)
当介质与距离无关即n≡n(z),算符p和q可以相互交换位置,因此式(19)最后一项为零,只考虑输出波分量,就得到
为了后续求解方便,这里令
从而式(17)给出的平方根算符可以写为
随着研究的不断深入,对根式(22)采取不同的近似处理,可以得到不同精度要求和不同传播角度的抛物方程,常用的处理方式是一种基于帕德级数展开的甚大角度范围的广义pe,这种根式的近似处理方法使抛物方程方法能够处理的传播角度几乎达到了90°,对算符q进行帕德级数展开,可得到
其中,
将式(23)带入方程(20),便可得到以水平算符的帕德级数展开式为基础的大角度抛物线方程
本发明采用claerbout得出的大角度pe,需用有限差分法(finitedifferencemethod,简称fdm)进行求解,数值求解过程中需要对流体域进行离散,海洋信道下离散各方向满足的步距为
式中:
可知,抛物线有限差分法是一个步进求解的过程,由前一个场信息便可求解下一个场的信息,采用该形式便可求解在初始距离r0处沿深度分布的声场。在声场求解过程中需要对海洋环境的初始条件和边界条件做出详细规定,通常,海面自由表面采用声场软边界处理,因而在表面满足ψ(r,0)=0,海底采用半无限均匀液态/弹性空间处理,把海底延续部分满足的辐射条件通过用厚度为几个波长的人工吸收层来限定。关于建立pe初始场数据的问题,可指定为试验测得或数值计算获取的具有一定指向性声源所产生的沿深度方向的声场(包括幅度和相位),本发明通过有限元数值法获取多物理场局域环境下复杂结构体声源在深度方向的声场,作为pe的计算的初始场条件。
根据有限元计算在深度方向提取声场结果p,可建立深度坐标zf与声压pf的关系,在pe计算域深度方向上水平距结构中心ro(ro=λ)处每个深度格点(r0,z)声场可通过以下分段三次艾尔米特插值法获取
其中,zp为最小区间[zk,zk+1]上一点,p'f,k为有限元提取结果函数在节点zk处的导数值,k=0,1,2,3···n,n为fem提取结果离散个数。
通过式(27)便可获取pe在深度方向的每个格点的声场值
方法准确性验证
算例1:pekeris浅海信道下点源声传播计算
建立如图2所示的轴对称条件下pekeris信道下点源声传播模型,信道深度h=200m,信道上界面为dirichlet边界,下界面为无限半空间等声速液态层,海水参数:海水密度为ρa=1024kg/m3,海水声速ca=1500m/s;海底密度ρb=1800kg/m3,海底声速cb=2500m/s。声源为单极点源,单极幅值为1,点源深度为z0=100m。采用多物理场有限元数值计算点源在信道环境下的声场传播,计算距离为r=2000m,并提取距离点源一个波长λ处有限元(fem)在整个浅海深度方向上的声压信息(幅值与相位),作为fem-pe耦合条件,即pe声场计算的初始条件,设置pe计算浅海信道环境参数后,便可计算在浅海信道下任意场点的声场信息。
如图3a所示,通过多物理场耦合fem数值法计算获取在距离点源为λ=75m处深度方向上频率为20hz的复声压值,并画出在深度方向上的声压模值分布曲线,并与相同信道条件下波数积分法计算程序ffp计算结果进行对比。然后把fem在深度方向准确的计算结果进行分段三次艾尔米特插值处理后,计算了频率为20hz所对应的声传播曲线(各场点深度为50m),并与相同条件下fem计算结果进行了对比,如图4a所示。同样,计算了30hz、40hz和50hz条件下结构辐射声场的辐射声压级随传播距离曲线图,并与有限元计算结果进行了对比,如图4b、4c和4d所示。可以看出20hz、30hz、40hz和50hz频率下本发明方法计算结果与fem计算结果吻合性很好,验证了本发明声场计算的准确性。
方法实用性和高效性说明
算例2:pekeris浅海信道下弹性结构声辐射计算
如图5a所示,建立pekeris浅海信道下弹性圆柱壳辐射声场fem-pe计算理论模型,该计算理论模型主要分为数值fem计算域和pe计算域,提取数值fem计算域在深度方向上的复声压信息并做插值后作为fem-pe耦合条件。其中s为弹性结构声源,zs结构中心在z轴上的位置,ρa,ca分别为海水密度声速,ρb,cb分别为无限大液态层密度和声速,pl为fem-pe耦合场,λ为pl距离结构中心一个波长距离,rp为pe计算域下任意场点距离初始场的距离。
如图5b所示,采用多物理场耦合数值法建立pekeris浅海信道下弹性圆柱壳声辐射模型,信道环境和弹性结构参数见图9的表1,利用有限元法计算局域流体下结构声辐射信息,如图6a-图6d所示,数值计算了pekeris浅海信道局域流体下弹性圆柱壳在z-o-y平面内100hz、200hz、400hz和800hz频率下所对应的声场分布。并提取了距离圆柱壳中心一个波长处截线上辐射声压(幅值与相位)信息,为了与pe离散格点(δr,δz)进行空间完美匹配,需要对fem计算提取的复声压进行分段三次艾尔米特插值后,作为fem与pe耦合条件即pe声场计算初始条件,然后采用声场pe外推计算了远场条件下弹性结构辐射声场信息,如图7a到图7d,分别计算了100hz、200hz、400hz和800hz频率下弹性圆柱壳辐射声压级随距离的变化曲线,各场点深度为3m;如图8a到图8d,分别计算了100hz、200hz、400hz和800hz频率下弹性圆柱壳在远场辐射声场分布伪彩图。可看出,浅海信道下弹性结构在远场辐射声场基本满足柱面波传播规律,其声场波动细节随着频率的上升而增加。
为了说明本发明在进行海洋信道下弹性结构声辐射计算速度方面的巨大优势,分别对不同距离和不同频率下声场计算时间进行测试,如图10的表2和图11的表3所示,测试的流体环境和弹性结构与上述相同,测试硬件条件参数lenovothinkstatione30(英特尔八核3.2ghzcpu,32gb内存)。
从图10的表2和图11的表3可看出,不同距离对计算时间的影响比不同频率对计算时间的影响显著,当计算距离为0-10km时,计算频率从100hz到800hz变化过程中,计算时间略小增加,因为海洋信道下离散各方向满足的步距为