实现水下目标声散射精确建模是检验反潜声呐检测算法的重要环节。针对水下目标进行有限元数值仿真可以精确模拟水下目标声散射现象，从而有效仿真水下声场。在有限元计算中，大约60%的工作在于前期处理。其中，有限元网络划分是非常关键的，网络剖分的质量将直接影响最终的计算结果准确性。目前针对二维平面剖分，Delaunary三角剖分方法由于具有较好的理论基础和数学属性，已经成为有限元分析软件前处理的广泛选择。Delaunary剖分算法的主要优势体现在空圆特性以及最大化最小角理论，其避免了狭长三角形单元的产生[1]。可以极大地避免低质量单元格的产生。但是在实际应用中仍有很多问题，比如在复杂狭长区域的划分不均，域外三角单元问题等。近年来，针对Delaunary剖分，有很多学者提出了优化方案，比如L.A.Piegl提出的边界外法向法，凸分解法等，在解决域外三角单元方面都有各自的优缺点[2]。针对这种情况，文献[1]提出了一种算法优化方案，较好地解决了域外三角单元和划分不均的问题。但是在处理实际问题中，特别是舰艇水声问题下的有限元仿真建模，我们往往不要求求解域中的单元均匀划分，在流固耦合作用区域，我们希望加密网络得到较为精确的解，在水体等大范围环境参数变化不大的区域，我们希望减少单元格来提高运算速度。为了解决这个问题，本文中利用权重赋值建立舰艇舱室结构不同重点区域，在流固耦合界面、固体力学传递界面和声波传递界面建立不同的加密区域，利用控制单元节点插入次序在不同区域建立均等三角单元网络，提高网络整体质量[3]。
1 基本原理
1.1 经典Delaunary三角算法
先对Delaunary算法进行简单介绍。针对其方法的实现途径有很多，其中B-W算法因操作简单、易于实现等优点被普遍应用[4]。其核心思想是循环利用插入点与已知点性质来生成细化三角形[5]。其具体操作步骤如下：
步骤1 根据输入模型进行识别，生成初始离散点；
步骤2 在边界点的基础上构造包围划分区域的初始三角单元；
步骤3 根据B-W算法依次构建已知离散点与初始三角单元相互作用的小三角单元；
步骤4 删除包含有初始三角单元的三角形，得到已知点的三角剖分；
步骤5 根据优化准则进行局部优化；
步骤6 循环执行步骤3。其中B-W算法步骤如图1所示。
图1 B-W剖分步骤示意图
1.2 舰艇舱室结构的权重比例因子
在水声仿真主动探测领域，舰艇复杂的舱室结构特征会给回波信号带来不可忽略的影响，相互链接的肋状支撑结构从一定程度上使舰艇舱室的本征振动频率发生改变，从而出现意料之外的共振峰[6]，精确模拟舱室结构特征对于声呐回波信号仿真及目标识别具有重要意义。本文以潜艇结构为例，其结构模型也可以对舰船模型构建产生良好的借鉴。某型双壳潜艇构建内部舱室结构模型如图2所示，结构细化如图3所示。
图2 某型潜艇舱室结构模型示意图
图3 结构细化示意图
由于部分区域涉及流固耦合以及声固相互作用，需要进行网格加密来提高计算精度，[7]部分区域因为单元数值变化梯度低，可以采用稀疏网格单元增加计算速度，因此需要对不同网络作用区域进行权重调整。首先，以重要程度划分，对水体和大型舱室空气作用权重赋初始值[8]，其余结构依次是肋板接触区域、单一流固接触界面、声固耦合作用再辐射区域。通过全加密剖分下的有限元应力分析，可以大致指定权重比例。采用模型整体最下网络剖分下，潜艇舱室结构横切面模型对1 000 Hz平面波水平激励下的结构应力分析如图4所示。依据此可以大体指定舱室内部的权重值配比。
针对应力作用区域，划分权重值，如表1所示。
图4 应力分析示意图
表1 权重划分
结构权重分级水体1舱室空气区域2外壳域3内壳域4舱室连接域5
1.3 区域优化算法
由于Delaunary剖分算法在局部区域进行细化插值时，插值点选取不同而导致低质量网格出现，如图5所示。
图5 随机点选取示意图
当△ABC和△ABE为待细分三角单元时，传统算法可能会在△ABC中插入点，这样会造成图6中的情况，产生了△ADC与△BDE两种单元，这是有限元分析中不想见到的。
图6 传统算法生成网络示意图
在相同类型区域计算时，希望尽可能用一组等边三角形去拟合网络。[9]因此，在相同权重区域进行剖分时，可以在待细分单元中选择可以生成最佳三角单元的区域插入新节点，并实时检验有效性。其算法实施思路总结如下：
步骤1 识别模型边界，并进行均匀离散，生成初始的Delaunary三角单元集合p1。
步骤2 判断域外三角形。计算所有边界三角形的重心D，设其三角形的边界边分别为ci和cj，按照右手定则判断域外三角形，若方向为正向，则删除该单元。
步骤3 针对权重比例，设置不同的网格加密区域并设置相关期望网格尺寸{αi}，其中i表示权重等级，αi表示不同等级下的期望尺寸。在权重区域判断上，可以重点区域的边界点所构造的初始三角形为区域细化目标，针对不同区域设置不同集合，下面的细化操作均在不同区域的集合内进行操作。
步骤4 首先针对最低等级进行全模型网格细化。通过for循环，对所有初始网格进行尺寸对比。设目标三角形边长最大值为M，对满足M>α1的所有三角单元进行分层。将(α1,M)区间划分为小型区间，其过程可以表示为：
(1)
式(1)中：Z表示动态变量；I是随着每次三角形质量检查时的Z的递变量；L表示对单元质量进行检查的次数。
步骤5 将满足条件M>Z=(Z-I·i)的三角形单元确定为有限插入单元，并添加至统一的集合p2中。
步骤6 当一次循环中有多个三角形满足步骤4中条件时，需要确定插入节点的顺序。这里利用边长排序进行区分。将这些三角形按最长边降序排列，如果出现最长边相等，就以边长之和大的靠前排序；若三边之和也相等，就比较其最长边与最短边的比值，较大的靠前。
步骤7 将新节点插入排序后的三角形单元的形心处，构造新的三角单元并删除原有单元。之后利用Delaunary原理进行优化并更新集合p1和p2。
步骤8 当满足条件的单元节点插入结束后，更新i的值，进行下一阶段的for循环，在下个阶段的循环中，Z的值也随之变换。当i>L时，结束这一区域的加密。跳转至下一区域重复步骤4之后的计算处理。
还是以图4中三角形为例，在细分算法优化后，点会先在△ABE中产生，从而得到较高质量单元，相对分布较为均匀，如图7所示。
图7 均分后网络
2 算例演示
利用Matlab对Delaunary剖分进行算法优化。先利用Comsol对潜艇初始横截面模型进行构建，保存成matlab可以直接读取的dxf文件。进行初始网格划分，如图8所示。
图8 Matlab初始网格
随后根据设置的阈值进行细化插值，其改善前后舱室细化结构如图9、图10所示。
图9 传统局部细化
图10 均分局部细化
由于模拟实际尺度计算，仅给出局部来清晰展示剖分效果。由图9、图10可以发现，优化算法对于由Delaunary剖分带来的网格不均匀现象起到了明显的改善效果。
3 算法检验
3.1 网格质量
针对三角形网络质量，仿照文献[1]，引入网络平均质量系数qm进行判定，即：
(2)
(3)
式(3)中：ri为三角形内切圆半径；Ri为三角形外接圆半径，当其趋于1时，为最优情况。对比结果如表2所示。
表2 检验对比
算法名称网格数量/个平均质量系数优化算法851670.9373传统算法846720.7632
结果表明，WOA优化算法可在整体上提高了Delaunary剖分的网格质量[10]。绘制网络质量直方图，如图11所示。
图11 传统算法(上)与优化算法(下)网络质量统计直方图
由图11可知，WOA优化算法相对传统算法无论是总体质量还是质量分布上均有改善，有效满足了有限元仿真计算的需求。同时发现，优化算法较传统算法在网格数量上有所增加，原因可能是：为了保证计算精度，在部分区域设置的较低期望值使得优化算法在部分区域产生了比传统算法更多的细化网络所致。下面讨论两者在计算精度以及计算时间上的差别。
3.2 计算精度及效率
为了检验优化算法剖分后的计算精度，以权重最高的区域期望值为最终期望值，进行全模型极细化剖分计算，以追求最大限度的准确性，并与利用优化算法后针对同一模型的相同位置处的散射声压级的仿真计算结果作图，如图12所示。
图12 仿真计算结果
由图12可以发现，传统有限元算法在常规密度下的计算结果与极细化的剖分结果在低频上相似度较高，但是在高频处有所偏离。原因是没有考虑不同结构带来的计算影响，部分重点计算区域没能细化处理，使得结果发生偏离。而改进后的算法在计算结果上与极细化剖分后的的算法在全频段保持一致，其网格数以及计算时间如表3所示。
表3 计算效率统计
算法名称网格数量/个计算时间/h极细化78652317.5传统算法753921.2权重优化846721.5
由表1可以发现，在达到极细化计算精度的情况下，通过建立相应的权重等级优化可以极大减少构建网格数量，降低仿真计算时间。另外，优化算法虽然在网格数量上较传统算法有一定的劣势，但是在计算时间上相差不大，综合计算精度考虑，优化算法在计算复杂结构体方面是可取的。
4 结论
本文提出了针对舰艇复杂舱室结构水声有限元仿真的WOA算法优化，也可以拓展至其他领域的复杂结构计算。针对实际模型仿真涉及流固耦合以及声固相互作用，需要进行网格加密提高计算精度的需求设计权重配比，再对相同权重区域进行有效优化剖分，WOA算法可以在保证计算精度的情况下降低畸形网络生成概率，提高网络质量和计算速度。本文中的优化算法仅仅涉及二维平面计算，如何将其拓展到三维区域进行有效构建是下一步工作的重点。
[1] 张晶飞,李射,崔向阳.Delaunay三角剖分的最优化网格节点生成算法研究[J].电子设计工程,2019,27(09):10-16.
[2] 王群,李爱平,马淑梅.局部网格狭长三角形的品质改善及实现[J].同济大学学报(自然科学版),2004(11):1508-1511.
[3] 台淼,王佩贤,姜东.利用Delaunary三角网优化CORS网构的方法研究[J].测绘科学,2011,36(05):134-135+124.
[4] 徐辉,张树有.基于弦内角映射参数的多边形三角剖分快速优化[J].中国图象图形学报,2004(04):63-68.
[5] 赵晓东,晏小宝,沈永明,王亮.基于GIS的高质量约束Delaunay三角网格剖分[J].地理与地理信息科学,2010,26(05):24-28.
[6] 李克清.基于Delaunary三角剖分的目标穿越算法[J].常熟理工学院学报,2009,23(02):98-100.
[7] 丁永祥,夏巨谌,王英,肖景容.任意多边形的Delaunay三角剖分[J].计算机学报,1994(04):270-275.
[8] 林晓穗,李桂清,尹梦晓,等.优化变形能的平面多边形同构剖分[J].计算机辅助设计与图形学学报,2014,26(12):2115-2125.
[9] 杜贻群,张恺.水声换能器发射性能的有限元分析[J].四川兵工学报,2010,31(12):8-10.
[10] 韩晶,黄建国,苏蔿.一种基于声线理论计算非相干平均声场的方法[J].系统仿真学报,2007(16):3653-3656.