0 引　言

为提升潜艇声隐身性能，一般需在均匀粘弹性覆盖层材料中嵌入空腔、金属硬核、局域共振单元等散射体^{[1 − 5]}，从而构成周期结构，以实现低频禁带效应。为便于设计和优化声学覆盖层降噪性能，如何准确高效的获取声学覆盖层的声学特性（如传递矩阵、吸声系数等）是一项重要挑战。

覆盖层声学特性的获取主要包括试验法、解析法和有限元法3类，对小样材料，一般可在声波导管中直接测量吸声系数、传递阻抗等参数^[6]，这种测量的前提是假设波导管中只有平面波传播。解析法的优点是物理机理清晰且计算效率高，但只适用于较为规整的周期单元，主要是运用弹性波传播理论与散射体壁面处的边界条件建立方程通过传递矩阵求解声学系数。Liu等^[7]结合Bloch定理与反演技术对二维圆柱空腔覆盖层等效均匀，结果表明空腔剪切波散射、覆盖层与背衬之间的耦合效应是影响等效精度的重要因素。Ye等^[8]将变截面空腔分层划分为多层圆柱短管，采用传递矩阵研究其吸声性能，但这种方法将覆盖层近似为流体，仅考虑纵波而忽视了横波的影响。有限元可以求解任意形状周期单元。罗英勤等^[9]建立了轴对称有限元模型，分析了局域共振结构各项参数对吸声系数的影响规律并用遗传算法改进吸声结构拓宽吸声频带。Jia等^[10]以功能梯度材料为基体通过局域共振结构增强声耗散揭示了吸声机理。数值法虽可计算任意形状周期单元，但在进行多组参数分析和性能优化时求解耗时较长。

基于以上分析，本文提出一种数值加解析的联合算法，以声阻抗管三水听器双负载测试法为基础，建立有限元模型，获取含任意周期单元声学覆盖层的传递矩阵，再通过多层均匀介质声传播解析理论获取吸声系数，为高效计算覆盖层吸声性能以及降噪性能优化提供理论与方法支撑。

1 基本原理

首先参照声管中测量小样材料声学参数的三水听器双负载法^[11]，建立声学覆盖层传递阻抗的有限元计算方法。在此基础上，联合分层介质声传播解析理论获取材料的吸声系数。如图1所示，声管最左端为平面波声源，在待测样品右侧有一个可移动的刚性平面，整个系统需布放3个接收水听器。

1.1 声管中传递矩阵测试基本原理

以ASTM E2611-19标准为参考^[11]，要得到传递阻抗矩阵中的4个元素T₁₁、T₂₂、T₂₁、T₂₂，需构造式 （1）和式（2），其中a，b代表声管右侧2种不同的声场条件（可通过调整距离D_i来实现），p₀、u₀分别代表样品前表面x=0处的声压与振速，p_d、u_d则代表样品后表面x=d处的声压与振速。

$ \left[\begin{array}{l} p_{0} \\ u_{0} \end{array}\right]_{x=d}=\left[\begin{array}{cc} T_{11} & T_{12} \\ T_{21} & T_{22} \end{array}\right]\left[\begin{array}{l} p_{\Delta} \\ u_{d} \end{array}\right]_{x=s}，$

(1)

$ \left[\begin{array}{l} p_{0} \\ u_{0} \end{array}\right]_{x=6}=\left[\begin{array}{cc} T_{11} & T_{12} \\ T_{21} & T_{22} \end{array}\right]\left[\begin{array}{l} p_{d} \\ u_{d} \end{array}\right]_{x=6}。$

(2)

可得到样品传递矩阵的表达式为：

$\begin{split} &T = \left[ {\begin{array}{*{20}{c}} {{T_{11}}}&{{T_{12}}}\\ {{T_{13}}}&{{T_{14}}} \end{array}} \right] = \\ &G\left[ {\begin{array}{*{20}{c}} {{p_a}(0){u_b}(d) - {p_b}(0){u_a}(d)} & {{p_b}(0){p_a}(d) - {p_a}(0){u_b}(d)}\\ {{u_a}(0){u_b}(d) - {u_b}(0){u_a}(d)} & {{p_a}(0){u_b}(d) - {p_b}(d){u_a}(0)} \end{array}} \right]。\end{split} $

(3)

式中：$ G=1 /\left(p_{a}(d) u_{b}(d)-p_{b}(d) u_{a}(d)\right) $，下标a、b分别代表D_i取2种不同的距离。因此要得到传递矩阵需知道样品前后表面的声压与振速。式（3）主要针对周期单元沿厚度方向非对称分布的形式，如果周期单元沿厚度方向对称分布，式（3）可简化为：

$\begin{split} & T_{a m}=\\ & G\left[ \begin{array}{cc} p_{a}(d) u_{a}(d)+p_{a}(0) u_{a}(0) & p_{a}(0)^{2}-p_{a}(d)^{2} \\ u_{a}(0)^{2}-u_{a}(d)^{2} & p_{a}(d) u_{a}(d)+p_{a}(\mathrm{C}) u_{a}(0) \end{array} \right]。\end{split} $

(4)

1.2 有限元计算流程

首先，在有限元模块中建立如图2所示的模型，计算覆盖层样品前后端面处的声压和振速，利用式（3）即可得到传递矩阵。因为在覆盖层样品端面近场可能存在高阶非均匀平面波干扰，所以不能直接采用有限元法计算样品前后表面的声压与振速，而是先计算距离样品表面一定距离处水听器1、水听器2、水听器3的声压值，然后利用平面波声传播规律进行计算。

图2给出了利用有限元法计算待测样品传递矩阵的示意图，待测样品左侧$ x < 0 $任意位置处的声压与振速表达式分别为：

$ p_{i}(x)=e^{-j k x}+R_{l}(0) e^{k x}, u_{i}(x)=\frac{e^{-j k x}-R_{i}(0) e^{j k x}}{z_{f}}。$

(5)

待测样品右侧$ x>d $任意位置处的声压与振速分别为：

$\begin{aligned} & p_{i}(x)=p\left(L_{1}\right) \cos \left[k\left(x-L_{1}\right)\right], u_{i}(x)=\\ & -\frac{j p\left(L_{1}\right)}{z_{f}} \sin \left[k\left(x-L_{1}\right)\right]。\end{aligned}$

(6)

式中：$ R_{i}(0) $为覆盖层表面$ (x=0) $处的反射系数，$ p\left(L_{i}\right) $为刚性背称处声压，$ z_{f} $为流体介质的特性阻抗。若能确定$ R_{i}(0) $、$ p\left(L_{i}\right) $，则能确定波导管中任意位置处的声压与振速。

反射系数$ R_{i}(0) $可通过水听器1和水听器2之间的传递函数$ H_{12}\left(D_{i}\right) $得到：

$ R_{l}(0)=\frac{H_{12}\left(D_{i}\right) e^{jks}-1}{1-H_{12}\left(D_{i}\right) e^{-jks}} e^{2 jks}。$

(7)

式中：S代表水听器1和2之间的距离，$ l $为水听器2到样品左侧端面的距离，$ D_{i} $为覆盖层右侧端面到刚性壁面的距离。

联合（5）～式（7）式可得为覆盖层前、后表面的声压振速表达式分别为：

$ p_{i}(0)=1+R_{i}(0) \text{，} u_{i}(0)=\left(1-R_{i}(0)\right) / z_f ，$

(8)

$ p_{i}(d)=p\left(L_{i}\right) \cos \left(k D_{i}\right) \text{，} u_{i}(d)={i} p\left(L_{i}\right) \sin \left(k D_{i}\right) / z_f 。$

(9)

将式（8）、式（9）带入式（3），即可得到传递矩阵$ {\boldsymbol T} $。至此，将原来含复杂结构的声学覆盖层单元转化为均匀层处理。

1.3 传递矩阵法计算吸声系数

通过有限元模块计算得到式（3）中的传递矩阵元素后，利用分层介质波传递理论可得到声学覆盖层的各类声学参数。

反射系数：

$ R=\dfrac{ {\boldsymbol T}_{11}+\dfrac{ {\boldsymbol T}_{12}}{z_{f}}-z_{f} {\boldsymbol T}_{21}- {\boldsymbol T}_{22}}{ {\boldsymbol T}_{11}+\dfrac{ {\boldsymbol T}_{12}}{z_{f}}+z_{f} {\boldsymbol T}_{21}+ {\boldsymbol T}_{22}}，$

(10)

透射系数：

$ t=2\left[ {\boldsymbol T}_{11}+\frac{ {\boldsymbol T}_{12}}{z_{f}}+z_{f} {\boldsymbol T}_{21}+ {\boldsymbol T}_{22}\right]^{-1}，$

(11)

吸声系数：

$ \alpha=1-|R|^{2}-|t|^{2}。$

(12)

对于实际中的问题，覆盖层是与背称底板一起作用产生吸声或者隔声性能的，因此要覆盖层与底板的整体传递矩阵，整体矩阵T_e=TT_s，其中，T_s为底板的传递矩阵，具体形式为：

$ {\boldsymbol T}_{s}=\left[\begin{array}{ll} \cos \left(k_{x} h_{s}\right) & \bar{z}_{x} \sin \left(k_{y} h_{s}\right) \\ \dfrac{i}{Z_{x}} \sin \left(k_{s} h_{s}\right) & \cos \left(k_{s} h_{s}\right) \end{array}\right]。$

(13)

式中：$ k_{x} $、$ Z_{x} $分别为底板弹性波数、力阻抗；$ h_{s} $为底板厚度。

图3为该算法的流程。

2 算法验证及性能分析

为了验证本文所提算法的有效性，针对如下4种常见类型的声学覆盖层进行分析：1）含周期空腔的声学覆盖层；2）含周期硬核的声学覆盖层；3）含周期局域共振单元（硬核外包裹柔性材料）覆盖层；4）多层组合式吸声单元。

首先利用多物理场耦合有限元软件Comsol建立有限元模型（见图4）并计算声波垂直入射时的吸声系数（下文统称“FEM”），将计算结果与传递矩阵法（下文统称“本文算法”）计算的吸声系数对比。

计算中所采用的材料参数如表1所示，$ E $为表杨氏模量，v为表泊松比，$ \eta $为表损耗因子，c为表声速。

2.1 含周期空腔声学覆盖层

不失一般性，以常用的圆柱、圆锥和圆台空腔声学覆盖层为例进行验证。基体橡胶厚度为50 mm，晶格常数为30 mm，空腔高度为30 mm，钢背衬厚度为10 mm。圆柱空腔为沿厚度方向对称分布形式，可利用式（4）计算，圆锥空腔和圆台空腔则为沿厚度渐变的非对称形式，必须采用式（4）进行计算。当声波垂直入射时，利用非均匀层等效传递矩阵法计算含不同类型空腔的声学覆盖层吸声系数，并与直接有限元计算结果对比。

图5所示为底边直径为8 mm的圆柱、圆锥和圆台空腔在平面波垂直入射情况下的吸声系数（实线为FEM结果，虚线为本文算法结果），可发现不同腔型下2种方法计算结果基本趋于一致，在100～3000 Hz频率范围内的吸声系数曲线完全吻合，但是在 3000～6000 Hz频率范围内均产生一定偏差。这跟覆盖层在不同频率吸声机理不同有关，在第一阶吸声峰值对应的频率附近，主要能量损耗机理为背称与覆盖层形成的整体“质量-弹簧”共振有关，本文算法的前提是声管内的声场处处为平面波，“质量-弹簧”共振不会对这种前提产生影响。当频率大于第一阶吸声峰值频率后，覆盖层的声吸收逐渐以腔体形变导致的波形转换有关，这时候覆盖层端面会产生非同相位移，与平面波传播假设冲突。

以圆锥空腔为例（曲线），图5中 4000 Hz频率附近2种算法结果产生明显差异，从图5子图可看到，此时圆锥空腔产生了明显的挤压形变，入射的压缩波产生波形转化变成了剪切波，这与本文算法前提不相符，从而产生计算误差。

2.2 含周期硬核单元吸声系数验证

本节以含周期硬核单元覆盖层为研究对象计算吸声系数，基体材料参数参考表1，铁硬核高度均为30 mm。

图6中3种不同形状硬核覆盖层，采用有FEM直接计算的结果与传递矩阵计算结果在绝大部分频率范围都有较好一致性，仅在 3500～4500 Hz频率范围内有较小误差，误差值要比空腔覆盖层工况下的小。

含周期硬核覆盖层吸声机理主要包括2个方面，一是硬核在基体中的摆动共振（偶极子共振），另一个是多重散射^[12]。这两类模式均不会对覆盖层端面位移同相性产生大的影响，因此声管内声传播更接近平面波形式，计算误差要小于周期空腔结构。

2.3 含周期局域共振单元吸声系数验证

局域共振型声学覆盖层，由于能以“小尺寸控制大波长”，被广泛应用于低频吸隔声中。图7为典型局域共振单元声学覆盖层（金属硬核外包裹柔性层嵌入基体中），铁芯体半径和高度分别为8 mm、20 mm，包覆层半径和高度分别为12 mm、30 mm，其余结构参数同上文。

图8中，本文算法与有限元法在大部分频点处计算结果趋于一致，仅在3000～4500 Hz范围有一定差异，由文献分析可知^[13]，局域共振型声学覆盖层在第一阶吸声峰值附近的能耗主要取决于金属内核的共振，此时不会对覆盖层端面运动产生明显效果，因此2种算法有较好的一致性。随着分析频率增加，以3500 Hz为例，从图8子图可看出，此时的能量损耗主要集中在柔性包覆层上，通过包覆层的剪切形变实现吸声，这也导致覆盖层端面会产生非同相运动，使声管内不满足平面波传播特点。

2.4 多层组合式声学覆盖层

近些年，为了获得更低频的吸声能力，拓宽声学覆盖层的吸声带宽，研究人员提出了多层组合式周期单元，这种组合方式也使能量损耗机理更复杂。故针对3类常见组合展开分析：空腔-空腔、硬核-空腔、局域共振单元-空腔。

图9对应的单元基体厚度为30 mm、空腔与硬核高度均为20 mm，钢背称厚度为10 mm。由图可知，对于空腔-空腔组合形式的结果对比，与单层空腔（图5）的类似，这里不再展开分析。而硬核-空腔组合则与单层的情况有一些区别，从图9可看出 4000 Hz时，2种方法计算结果存在差异，由子图可知，4000 Hz时硬核有明显的大幅值运动，影响了覆盖层端面的运动情况。另外双层单元之间运动耦合也会对最终结果产生一定影响。

图10为局域共振单元与空腔的组合式声学覆盖层，2种算法对比结果如图11所示，在大部分频率范围，2种算法有较好一致性，仅在第二阶吸声峰值处有少许偏差，原因可能是局域共振单元柔性包覆层的剪切形变，以及双层单元的运动耦合作用导致的。

3 算法适用性分析 3.1 填充率的影响

当覆盖层内部散射体填充率增加时，会影响等效传递矩阵的计算结果。以圆柱空腔为研究对象，基体厚度为30 mm，空腔高度为20 mm，钢板厚度为10 mm，固定其余参数不变，讨论填充率的影响。

图12所示为3种不同直径的圆柱空腔单元吸声系数计算结果，可发现随着圆柱空腔直径增加（也就是填充率增加），2种算法计算结果误差越来越大，尤其是在第二阶吸声峰值附近。很显然大直径空腔更容易产生大幅挤压形变，导致本文算法不再适用。

3.2 入射角度的影响

在第2节的分析中，所有计算均假设声波垂直入射，本节将探讨斜入射下算法的适用性。图13分别为采用直接FEM和本文方法计算60°斜入射角下，不同周期单元覆盖层的吸声系数。

4种工况下，斜入射的吸声系数，分别为圆柱空腔、圆柱硬核、圆锥空腔、圆锥硬核单元。2种算法对比时，除了圆柱空腔工况误差较小外，其余工况均表现出较大误差，且空腔单元的效果要优于硬核单元。斜入射波更容易引起周期单元的剪切形变，从而影响覆盖层端面的位移同相性。另外，对于硬核单元，采用本文算法时会漏掉一个吸声峰值，这可能是因为覆盖层内的弹性波遇到硬核单元会产生多重散射，这时候多重散射会产生耦合共振，而本文算法的内在原理忽略了这一点。

4 结　语

为提高声学覆盖层吸声系数的计算分析效率，提出一种数值加解析的联合算法，以声管测试技术为基础，建立有限元数值模型，计算含任意复杂周期单元声学覆盖层的均匀化传递矩阵，结合多层均匀材料的解析理论求解吸声系数，得到如下结论：

1）本文算法能较好处理大部分典型声学覆盖层的声吸收计算，且有利于提高参数分析及性能优化效率；

2）覆盖层周期单元发生共振或是产生波形转换时，会影响本文算法的准确性，其原因大尺度形变会改变声管中平面波传播模式；

3）当覆盖层周期单元填充率较高，或声波斜入射时，采用本文算法计算误差较大，如何克服此类问题是后续研究的重点。