一种基于完美匹配层的无限声场仿真方法与流程

1.本发明涉及工业仿真技术领域，特别是涉及一种基于完美匹配层的无限声场仿真方法。


背景技术：

2.对于无限声场问题，如水下结构的声辐射、声散射等问题，以及无限域中声传播等问题，通常采用声学边界元，声学无限元，声学有限元等方法，以上数值方法均需要对计算域进行截断处理。
3.现有比较主流的截断处理方法为完美匹配层(pml)方法，它是在计算域的周围增加一个吸收层，模拟各向异性材料，吸收辐射到该区域的散射场，计算域边界不会产生虚假反射。但在实际应用中，由于pml必须在感兴趣区域的有限距离上被截断，因此它的外部边界产生了人工反射。此时，对问题进行离散化求解通常会导致误差变大，虽然通过增加pml厚度等方法可以降低计算误差，但是这种方式会大大增加计算成本，导致有时候并不可行。


技术实现要素：

4.本发明所要解决的技术问题是提供一种基于完美匹配层的无限声场仿真方法，能够有效地解决虚假反射，计算成本高，难以实现等问题。
5.本发明解决其技术问题所采用的技术方案是：提供一种基于完美匹配层的无限声场仿真方法，包括以下步骤：
6.(1)在不考虑流动情况下，构建无源声场的亥姆霍兹方程；
7.(2)在计算域周围赋予完美匹配层区域，并使用复数坐标形式表示增加完美匹配层区域后的亥姆霍兹方程；
8.(3)对所述计算域采用三角形非结构网格进行分割，对所述完美匹配层区域采用四边形结构网格进行分割；
9.(4)对分割后的所述计算域和完美匹配层区域采用有限元方法进行离散，得到弱形式的亥姆霍兹方程，并进行求解得到声压分布曲线。
10.所述步骤(1)中构建的无源声场的亥姆霍兹方程表示为：δp+k2p＝0，其中，δ为拉普拉斯算子，p为声场中的声压，k为波数。
11.所述步骤(2)中使用复数坐标形式表示增加完美匹配层区域后的亥姆霍兹方程具体为：其中，表示x方向映射系数，表示y方向映射系数，σ
x
表示完美匹配层区域在x方向的吸收函数，σy表示完美匹配层区域在y方向的吸收函数，v
x
表示网格点在x方向到边
界上的法向距离，vy表示网格点在y方向到边界上的法向距离，a表示完美匹配层区域在x方向的最小坐标值，b
*
表示完美匹配层区域在y方向的最大坐标值，b表示完美匹配层区域在y方向的最小坐标值。
12.所述弱形式的亥姆霍兹方程表示为：其中，ωf为计算域，ωa为完美匹配层区域，表示计算域内网格点的声压，表示有限元中型函数，表示pml区域内网格点的声压，γ表示计算域边界，g表示物理域的声波振幅对平面外法向量的求导，s表示边界长度离散单元。
13.所述完美匹配层区域在x方向的吸收函数所述完美匹配层区域在y方向的吸收函数c为声速。
14.有益效果
15.由于采用了上述的技术方案，本发明与现有技术相比，具有以下的优点和积极效果：本发明只需要选定边界，简单设置后可自动划分网格，外推出pml区域，构造的pml区域能够实现对声波的完全吸收，无虚假反射，且在离散化求解是不会令误差变大。
附图说明
16.图1是本发明实施方式的流程图；
17.图2是本发明实施方式中pml区域示意图；
18.图3是本发明实施方式中计算域和pml混合网格示意图；
19.图4是本发明实施方式的pml方法和特性阻抗方法计算结果对比图；
20.图5-图8是本发明实施方式的pml方法与声学商业软件comsol的结果对比图。
具体实施方式
21.下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本技术所附权利要求书所限定的范围。
22.本发明的实施方式涉及一种基于完美匹配层的无限声场仿真方法，该方法应用pml方法求解线性声学中的波传播谐散射问题，即标量亥姆霍兹方程的求解。如图1所示，具体包括以下步骤：
23.步骤1，在不考虑流动情况下，构建无源声场的亥姆霍兹方程。该亥姆霍兹方程表示为：δp+k2p＝0，其中，δ为拉普拉斯算子，p为声场中的声压，k为波数。
24.步骤2，在计算域周围赋予完美匹配层区域，并使用复数坐标形式表示增加完美匹配层区域后的亥姆霍兹方程。具体地说：
25.采用复数坐标形式，
26.对于二维声场问题，如图2所示，在计算域ωf周围赋予厚度为δ
pml
的pml区域ωa，则亥姆霍兹方程可表示为：其中，表示x方向映射系数，，表示y方向映射系数，σ
x
表示完美匹配层区域在x方向的吸收函数，σy表示完美匹配层区域在y方向的吸收函数，v
x
表示网格点在x方向到边界上的法向距离，vy表示网格点在y方向到边界上的法向距离，a
*
表示完美匹配层区域在x方向的最大坐标值，a表示完美匹配层区域在x方向的最小坐标值，b
*
表示完美匹配层区域在y方向的最大坐标值，b表示完美匹配层区域在y方向的最小坐标值。
27.综上可知，只要吸收函数为正值，任何传播波进入pml区域后都不会在边界出现反射，同时会迅速衰减。当pml区域足够厚或吸收函数值足够大时，传播波衰减非常快，以至于它们在pml区域外边界的反射非常弱；此外，反射波在通过pml区域返回计算域时进一步衰减，导致从pml区域外边界到计算域外边界的人工反射最小。
28.但是，当采用离散数值方法对问题进行求解时，这些优点没有一个被完全保留，因此必须对pml方法进行优化，以实现给定计算成本最佳可能解决方案。
29.步骤3，对所述计算域采用三角形非结构网格进行分割，对pml区域采用四边形结构网格进行分割。
30.本步骤中将计算域和pml区域采用不同的分割方式，形成了混合网格(见图3)。对于混合网格，计算域采用三角形非结构网格，pml区域采用四边形结构网格，计算精度取决于以下因素：1)单元函数的阶数；2)计算域内的单元大小；3)pml的单元层数；4)pml区域的厚度；5)吸收函数的形式。
31.步骤4，对分割后的所述计算域和完美匹配层区域采用有限元方法进行离散，得到弱形式的亥姆霍兹方程，并进行求解得到声压分布曲线。
32.本步骤中，对计算域和pml域采用有限元方法进行离散，亥姆霍兹方程方程可写为弱形式：
[0033][0034]
其中，表示计算域内网格点的声压，表示有限元中型函数，表示pml区域内网格点的声压，γ表示计算域边界，g表示物理域的声波振幅对平面外法向量的求导，s表示边界长度离散单元。
[0035]
通过对比计算发现，完美匹配层区域在x方向的吸收函数为在y方向的
吸收函数为时，计算精度最高，其中，c为声速。
[0036]
下面通过两个实施例进一步说明本发明。
[0037]
实施例1：二维管道
[0038]
本实施例中的二维管道长1m，宽0.05m，声速c＝343m/s。计算频率为343hz，采用四边形单元，网格尺寸0.01m，pml厚度为0.2m，划分15层。对比本实施方式的方法和特性阻抗方法的计算结果，取管道轴线处的声压分布曲线，如图4所示。从图4可以看出，两条声压曲线完全重合，出口处pml边界与特性阻抗边界条件计算结果一致，波长均为1m，与理论解一致，因此本实施方式的方法应用至2d管道声传播计算准确有效。
[0039]
实施例2：双声源干涉
[0040]
采用双声源干涉测试本实施方式的计算结果准确性，计算域为1m
×
1m二维平面，坐标原点在正方形中心，声源半径0.05m，声源距离0.1m，pml厚度设置0.1m，层数为10层。计算结果与成熟声学有限元软件comsol进行对比，图5为3000hz的实部声压云图对比，图6为3000hz的虚部声压云图对比，图7为x＝0处截线的实部声压分布曲线(虚线部分为本实施方式，实线部分为声学商业软件comsol)，图8为x＝0处截线的虚部声压分布曲线(虚线部分为本实施方式，实线部分为声学商业软件comsol)。由图5～图8可知，本实施方式的计算结果与声学商业软件comsol计算结果一致。