一种水下声辐射预报方法、系统、计算机设备、存储介质

1.本发明属于水下声辐射预测技术领域，尤其涉及一种水下声辐射预报方法、系统、计算机设备、存储介质。


背景技术：

2.目前：水下声问题是近几十年来非常活跃的研究领域。然而，大多数问题无法解析求解。因此，各种数值方法或技术被提出用于处理水下声复杂问题。到目前为止，基于网格的方法仍然是计算水下声辐射的主要手段，如边界元法(bem)和有限元法(fem)。一般来说，声学问题分为两类：内部声学问题和外部声学问题。边界元法和有限元法均可直接用于求解内部声学问题。对于外部声问题，边界元法采用自由空间格林函数作为基本解，在无穷远处自动满足sommerfeld辐射条件，是求解无限域声学问题的有效工具。然而，边界元法在分析大型问题时，由于系统矩阵为非对称性和非稀疏性，导致计算效率低。在过去的几十年中，边界元法的取得了显著的进展，这推动了边界远法在计算声学中的应用。例如，采用快速多极边界元法可显著提高计算效率。关于边界元法的更详细的介绍见参考文献。目前，声学边界元法仍在不断改进完善中。
3.有限元法中的系统矩阵为稀疏阵，使得有限元法已成为求解大型问题的一种非常有效的方法。然而，标准有限元法不能直接用于处理无限域的外部声学问题，需要对原始问题域进行特殊处理。通常通过人工边界截断无限域来产生有限计算域，并在人工边界上施加无反射边界条件，以保证声波在向外行进过程中自由衰减。dirichlet
‑
to
‑
neumann(dtn)边界条件、完全匹配层、无限元技术和吸收边界条件均可作为无反射边界条件。对于二维外部声问题，提出的r
‑
1/2
衰减特性的映射无限元可有效模拟声波的自由衰减过程，并给出了相应的数值积分形式。采用这些改进的映射无限元的一个优点是不需要人为移动原点。与具有指数衰减的无限元和其他阻尼单元相比，使用该类型无限元可获得更高的精度。在标准有限元法中，离散模型的刚度总比原始模型更硬。因此，数值色散误差不可避免，且该误差随着声波数的升高而急剧增加。然而，细化网格并不能有效降低色散误差。为抑制这种误差，近年来提出了一些基于光滑技术或无网格的新方法来软化有限元模型的“过硬”的刚度，如基于g空间理论的弱式，基于梯度光滑技术和w2形式的方法，包括混合光滑有限元法(hs
‑
fem)，单元基光滑α径向点插值法(cs
‑
αrpim)，边基光滑有限元法(esfem)、稳定点基光滑有限元方法(sns
‑
fem)和超收敛α有限元法。借助于这些新方法或技术，数值色散误差得以显著减小。
4.用无穷傅里叶级数解析表示的dtn边界条件是一种精确的无反射边界条件，该边界条件通过积分算子建立dirichlet变量与neumann变量解析关系。虽然dtn边界条件是非局部边界条件，但对于大型声学计算问题，计算成本可忽略。然而，在实际应用dtn边界时，必须将无穷级数截断为有限项数n。当n不够大时，无法保证外部声学问题解的唯一性和可解性。因此，截断形式的dtn边界条件不再精确。虽然可以包含足够的项数来确保唯一性和可解性，但计算量过大，尤其是对高声波数或大尺度人工边界。为了克服这些困难，提出了
改进的dtn(mdtn)边界条件，对于任意的n值，可确保声场求解的唯一性和可解性。mdtn边界条件通过在标准的dtn映射中添加一个局部微分算子，同时减去该算子来表达。与截断的dtn边界条件相比，mdtn边界条件使用少量的项数即可获得更高的精度。而在声场预报中事先不知道须包含多少项数才能达到所期望的精度，因此mdtn边界条件在实际计算中更具优势。mdtn边界条件已应用于实际声学计算问题和处理其他类型问题。
5.在数值计算领域，与基于网格的方法相比，无网格法是相对新兴的方法。无网格法采用任意分布的节点来表示问题域及其边界，无需使用节点连接信息或网格来进行变量插值或近似。根据公式类型，无网格方法一般分为三类：如无网格强式法(如：无网格配点法)，无网格弱式法(如：无单元伽辽金法，efg)和无网格强
‑
弱(mws)式法。近年来，无网格法被广泛应用于处理一些具有挑战性的数值和工程问题，并取得了巨大的成功。径向基函数(rbf)法、移动最小二乘法(mlsa)和单位分解法是无网格方法中常用的形函数构造方法。
6.无网格强式法是真正的无网格方法，该类方法变量插值和数值积分均不需要网格。但是，无网格强式法对具有导数边界条件的问题往往是不精确，不具备稳定和鲁棒性。无网格弱式法对近似函数的一致性要求较弱，需要使用全局或局部背景单元进行数值积分。与强式法相比，能够提供更稳定的离散系统方程，获得更精确的结果。基于全局galerkin弱式的无网格点插值方法(pim)是一种无网格全局弱式法。在pim中，采用多项式点插值法(ppim)在计算点的局部支持域中使用一组节点构造形函数，该形函数具有kronecker
‑
delta函数性质，可以像有限元法一样直接施加本质边界条件可。然而，ppim形成的矩矩阵可能是奇异的。使用两阶矩阵三角化算法可以消除该奇异问题，但由于ppim形状函数的不相容性，使得pim对不规则分布节点的鲁棒性较差。此外，当支持域中含有过多的节点时，会产生过高阶的插值多项式，导致ppim形函数剧烈振荡。为消除这种奇异性问题，基于径向基函数(rbf)提出了径向点插值方法(rpim)。rpim对任意分布的节点具有良好的稳定性和鲁棒性，具有比fem更快的收敛性。得益于rpim的稳定性性和有效性，该方法已被广泛应用于处理诸多数值和工程问题，并取得了巨大的成功，如智能材料问题、二维和三维固体力学问题、板壳结构问题，土木工程中的材料非线性问题等。
7.在无网格全局弱式法中，需要生成全局背景单元进行数值积分，这类全局弱式法不是真正的无网格方法。基于局部弱式，无网格局部petrov
‑
galerkin(mlpg)方法被提出，数值积分只在场节点对应的局部域上进行，只生成对应的局部背景单元，不使用全局背景单元。mlpg与无网格强式法非常相似，实现过程与无网格强式法一样简单。因此，mlpg可以说是一种真正的无网格方法。在mlpg中，采用mlsa构造形函数，导致形状函数不具备kronecker delta函数性质，需要采用特殊的技术来施加本质边界条件。基于mlpg，提出了一些改进型方法。然而，局部弱式法的这些优点是以一定代价换取的。在mlpg中，引入了附加参数，形成的系统矩阵通常是非对称阵，局部积分区域尺寸和试函数形式对计算精度影响明显；与efg相比，由于mlpg矩阵为非对称阵，计算耗时更长。
8.为了解决节点积分的困难，积分局部无网格(ilmf)方法被提出，以提高二维弹性力学和线弹性断裂力学求解的精度和效率。研究表明，应用ilmf方法，节点紧支撑域的大小和局部积分域分别直接决定计算精度和效率，这两个参数可以通过多目标优化自动获得，使得该类无网格局部法具有极佳的可靠性和鲁棒性。除了这些域基无网格方法外，边界型无网格方法被提出，该方法仅使用一组节点表示问题域的边界，其已用于处理许多具有挑
战性的问题，如边界型无网格方法已被应用于求解二维非线性问题。在边界型无网格方法中，通常需要通过数值实验预先确定紧支撑域的尺寸，并且问题域的边界不能自动离散。为了优化确定相应的参数和自动离散问题域边界，一种能够自动生成pareto最优节点排列和确定紧支撑域尺寸的多目标优化策略被提出，显著提升了边界型无网格方法的求解精度、稳定性和计算效率。无网格局部弱式法和边界型无网格方法仍在改进和发展中。
9.rpim已被用于室内声学计算问题，可以产生比fem更精确的结果。如果直接用rpim来处理自由场声辐射问题，则需要一组分布的场节点来表示无限大问题域。除此之外，还必须生成一组背景单元来进行积分运算。因此，无法直接采用rpim来实现自由场水下声辐射预报。一种无网格径向点插值与无限声波包络元法(wee)相结合的自由场声辐射预测方法被提出，该方法比有限元方案具有更高的计算精度。然而，该方法严格意义上不是无网格弱式方法，而且必须通过复杂的数学处理来建立rpim与wee之间的耦合关系。考虑到rpim和mdtn边界条件的各自优点，提出了一种基于rpim和mdtn边界条件相耦合的无网格弱式方法来精确预报无限域水下声辐射。
10.通过上述分析，现有技术存在的问题及缺陷为：现有的基于网格方法的水下声辐射预测精度不高，且计算效率不高。
11.解决以上问题及缺陷的难度为：精确计算无限域水下声辐射首先需对无限域进行截断，对人工边界进行处理，防止产生反射声波，确保声波自由衰减，满足无穷远处的辐射条件，获取等效的有限计算域；现有的基于网格的数值方法和技术计算水下声辐射时，在中高频段产生数值色散误差，如何有效抑制这种随频率升高而急剧增加的色散误差一直没有得到很好地解决；采用基于网格的方法计算水下声辐射，网格质量对计算精度影响很大，导致网格前处理过程繁琐、耗时长，需采用新兴的不依赖网格的数值方法来克服此问题；无网格法作为一种新兴的数值算法，在无限域水下声辐射计算领域还未得到应用，采用无网格方法来计算水下声辐射具有相当的挑战和难度。
12.解决以上问题及缺陷的意义为：基于伽辽金全局弱式，将径向基点插值无网格法与mdtn边界条件耦合，使得无限域水下声辐射计算问题等效为有限域声学计算问题，在有限计算域内采用径向基点插值无网格法提高了插值精度，并能有效抑制中高频色散误差，同时消除了采用dtn边界条件引起的声学非唯一性问题。该发明方法不依赖网格和场节点的连接信息进行声场量计算，消除了传统基于网格方法中网格质量对计算精度的影响，减少了前处理耗时；采用径向基点插值法构造声形函数，能同时提高插值精度和抑制中高频色散误差，全面提升了声场计算精度。因此，采用无网格方法来计算水下声辐射是一种思路创新和方法创新，突破了现有基于网格方法的理论局限。


技术实现要素：

13.针对现有技术存在的问题，本发明提供了一种水下声辐射预报方法、系统、计算机设备、存储介质。
14.本发明是这样实现的，一种水下声辐射预报方法，所述水下声辐射预报方法包括：
15.基于径向点插值法和改进的dirichlet
‑
to
‑
neumann边界条件进行纯无网格弱式的水下结构声辐射预报。
16.进一步，所述水下声辐射预报方法包括以下步骤：
17.步骤一，采用人工边界截断结构外部无限问题域获得有限计算域，将无限域水下声场计算问题等效转换为有限空间声场计算问题；
18.步骤二，在rpim
‑
mdtn中，采用rpim构造声形函数，在人工边界处施加mdtn边界条件，进行水下声辐射预报，联合rpim与mdtn消除人工边界处的虚假反射声波，确保声波向外自由衰减，提高插值精度和抑制中高频色散误差。
19.进一步，所述采用rpim构造声形函数包括：
[0020][0021]
n
t
(x)＝[n1(x),n2(x),
…
,n
n
(x)]；
[0022]
其中，n
t
(x)表示采用rpim构造的n个场节点的形函数；p
s
＝[p1,p2,
…
,p
n
]
t
表示声压值组成的向量；p
h
(x)表示声压p(x)的试函数；p
t
＝[p1(x),p2(x),
…
,p
m
(x)]。
[0023]
进一步，所述在rpim
‑
mdtn中，采用rpim构造声形函数，在人工边界处施加mdtn边界条件，进行水下声辐射预报还包括：
[0024]
在加权残差公式中，采用rpim构造的形函数n(x)作为检验函数w(x)，则rpim构造声形函数表示为：
[0025][0026]
建立如下矩阵形式的离散系统方程
[0027]
[k
‑
k2m+k
b
]{p}＝f；
[0028][0029][0030][0031][0032]
式中：{p}表示声压值的矢量，并且c
n
(x)＝[cosnθsinnθ]；
[0033]
基于得到的人工边界处的声压值，利用下式计算得到v
o
中任意场点处的声压值：
[0034][0035]
其中，h
n(1)
表示第一类n阶的汉克尔函数，p(r,θ)表示γ
r
上的声压，求和符号后的“'”表示对应n＝0的项乘0.5。
[0036]
进一步，所述mdtn边界条件包括：
[0037][0038]
其中，
[0039]
其中，上“'”表示对kr求导。
[0040]
本发明的另一目的在于提供一种接收用户输入程序存储介质，所存储的计算机程序使电子设备执行所述水下声辐射预报方法包括下列步骤：
[0041]
步骤一，采用人工边界截断结构外部无限问题域获得有限计算域；
[0042]
步骤二，在rpim
‑
mdtn中，采用rpim构造声形函数，在人工边界处施加mdtn边界条件，进行水下声辐射预报。
[0043]
本发明的另一目的在于提供一种存储在计算机可读介质上的计算机程序产品，包括计算机可读程序，供于电子装置上执行时，提供用户输入接口以实施所述水下声辐射预报方法。
[0044]
本发明的另一目的在于提供一种实施所述水下声辐射预报方法的水下声辐射预报系统，所述水下声辐射预报系统包括：
[0045]
有限计算域获得模块，用于采用人工边界截断结构外部无限问题域获得有限计算域；
[0046]
水下声辐射预报模块，用于在rpim
‑
mdtn中，采用rpim构造声形函数，在人工边界处施加mdtn边界条件，进行水下声辐射预报。
[0047]
结合上述的所有技术方案，本发明所具备的优点及积极效果为：本发明基于径向点插值法和改进的dirichlet
‑
to
‑
neumann边界条件提出了一种预报水下结构声辐射的纯无网格弱式方法，相比传统基于网格的方法，本发明方法不需要使用网格和场节点连接信息来计算声场量，能有效抑制中高频数值色散误差，消除了网格质量对计算精度的影响，全面提升了水下声场计算精度，同时简化了前处理过程，该发明方法突破了现有基于网格方法的局限。通过数值算例研究了该方法的性能。主要结论总结如下：
[0048]
1.rpim
‑
mdtn方法能够精确且稳定地计算无限域水下声辐射，而无需使用网格或节点连接信息来进行场变量插值；
[0049]
2.与传统的有限元方法相比，rpim
‑
mdtn方法可以产生更精确的结果，并且在收敛性方面表现得更好。另外，本发明受声波数的影响小。
[0050]
3.由于prim形状函数构造的原因，在相同配置下，rpim
‑
mdtn方法比标准的fem需要更多的计算时间。而当考虑计算精度时，本发明方法实际上具有更高的计算效率；
[0051]
4.相比exp
‑
rbf，mq
‑
rbf和tps
‑
rbf具有更高的精度。tps
‑
rbf具有稳定的收敛性，而exp
‑
rbf和mq
‑
rbf的收敛过程不稳定。
[0052]
5.通过数值试验获得了各参数的可用取值，为水下声辐射预报提供了数据参数支撑。
附图说明
[0053]
图1是本发明实施例提供的水下声辐射预报方法流程图。
[0054]
图2是本发明实施例提供的使用一组离散分布场节点表示的有限计算域示意图。
[0055]
图3是本发明实施例提供的支持域与影响域之间的区别示意图。
[0056]
图4是本发明实施例提供的浸没式无限长圆柱周向简谐声辐射模型示意图。
[0057]
图5是本发明实施例提供的使用任意分布的节点表示的计算域及用于积分的背景单元示意图。
[0058]
图6(a)是本发明实施例提供的第二阶辐射模式下ε
real
随α
s
的变化示意图。
[0059]
图6(b)是本发明实施例提供的第二阶辐射模式下ε
imag
随α
s
的变化示意图。
[0060]
图7(a)是本发明实施例提供的第三阶辐射模式下ε
real
随α
s
的变化示意图。
[0061]
图7(b)是本发明实施例提供的第三阶辐射模式下ε
imag
随α
s
的变化示意图。
[0062]
图8(a)是本发明实施例提供的对应mq
‑
rbf的ε
real
随α
c
的变化示意图。
[0063]
图8(b)是本发明实施例提供的对应mq
‑
rbf的ε
imag
随α
c
的变化示意图。
[0064]
图9(a)是本发明实施例提供的第二阶辐射模式下对应mq
‑
rbf的ε
real
随q的变化示意图。
[0065]
图9(b)是本发明实施例提供的第二阶辐射模式下对应mq
‑
rbf的ε
imag
随q的变化示意图。
[0066]
图9(c)是本发明实施例提供的第三阶辐射模式下对应mq
‑
rbf的ε
real
随q的变化示意图。
[0067]
图9(d)是本发明实施例提供的第三阶辐射模式下对应mq
‑
rbf的ε
imag
随q的变化示意图。
[0068]
图10(a)是本发明实施例提供的对应exp
‑
rbf的ε
real
随α
c
的变化示意图。
[0069]
图10(b)是本发明实施例提供的对应exp
‑
rbf的ε
imag
随α
c
的变化示意图。
[0070]
图11(a)是本发明实施例提供的第二阶辐射模式对应tps
‑
rbf的ε
real
随η的变化图。
[0071]
图11(b)是本发明实施例提供的第二阶辐射模式对应tps
‑
rbf的ε
imag
随η的变化示意图。
[0072]
图11(c)是本发明实施例提供的第三阶辐射模式对应tps
‑
rbf的ε
real
随η的变化图。
[0073]
图11(d)是本发明实施例提供的第三阶辐射模式对应tps
‑
rbf的ε
imag
随η的变化示意图。
[0074]
图12(a)是本发明实施例提供的与mdtn和dtn边界条件对应的误差随n的变化图。
[0075]
图12(b)是本发明实施例提供的与mdtn和dtn边界条件对应的误差随n的变化示意图。
[0076]
图13(a)是本发明实施例提供的第四阶辐射模式下在ka＝π/2时的声压幅频响应示意图。
[0077]
图13(b)是本发明实施例提供的第五阶辐射模式下在ka＝π/2时的声压幅频响应示意图。
[0078]
图13(c)是本发明实施例提供的第六阶辐射模式下在ka＝π/2时的声压幅频响应示意图。
[0079]
图14(a)是本发明实施例提供的第五阶辐射模式下观测点处的声压幅频响应图。
[0080]
图14(b)是本发明实施例提供的第五阶辐射模式下观测点处的声压幅频响应示意图。
[0081]
图15(a)是本发明实施例提供的方法与fem的收敛性比较图。
[0082]
图15(b)是本发明实施例提供的方法与fem的收敛性比较示意图。
[0083]
图16(a)是本发明实施例提供的方法与fem之间的cpu时间随平均节点间距的变化
示意图。
[0084]
图16(b)是本发明实施例提供的方法与fem之间的cpu时间随声压实部误差的变化示意图。
[0085]
图17是本发明实施例提供的浸没无限长圆柱扇形声辐射示意图。
[0086]
图18(a)是本发明实施例提供的k＝4π时人工边界处的声压分布图。
[0087]
图18(b)是本发明实施例提供的k＝4π时人工边界处的声压分布示意图。
[0088]
图18(c)是本发明实施例提供的k＝8π时人工边界处的声压分布图。
[0089]
图18(d)是本发明实施例提供的k＝8π时人工边界处的声压分布示意图。
[0090]
图19是本发明实施例提供的刚性舵形物体的几何形状示意图。
[0091]
图20(a)是本发明实施例提供的kl＝2π时人工边界处的声压分布图。
[0092]
图20(b)是本发明实施例提供的kl＝2π时人工边界处的声压分布示意图。
[0093]
图21(a)是本发明实施例提供的kl＝4π时人工边界处的声压分布图。
[0094]
图21(b)是本发明实施例提供的kl＝4π时人工边界处的声压分布示意图。
[0095]
图22(a)是本发明实施例提供的kl＝8π时人工边界处的声压分布图。
[0096]
图22(b)是本发明实施例提供的kl＝8π时人工边界处的声压分布示意图。
[0097]
图23是本发明实施例提供的二维海底形结构的几何形状和尺寸示意图。
[0098]
图24(a)是本发明实施例提供的f＝100hz时人工边界上的声压分布图。
[0099]
图24(b)是本发明实施例提供的f＝100hz时人工边界上的声压分布示意图。
[0100]
图25(a)是本发明实施例提供的f＝200hz时人工边界上的声压分布图。
[0101]
图25(b)是本发明实施例提供的f＝200hz时人工边界上的声压分布示意图。
[0102]
图26(a)是本发明实施例提供的f＝400hz时人工边界上的声压分布图。
[0103]
图26(b)是本发明实施例提供的f＝400hz时人工边界上的声压分布示意图。
[0104]
图27(a)是本发明实施例提供的观测点a处的声压频率响应图。
[0105]
图27(b)是本发明实施例提供的观测点a处的声压频率响应示意图。
[0106]
图28(a)是本发明实施例提供的观测点b处的声压频率响应图。
[0107]
图28(b)是本发明实施例提供的观测点b处的声压频率响应示意图。
具体实施方式
[0108]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0109]
针对现有技术存在的问题，本发明提供了一种水下声辐射预报方法，下面结合附图对本发明作详细的描述。
[0110]
本发明实施例提供的水下声辐射预报方法包括：
[0111]
基于径向点插值法和改进的dirichlet
‑
to
‑
neumann边界条件进行纯无网格弱式的水下结构声辐射预报。
[0112]
如图1所示，本发明实施例提供的水下声辐射预报方法包括以下步骤：
[0113]
s101，采用人工边界截断结构外部无限问题域获得有限计算域；
[0114]
s102，在rpim
‑
mdtn中，采用rpim构造声形函数，在人工边界处施加mdtn边界条件，
进行水下声辐射预报。
[0115]
本发明实施例提供的采用rpim构造声形函数包括：
[0116][0117]
n
t
(x)＝[n1(x),n2(x),
…
,n
n
(x)]；
[0118]
其中，n
t
(x)表示采用rpim构造的n个场节点的形函数；p
s
＝[p1,p2,
…
,p
n
]
t
表示声压值组成的向量；p
h
(x)表示声压p(x)的试函数；p
t
＝[p1(x),p2(x),
…
,p
m
(x)]。
[0119]
本发明实施例提供的在rpim
‑
mdtn中，采用rpim构造声形函数，在人工边界处施加mdtn边界条件，进行水下声辐射预报还包括：
[0120]
在加权残存公式中，采用rpim构造的形函数n(x)作为检验函数w(x)，则rpim构造声形函数表示为：
[0121][0122]
建立如下矩阵形式的离散系统方程
[0123]
[k
‑
k2m+k
b
]{p}＝f；
[0124][0125][0126][0127][0128]
式中：{p}表示声压值的矢量，并且c
n
(x)＝[cosnθsinnθ]；
[0129]
基于得到的人工边界处的声压值，利用下式计算得到v
o
中任意场点处的声压值：
[0130][0131]
其中，h
n(1)
表示第一类n阶的汉克尔函数，p(r,θ)表示γ
r
上的声压，求和符号后的“'”表示对应n＝0的项乘0.5。
[0132]
本发明实施例提供的mdtn边界条件包括：
[0133][0134]
其中，
[0135]
其中，上“'”表示对kr求导。
[0136]
下面结合具体实施例对本发明的技术方案做进一步说明。
[0137]
实施例1：
[0138]
1.水下声辐射计算理论背景
[0139]
1.1.无限域水下声辐射计算理论
[0140]
在无限大自由场v中，结构做稳态简谐振动辐射的声压p(x)满足如下形式的helmholtz方程
[0141][0142]
式中：表示拉普拉斯算子，k＝ω/c为声波波数，ω为结构边界振动圆频率，c为声音在流体介质中的传播速度。将结构边界γ分为dirichlet型边界γ
d
和neumann边界γ
n
，γ
d
和γ
n
满足γ
d
∪γ
n
＝γ和两种类型边界定义如下：
[0143]
p＝p
d
onγ
d
(2)
[0144][0145]
式中：p＝p
d
为边界γ
d
的已知声压，v(x)为边界γ
n
上给定的法向振速，n
b
为边界单位外法向，ρ为流体介质密度，i2＝
‑
1。
[0146]
为确保声波向外传播过程中自由衰减和声场求解的唯一性，声压应满足sommerfeld辐射条件
[0147][0148]
式中：r为场点到坐标原点的距离；β为维度参数，对二维和三维问题分别取1和2。若p(x,t)＝p(x)exp(
‑
iωt)，则通常只用第二项表示sommerfeld辐射条件，在无穷远处，此辐射条件与ρc阻抗条件等效，即
[0149]
1.2.修正型dirichletto neumann边界条件
[0150]
为了获得有限计算域，采用人工边界γ
r
(对应2d或3d情况分别为半径为r的圆或球)截断无限大外部问题域，如图2所示。有限计算域以γ为内边界，以γ
r
为外边界。在人工边界上施加dtn边界条件式中m称为dtn映射算子，其作用是在边界γ
r
上建立dirichlet量声压p与neumann量的关系。从而，最初的声场求解问题可以等效为如下两个声场求解问题。
[0151][0152][0153]
为了完整起见，这里简要介绍dtn边界条件理论，更为详细的理论推导可以在参考文献。对二维声问题，采用极坐标(r,θ)可解析求解出方程(6)对应问题的解
[0154][0155]
式中：h
n(1)
表示第一类n阶的汉克尔函数，p(r,θ)为γ
r
上的声压，求和符号后的“'”表示对应n＝0的项乘0.5。求方程(7)关于r的偏微分，并取r＝r得到：
[0156][0157]
式中
[0158][0159]
其中h
n(1)
(kr)上“'”表示对kr求导。dtn内核m
n
(θ
‑
θ')可表示为
[0160][0161]
在计算中，将式(8)中的无穷级数截断为n个有限项，从而式(8)变为：
[0162][0163]
式中：m
n
表示截断后的dtn映射。从而，截断后的dtn边界条件不在精确。当n取值较小时，无法确保问题的可解性和解的唯一性。为克服这个缺陷并提高求解精度，mdtn边界条件被提出，以确保n取任意值时具有唯一的解。通过在式(8)右边加上和减去bp(r,θ)，同时将减去项bp(r,θ)求和并截断，得到如下mdtn边界条件
[0164][0165]
式中：代表mdtn映射，b代表局部微分算子。关于mdtn边界条件的更为详细的介绍可参考文献。本发明采用如下局部微分算子。
[0166][0167]
从而，mdtn边界条件可确定为：
[0168][0169]
其中
[0170][0171]
2.所提方法的构建公式
[0172]
2.1.形函数构造
[0173]
采用p
h
(x)表示声压p(x)的试函数，在v
i
中运用加权残存公式和格林理论得到如下弱式方程：
[0174][0175]
式中：w(x)表示加权残存公式中的检验函数。用散布在v
i
及其边界上的一组离散节点表示该有限计算域及其边界。v
i
中任意场点x处的声压可使用场节点声压值插值确定。
[0176]
采用rpim构造声压形函数，带有多项式的rpim插值可表示为：
[0177][0178]
式中：r
i
(x)表示径向基函数(rbf)，p
j
(x)为基于空间坐标x＝[x,y,z]
t
表示的单项式，n和m分别为rbf和多项式基函数的项数，a＝[a1,a2,
…
,a
n
]
t
和b＝[b1,b2,
…
,b
m
]
t
为未知的常数系数向量
[0179]
r
t
＝[r1(x),r2(x),
…
,r
n
(x)]
ꢀꢀꢀ
(18)
[0180]
p
t
＝[p1(x),p2(x),
…
,p
m
(x)]
ꢀꢀꢀ
(19)
[0181]
在式(17)中，并不总是使用多项式项，当m＝0时，只使用rbf基函数构造形函数。通常可以通过添加多项式项来改善rpim的性质，研究表明：使用这些多项式可以提高插值的稳定性，获得更精确的结果，并降低对形状参数的灵敏性，在更宽范围内选择形状参数取值。在rpim中，multi
‑
quadrics(mq)，gaussian(exp)和thinplate spline(tps)是最为常用的基函数：
[0182][0183]
式中：r
i
表示2d问题的节点x
i
与计算点x之间的距离
[0184][0185]
d表示x点附近的节点平均间距，α
c
≥0，q及η为rbf中的形状参数。由以上公式可知，r
i
(x)是r
i
的函数。在计算时。需要预先确定这些形状参数的取值，通常可通过数值实验获得这些形状参数的取值。
[0186]
为求解式(17)中的未知系数a和b，在计算点x的支持域ω
s
中的n个场节点处满足式(17)来构造n个方程，并使用m个附加方程
[0187][0188]
得到如下采用矩阵形式表示的系统方程
[0189][0190]
式中：p
s
＝[p1,p2,
…
,p
n
]
t
表示声压值组成的向量，r为
[0191][0192]
p
m
为
[0193][0194]
对于任意分布的节点，r
‑1通常存在，式(17)通常使用低阶多项式，则rpim中不会出现任何奇异问题。求解式(23)得到系数a和b。
[0195][0196]
将式(26)代入式(17)得到
[0197][0198]
式中：n
t
(x)＝[n1(x),n2(x),
…
,n
n
(x)]为采用rpim构造的n个场节点的形函数。
[0199]
在rpim形函数的构造中，局部支持域的确定至关重要。通常，采用影响域来确定计算点x的支持域包含的节点。影响域ω
f
被定义为一个场节点的影响区域，影响域的中心为该场节点。如计算点位于某场节点的影响域中，则该场节点属于参与构造该计算点的形函数，如图3所示。在后续部分，为了简化起见，使用圆形影响域。根据建议，采用r
w
＝α
s
d确定影响域半径r
w
，式中α
s
表示场节点x
i
的影响域的无量纲尺寸，该参数对计算精度至关重要。对于给定类型的问题，可以通过数值试验来确定α
s
取值。
[0200]
2.2.rpim
‑
mdtn离散系统方程
[0201]
在加权残差公式中，采用rpim构造的形函数n(x)作为检验函数w(x)，则式(16)表示为：
[0202][0203]
从而，可建立如下矩阵形式的离散系统方程
[0204]
[k
‑
k2m+k
b
]{p}＝f
ꢀꢀꢀ
(29)
[0205][0206][0207][0208][0209]
式中：{p}表示声压值的矢量，并且
[0210]
c
n
(x)＝[cosnθ sin nθ]
ꢀꢀꢀ
(34)
[0211]
以上这些矩阵是基于节点i和j的全局索引组装而成，计算这些系统矩阵需要一组独立于用于场变量插值的场节点的全局背景单元，这些背景单元对插值精度没有影响。注意，只有当场节点i和j同时属于至少一个积分点的支持域时，k
ij
≠0和m
ij
≠0，否则，k
ij
＝0
且m
ij
＝0。因此，全局系数矩阵k和m为稀疏矩阵。式(29)表明，mdtn边界条件对标准rpim的影响是包含了k
b
。与k
ij
和m
ij
相比，由于mdtn映射的非局部特征，对γ
r
上的任意节点i和j，k
ijb
不为零。研究表明，在大型问题上使用这种非局部边界条件产生的额外计算耗时可忽略不计。利用求解得到的人工边界处的声压值，通过式(7)可计算得到v
o
中任意场点处的声压值。
[0212]
3.数值验证
[0213]
本发明采用数值算例对所提方法进行验证，通过使用具有解析解的浸没式无限长圆柱数值实验获得了用于水下声辐射计算的影响域无量纲尺寸和形状参数的取值，并对该方法的精度、收敛性以及计算效率开展了研究。通过浸没式无限长圆柱扇形辐射模型、舵形结构和二维潜艇形结构的声辐射算例，进一步对所提方法进行了验证。
[0214]
3.1.浸没式无限长圆柱周向简谐声辐射算例
[0215]
如图4所示，研究了半径为a＝0.25m的无限长圆柱的声辐射，声介质为水，密度取ρ＝1000kg/m3，声速c＝1500m/s。以坐标原点为中心，构建r＝0.5m圆形人工边界。在圆柱的湿表面上(r＝a)指定周向dirichlet边界条件p(θ)＝cos(nθ)，对应的归一化解析解为：
[0216][0217]
在本发明中，分别采用ε
real
和ε
imag
来表示实部和虚部数值声压值相对误差的均方根误差，定义为：
[0218][0219][0220]
式中：p
iexact
和p
inume
分别表示在节点x
i
处的声压解析解和数值解。接下来，采用几种辐射模式对所提方法开展研究。
[0221]
3.1.1.影响域无量纲尺寸的影响
[0222]
通常，对一类基准问题，需进行数值实验来预先确定影响域无量纲尺寸α
s
的可用值。通过给出声压预报误差随α
s
的变化，研究了无量纲尺寸对所提方法性能的影响，使用平均节点间距为d＝0.023m的一组不规则分布节点表示有限计算域，采用496个二次元单元作为全局背景单元进行积分运算，如图5所示。声波波数满足ka＝π/2，以确保半径a恰好等于一个波长，构造形函数时在式(17)添加线性多项式。为了确保无网格全局弱式方法的积分精度，对于二维问题，积分点的数量n
q
应满足n
q
>2n
f
/3，其中n
f
为场节点的数量。由于n
f
＝1171，在全局背景单元的每个方向上使用5个gauss
‑
legendre积分点，以确保系统矩阵系数的计算精度。图6和图7分别给出了第二阶和第三阶辐射模式(n＝1,2)下人工边界处的ε
real
和ε
imag
随α
s
的变化曲线，其中a
c
＝0.4和q＝1.03(mq)，a
c
＝2(exp)，η＝3(tps)，n＝7(mdtn)。
[0223]
由图中曲线可知，声压预报误差随无量纲尺寸取值的变化而变化，无量纲尺寸取值对所提方法的精度具有显著的影响。为了获得精确的计算值，必须仔细确定α
s
的取值。当影响域尺寸过小(α
s
≤1.5)时，对应的rpim中的矩阵g为奇异矩阵，导致计算失败。因此，图中没有给出对应α
s
≤1.5计算结果。对于exp
‑
rbf，用所提方法获得的结果在2≤α
s
≤6范围是
可接受的；当α
s
<2或α
s
>6时，精度会降低。对于mq
‑
rbf和tps
‑
rbf，在α
s
≥2范围内，所提方法可以产生精确的结果。可以看出，随着无量纲尺寸α
s
的增加，并不能有效地提高计算精度。但是，较大的影响域将显著增加计算耗时，尤其是对于大规模问题。因此，需保持计算精度和计算耗时之间的平衡。结果还表明，使用mq
‑
tpf和tps
‑
rbf作为基函数具有比exp
‑
rbf更好的精度。为了不增加计算耗并获得较好的计算精度，在后续算例中，对exp
‑
rbf取α
s
＝3，对mq
‑
rbf和tps
‑
rbf取α
s
＝4.5。
[0224]
3.1.2.形状参数的影响
[0225]
为研究形状参数对所提方法的影响，给出了人工边界处的ε
real
和ε
imag
随形参数的变化，相应的结果如图8～图11所示。由图中曲线可知，计算精度受形状参数的影响很大，计算结果对形状参数取值敏感。因此，必须仔细选择这些rbf中包含的形状参数的值，以获得精确的计算结果。从图8可以看出，当采用exp
‑
rbf时，α
c
可取值范围为0≤α
c
≤5。图9中的曲线表明，形状参数q对计算结果的精度有显着影响，q＝j(j＝1,2,3,
…
)时，rpim中的g矩阵为奇异性矩阵，可能导致错误的结果。对于exp
‑
rbf，从图10可以看出，α
c
在0.07～2.5可以产生可接受的结果。当使用tps
‑
rbf时，从图11给出的结果可知，η应在1≤η≤9范围取值；注意，η取值接近2j(j＝1,2,3,
…
)会使矩阵g的条件数变大，从而导致精度变差。根据以上结果，在后续研究中，形状参数取值为：对mq
‑
rbf取α
c
＝0.5和q＝1.1，对exp
‑
rbf取α
c
＝1.9，对tps取η＝3.5。
[0226]
3.1.3.有限项数n的影响
[0227]
研究了mdtn映射中的有限项数n对所提方法的影响，图12给出了第五阶辐射模式下人工边界处的误差随有限项数n的变化，其中波数k满足ka＝2π。为了比较，在同一图中还给出了使用标准dtn边界条件获得的数值结果误差。从结果可以看出，与截断的dtn边界相比，使用mdtn边界条件可以提高计算精度，尤其是对于较小的n值。对于截断的dtn边界条件，通过设定n≥kr，可以确保适定问题，但不能总是保证所需要的精度；另外，通常不知道应预先采用多少项来达到期望的精度。与截断的dtn边界条件不同，mdtn边界条件对积分算子内核中包含的低阶模态是精确的，而局部算子包含了与高阶模态对应的近似非反射边界条件。因此，采用mdtn时，误差在n＝5时就快速收敛到最小值，在后续算例中，取n＝6。
[0228]
3.1.4.声场精度
[0229]
为了进一步检验所提方法的准确性，图13给出了三种辐射模式(n＝3,4,5)下声压幅值沿人工边界的分布情况，其中k满足ka＝2π。图14给出了第七阶辐射模式下的声压频响，其中观测点设定在r＝0.5且θ＝π/2。图13所示的结果清楚地表明，用所提出的方法获得的数值解与精确解吻合良好。从图14所示的结果可以看出，在整个频率范围内，数值解与精确解具有很好的一致性。此外，还可以观察到，mq
‑
tpf和tps
‑
rbf比exp
‑
rbf具有更好的精度。
[0230]
3.1.5.收敛性和计算效率研究
[0231]
为研究所提方法的收敛性，使用不同的平均节点间距(d＝0.011,0.020,0.032,0.043)进行仿真。对于第三阶辐射模式，图15给出了三种rbf对应的收敛曲线，其中k满足ka＝2π。为了比较，对应fem的收敛曲线绘制在同一图中。对于所有方案，采用相同的配置，使用二次有限元单元和相同数量规则分布的节点。图15中的曲线表明，所提出的方法的收敛速度比fem的收敛速度更快，在相同的配置下，所提方法可以产生更精确的结果。当使用mq
‑
ebf或exp
‑
rbf时，收敛不稳定，一个可能原因是无量纲尺寸与形状参数不匹配。到目前为止，还没有可用的方法从理论上确定无量纲尺寸和这些形状参数的取值，只能通过数值实验获得，需要更详细的研究来求解它们的取值。
[0232]
通过与有限元法的计算耗时进行对比来开展对所提方法的计算效率的研究，所有这些程序都在相同的软件和硬件(w
‑
2255 cpu 3.7ghz和ram 64gb)环境中执行。在该研究中，采用tps
‑
rbf。图16(a)给出了在ka＝2π时，cpu计算时间与平均节点间距的关系。结果表明，在相同配置下，所提方法需要更多的计算时间。其原因是：对所提方法，必须为积分点支持域内的场节点重新构造rpim形函数，而对fem能预先确定覆盖积分点单元的fem形函数；此外，与有限元方法相比，所提方法中需使用更多的节点信息来组装系统矩阵，从而所得系统矩阵的带宽通常大于fem。众所周知，成功的数值方法应以较低的计算耗时实现高精度，因此应将计算耗时与精度结合来公平地评估计算效率。图16(b)给出了cpu计算时间随声压实部误差的变化。从曲线可以看出，当ε
real
＝17.5％时，该方法的计算耗时约为fem的39％。因此，考虑到精度，该方法实际上比fem具有更好的计算效率。
[0233]
3.2.浸没无限长圆柱扇形声辐射算例
[0234]
为进一步验证所提方法，采用浸没无限长圆柱的扇形声辐射进行仿真，对于该算例，在圆柱边界上指定p＝1，其余圆柱边界处为p＝0，如图17所示，几何尺寸，相关参数和平均节点间距与前面的算例相同，相应的归一化解析解为：
[0235][0236]
图18和19分别为对应于k＝4π和8π，在人工边界处的声压分布。可以看出数值解与解析解具有很好的一致性。在高波数下，相应的声压值与解析解略有偏离，计算结果仍可接受。这些结果再次表明，与exp
‑
rbf相比，使用mq
‑
rbf和tps
‑
rbf可以获得更准确的结果。
[0237]
3.3.舵形结构声辐射算例
[0238]
本发明开展浸没于水中的刚性舵形结构声辐射数值试验，其中水密度ρ＝1000kg/m3，声速c＝1500m/s，该结构的形状和尺寸见图20。在结构的表面施加v
n
＝10
‑5m/s的neumann边界条件，以r＝4m的圆作为人工边界，有限计算域使用1193个双线性四边形单元离散(1277个节点)。为进行对比，在相同场节点和单元配置下，fem结果一并给出。由于该问题不存在解析解，采用细化网格(69,357个节点，68,677个四边形单元)对应的fem解作为参考解。
[0239]
在rpim插值中采用tps
‑
rbf作为基函数。图20、图21和图22分别给出了kl＝2π,4π和8π时人工边界处的声压分布。从这些图中可以看出，在三个频点处，rpim
‑
mdtn获得的声压分布与参考值吻合良好。对于较小的波数值，fem
‑
mdtn结果与参考解具有良好的一致性，但是随着波数的增加，fem的结果与参考值相差明显，尤其是在kl＝8π时。从对比结果可以发现，所提方法对波数不敏感，在水下声辐射计算中fem具有更高的计算精度。
[0240]
3.4.潜艇形状结构声辐射算例
[0241]
本发明研究二维潜艇结构的真实声辐射，该潜艇沉浸在密度为ρ＝1000kg/m3，声速c＝1500m/s的水中，其形状和尺寸如图23所示。对于潜艇，需要评估辐射噪声水平。而通常，敌人声纳距离被探测物体很远。在所提的rpim
‑
mdtn方法中，可以在获得人工边界上的声压数据之后使用公式(7)计算任意远场点处的声压。实际上，在预报潜艇噪声时，应考虑
与潜艇的尾部和前部船体相对应的的振动模态，这取决于分析频率。为简单起见，仅考虑位于潜艇尾部的发动机引起的船尾振动，辐射的噪音主要来自船尾。在此算例中，在尾部施加振动速度为v
n
＝10
‑5m/s的neumann边界条件，如图23所示。以原点为中心，r＝50m的圆作为人工边界。在人工边界上选取两个点a(θ＝0)和b(θ＝π/2)作为观测点，采用10,545个双线性四边形单元(10,841个节点)离散计算域。由于该结构声辐射问题不存解析解，使用细化的网格(179,005个节点，177,820个四边形单元)对应的fem解作为参考解。
[0242]
图23～图25给出了在f＝100hz,200hz和400hz时，声压沿人工边界的分布。在计算声压时，采用tps
‑
rbf。为了进行比较，在相应的图中还给出了了使用相同配置(10,841个节点和10,545个双线性四边形单元)的fem结果。由结果可知，在低频f＝100hz时，rpim
‑
mdtn和fem
‑
mdtn均可产生良好的结果；随着频率的增加，rpim
‑
mdtn结果依然非常接近参考解，而fem
‑
mdtn结果明显偏离了参考解，特别是在f＝400hz时。为进一步检验所提方法的精度，研究观察点处的10hz至400hz的声压频率响应，其中频率步长为1.0hz。对于观察点a和b，rpim
‑
mdtn和fem
‑
mdtn的计算结果以及参考解分别绘制于图27和图28。显然，在全频率范围内，所提出的方法比fem
‑
mdtn具有更高的精度。对于高频，rpim
‑
mdtn的解与参考解之间仍然具有良好的一致性。然而，fem
‑
mdtn在高频处不能获得可接受的精度。此算例再次表明，与传统的fem相比，所提出的方法对波数的敏感性要低得多。相比fem，所提出的方法能在实际的水下声辐射预报中表现得更好。
[0243]
应当注意，本发明的实施方式可以通过硬件、软件或者软件和硬件的结合来实现。硬件部分可以利用专用逻辑来实现；软件部分可以存储在存储器中，由适当的指令执行系统，例如微处理器或者专用设计硬件来执行。本领域的普通技术人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器控制代码中来实现，例如在诸如磁盘、cd或dvd
‑
rom的载体介质、诸如只读存储器(固件)的可编程的存储器或者诸如光学或电子信号载体的数据载体上提供了这样的代码。本发明的设备及其模块可以由诸如超大规模集成电路或门阵列、诸如逻辑芯片、晶体管等的半导体、或者诸如现场可编程门阵列、可编程逻辑设备等的可编程硬件设备的硬件电路实现，也可以用由各种类型的处理器执行的软件实现，也可以由上述硬件电路和软件的结合例如固件来实现。
[0244]
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。