基于水平集的二维半空间声屏障优化设计

张 放，姜富航，陈海波

（中国科学技术大学中国科学院材料力学行为与设计重点实验室，合肥 230027）

近年来随着工业的发展，噪声控制已经逐渐成为产品设计、城市规划等领域重点考虑的课题。在众多降噪手段中，利用隔声设备隔断声波的传播路径是一种简单有效且被广泛使用的方法，比如都市中常见的声屏障。而隔声体的隔声效果与其形状、分布及材料属性有关，同时针对不同的声场也不尽相同。为了达到预期的降噪效果，可通过声学结构优化来设计出满足要求的隔声体。

以声屏障为例，现有的研究包括参数优化[1]、形状优化[2]及表面吸声材料的拓扑优化[3]等。将声屏障简化成二维结构，在给定声源的情况下，分别优化各尺寸参数[1]、边界形状[2]以及表面吸声材料分布[3]。但以上方法都需要给定初始结构，且仅能对结构进行局部调整，而无法充分利用设计空间。若不限制障碍物的具体形状，则可利用拓扑优化的方法对结构本身进行设计。拓扑优化的最终目的是得到材料的最佳分布，而不受结构外形的约束，因此可充分利用设计区间，具有更高的优化自由度。

水平集方法是一种用于追踪界面演化的数值技术，可以方便地表征物体的拓扑结构改变[4]。与传统拓扑优化相比，其优点在于避免了以单元密度作为设计变量而导致的如棋盘格、中间密度等数值问题。当前阶段，水平集理论经过发展也产生了多个分支，如标准水平集、参数化水平集和反应扩散方程水平集等方法。尽管求解手段各异，水平集方法的核心思想都是用高一维度的水平基函数的0等高线表示结构边界，也称为零水平集，通过迭代更新水平集函数实现结构边界的追踪。

反应扩散方程水平集由Yamada 等[4]提出，根据虚拟界面能的思想，利用反应扩散方程实现水平集函数的更新，并被应用于最小柔度优化问题[5]中。这一方法避免了标准水平集方法中求解Hamilton-Jacobi 方程中遇到的诸多数值困难，但需要使用Fréchet 导数作为拓扑导数进行计算。在声学问题中，Matsumoto等[6]基于伴随变量法给出了拓扑导数的高效计算方法。Ⅰsakari 等[7]利用快速多极边界元法进行了三维声学结构拓扑优化。Sugihara 等[8]则对Robin边界条件下二维结构进行优化。最近，Gao等[9]针对边界积分形式的目标函数，给出了二维结构声学拓扑导数的具体推导，并用于声子晶体的优化研究中。以上研究均是在全空间自由声场中进行，而实际生活中地面对声波的散射往往不容忽略。在计算声屏障的声学性能时，通常可用半空间格林函数模拟地面对声波的反射效果[2]。

本文将针对半空间的声学问题，以边界元法对声屏障的降噪性能进行计算分析，并借助Yamada等[4]提出的反应扩散方程水平集方法，对声屏障截面进行优化设计。为了避免优化中产生孔洞而破坏结构的拓扑关系，仅取零水平集的主体部分作为其截面的实际形状，最终达成形状优化的目的。同时为了区别于Gao等[9]的研究，将以域积分作为目标函数的定义形式，并给出对应的伴随方程。该方法适用于任何域积分形式的目标函数，也可用于更大规模的半空间声学优化问题中。

1 二维半空间声学计算

理想流体介质中频域声场由Helmholtz 方程及边界条件控制：

以二维声学边界元法对环境中的声场进行求解，入射波作用下的边界积分方程如式（2）所示：

式中：c（x）根据点x在域内、光滑边界上和域外分别取1、1/2和0；p（x）为声压，q（x）为声通量，且q（y）=为格林函数，也称基本解；q*（x,y）为基本解的法向导数，且q*（x,y）=为点y处的外法向矢量；pi（x）为入射声压。

区别于全空间格林函数，对于二维半空间问题，格林函数包含镜面项：

式中：k为波数；r=|x-y|代表场点x与源点y之间的距离；r′=|x′-y|代表场点x关于地面的镜像点x′与源点y之间的距离；H0（1）为第一类0 阶Hankel函数。

数值计算中，首先将结构边界Γ划分成一系列的边界单元，再对式（2）进行离散，并计算出每个边界单元上的积分系数后组装成：

式中：p和q分别是边界上节点对应的声压和声通量组成的列向量；H和G是对应的系数矩阵；pi是入射声压在边界节点上值的列向量。

2 水平集方法

2.1 水平集函数

水平集函数φ（x,t）是设计域内点x坐标（x，y）的函数。在二维问题中，φ（x,t）本身对应第三个维度的坐标，也即高度z。以其零水平集作为结构的边界，以区分材料域与声场区域：

式中：Ω为声场域为材料域；Γ为两域交界面。

2.2 优化模型

定义目标函数如式（6）所示：

式中：f是声压的函数，用于建立域积分形式的统一格式。优化的目的是得到最优的材料区域配置，使得目标函数达到全局或局部的最小值。水平集方法通过改变水平集函数φ（x,t），更新结构边界，进而影响到观测区域的声压分布。因此若希望使目标函数逐步减小，则需要使水平集函数沿着其负梯度方向演化。Yamada等[4]提出的反应扩散方程如下所示：

式中：T（x,t）是目标函数的拓扑导数，表明当前位置的拓扑改变对目标函数的影响。K＞0 是正则化参数，负号则满足了优化的负梯度要求。第二项里，Δφ是曲率项，用于确保水平集函数的光滑性，系数τ用于人为指定迭代过程中水平集函数变化的几何复杂度。

将时间项φ（x,t）/∂t以有限差分格式表示，φ（x,t）以有限元法进行离散，最终可建立起水平集函数φ（x,t）的迭代关系式，详细推导可参考Otomori等[5]对结构柔度优化问题给出的样例程序，但须注意他们使用的数值模型里每个单元长度固定，即最终模型会随单元数增大而变大，并不适用于指定区域内的优化计算。

2.3 拓扑导数

拓扑导数T（x,t）定义为，声场在位置x处产生一个小孔后，目标函数的改变量与小孔面积（二维情况下）比值的极限，如式（8）所示，用于表示目标函数受当前位置声场拓扑变化的影响。

式中：表示增设小孔后的目标函数值；JΩ则表示小孔产生前的目标函数值；ε为小孔半径；最终结果T（x,t）为小孔圆心x处的拓扑导数值。当T（x,t）＞0 时，在x处铺设材料将会使目标函数增大，反之则减小，因此可以作为φ（x,t）随时间的迭代依据。

依据Matsumoto等[6]给出的拓扑导数推导方法，计算域积分形式的目标函数，拓扑导数最终表达式与Gao 等[9]全空间问题给出的表达式一致，如式（9）所示。

式中：p（x）为声压，由微分方程（1）控制；λ（x）为伴随变量，由微分方程（10）控制；λ,j（x）与p,j（x）分别表示λ（x）和p（x）的梯度分量；Re 表示取复数的实部。此时，由于目标函数选取为特定区域的声压积分，伴随场的控制方程中存在域内的非齐次项，如式（10）所示。

声场和伴随场的微分方程均通过边界元方法求解。在积分计算系数矩阵时，尽管式（1）和式（10）并不完全相同，但由它们推导得到的边界积分方程系数矩阵H和G仍是一致的，只是λ（x）计算中需要将式（2）中的入射波自由项替换成与f有关的域积分而后根据边界条件组装系数矩阵及自由项，就可求得p（x）和λ（x）的边界值。

域内点的计算由积分方程给出，如式（11）和式（12）所示。p,j（x）与λ,j（x）则可通过将其对坐标方向求导得到。

最后根据式（9），可直接得到点x处的拓扑导数值，并结合反应扩散方程对φ（x,t）进行更新。

在电气工程中，首先要做的是确定电气工程中电压等级的选择，以获得节能效果。例如，当前通用电气工程中的中压配电电压为10kV或者20kV，低压配电电压为220V或380V[6]。对于电压的选择，需要进行全面考量。电压选择与许多因素有关：一是国家电网在当地的供电设备的特点；二是电气工程中的电源电路数量；三是电力用户的计算电压容量；四是电气工程所在地区的发展水平。高压电网和低压电网在选择时也有不同：选择低压配电网的情况是用户的电力用于低功耗；选用高压配电网的情况是用户的计算能力，电气设备，单电源仪表电源的距离和计算负荷超过一定的范围内。在电压选择过程中，需要参考上述要求选择正确的电压以避免浪费。

2.4 优化流程

整个优化的流程如图1所示。首先给出初始水平集函数，其零水平集对应结构的初始构型。在半空间边界元计算时，需要对结构进行单元离散，用于积分系数的计算及系统方程的组装，并结合边界条件计算出所有边界值。根据域内点积分表达式计算观测点处的声压值之后，即可由式（9）得到设计区域内有限元节点处的拓扑导数值。进而依据式（7）所示的反应扩散方程更新水平集函数，对比目标函数的改变量，以判断是否收敛。

图1 水平集优化流程图

若并未收敛，则需要根据更新后的水平集函数重新提取其零水平集，并划分网格。值得注意的是，在结构的声学拓扑优化中，很容易产生多个独立悬浮的微小散射体，尽管确有一定的隔声效果，但难以实际应用，因此可以考虑仅保留面积最大的主体部分作为实际声屏障的截面。这一行为限制了结构拓扑关系的改变，最终表现为截面的形状优化。

3 数值算例

3.1 Γ形声屏障优化模型

以Γ 形声屏障的优化为例，y=0 m 作为地面，用位于（0,1）处的单位点声源模拟公路上直行的汽车噪声。则环境中的入射波场可考虑为：

式中：k为波数，r代表当前位置x与点声源之间的距离。

声屏障宽为0.2 m，高为3 m，左下角置于（5,0）处。观测点设置为x∈（6.1,16），y∈（0.1,3）之间均匀分布的100 个点，具体模型如图2 所示，计算频率为400 Hz，声速为340 m/s。

图2 声屏障截面优化模型

半空间问题中，声屏障底部与地面贴合，而地面作为镜面，视作刚性边界条件。在进行边界元分析时，通过选取半空间格林函数作为基本解，模拟声波在镜面上的反射过程，可以抵消声屏障在地面上的边界积分，从而只需要对地面之上的部分进行建模计算[2]。

设计域高为4 m，宽为2 m。在求解反应扩散方程时，给定设计域边界处φ（x,t）=-1 作为有限元法的边界条件。

目标函数J可表示为：

则在微分方程（10）中有：

3.2 拓扑导数验证

首先为展示式（9）所示的拓扑导数在半空间情况下的正确性，分别选择设计域内y=3.2 m 和4.0 m直线上的有限元节点，利用式（8）以有限差分方法进行验证。小孔半径设置为ε=1×10-4m，对比结果绘图如下：

由图3 可见，式（9）计算结果与差分方法得到的拓扑导数值完全一致，证明了基于伴随变量法计算得到的拓扑导数的正确性，其在半空间问题中也能够很好地反映出当前位置的拓扑变化对目标函数的影响。

图3 半空间拓扑导数差分验证

3.3 半空间水平集优化

将式（9）得到的拓扑导数作为水平集函数更新的梯度依据，可以使目标函数逐步减小。每次迭代中基于面积对零水平集进行的过滤处理，则能够保证结构拓扑关系不发生改变，从而只调整其外形。本小节将采用了过滤处理的形状优化结果与不采用过滤处理的拓扑优化结果进行对比，并利用商业软件COMSOL计算出优化前后的声场分布。

水平集函数的初始配置由离散后的有限元节点指定，在此将其指定为二值函数，位于声屏障内部赋值为1，否则为-1。同时，为了固定声屏障的位置，将0.5 m 作为声屏障的基座尺寸，即y＜0.5 m 的有限元节点在迭代中固定不变，取为初始值。

依据式（9）得到的拓扑导数，在反应扩散方程式（7）中对水平集函数进行迭代更新。参数K依据Otomori 等[5]给出的样例程序取为面积积分与灵敏度积分之比，τ取为1×10-5。优化步长dt初值取为0.06，并在目标函数增大较为明显时适当减小。形状优化和拓扑优化的目标函数的迭代收敛曲线如图4 所示，纵轴以对数坐标表示。

图4 目标函数迭代历史

观察迭代历史，两种处理下目标函数均能够保持稳定下降直至收敛。最终收敛值相近，都相比初值降低了两个多数量级，对应的结构对比如图5 所示。同时与Liu 等[2]对声屏障顶端的等几何优化将近80%的降噪幅度相比，优化效果更为明显。

图5 优化结果对比

观察图5，拓扑优化结果中出现了多个微小且独立的孔洞，而声屏障主体结构仅发生了几处局部的变形。形状优化结果中，结构拓扑关系并未改变，但整体形状与初始结构有了很大的区别。

在实际背景下，每一块孤立的二维散射体都代表一个以此为截面的无限长杆状结构。拓扑优化中出现的多个不规则孔洞对应着各自的杆状体，其与主体结构分离，悬浮于空中，将给生产加工带来极大的不便。而形状优化后的结构外形更为复杂，为得到相当的隔声性能使原本平整的结构外形变得起伏不定。

理论上拓扑优化由于设计自由度更高而应该能得到更好的优化结果，但若是希望只调整结构的外形，那么进行面积过滤的形状优化则是一种值得尝试的方法。这两种处理并没有优劣之分，只是按照实际需求进行取舍。最后利用商业软件COMSOL对优化前后结构周围的半空间声场进行计算，对应的声压级云图如图6所示。

图6 声压级云图对比

图6（a）是初始结构声场分布，图6（b）和图6（c）分别代表拓扑与形状优化后的声场分布。与初始结构的声场对比，观测区域云图由以40 dB 为主的金黄色区域变成以0～20 dB为主的蓝绿色区域，同时也出现了负声压级的部分深蓝色区域，整体声压级有了显著的降低。

4 结 语

对于二维半空间问题而言，通过伴随变量法计算的拓扑导数表达式与全空间一致，但需要利用半空间格林函数求解声压、伴随变量及其梯度。在提取零水平集作为结构边界时，可以取面积最大的部分作为目标结构，以限制其拓扑关系的改变。最终随反应扩散方程更新结构外形，同样能实现相当的优化效果。

基于反应扩散方程的水平集方法根据拓扑导数对水平集函数进行调整，不依赖于迭代历史，仅针对当前结构进行优化，稳定可靠。对声屏障截面进行优化设计时，降噪效果极为显著，能够为实际工程中的优化设计提供清晰具体的理论指导。