基于有限元法的跨介质声传输数值模拟与分析

Numerical Simulation and Analysis of Cross-Media Acoustic Transmission Based on the Finite Element Method

摘要: 在由空气、薄冰层与水体构成的跨介质体系中，由于不同介质之间存在显著的声阻抗差异，声波在界面处易发生强烈反射，导致声能透射效率显著降低，从而制约了薄冰环境及实验室理想条件下跨介质声学测试与感知研究的发展。针对这一问题，本文从数学建模与数值分析的角度出发，建立了一种适用于毫米级薄冰条件的跨介质声传输等效介质模型，并在此基础上构建相应的有限元数值计算框架。通过对单元结构进行参数化建模，将其等效为亚波长介质层，从而将复杂的多介质界面声学问题转化为可求解的边值问题。基于有限元方法对该模型进行数值离散与求解，系统分析了单元几何参数及阵列结构对声能传输特性的影响。数值仿真结果表明，在单元共振频率附近，体系中的声能透射效率显著提升，且有限元计算结果与等效介质模型预测具有良好一致性。研究结果验证了所建立的数学模型及有限元数值方法在毫米级薄冰及实验室理想环境条件下分析跨介质低频声传输问题的有效性，为相关声学结构的设计与优化提供了一种可行的数值分析工具。

Abstract: In cross-media systems composed of air, thin ice layers, and water, the pronounced acoustic impedance mismatch among different media leads to strong wave reflection at the interfaces, resulting in a significant reduction in acoustic energy transmission efficiency. This limitation restricts the development of cross-media acoustic testing and sensing research under thin-ice and laboratory ideal conditions. To address this issue, this study approaches the problem from the perspectives of mathematical modeling and numerical analysis. An equivalent medium model applicable to millimeter-scale thin-ice conditions is established to describe cross-media acoustic transmission, and a corresponding finite element numerical framework is developed on this basis. By parameterizing the unit-cell structure and treating it as a subwavelength effective medium layer, the complex acoustic problem involving multiple media interfaces is transformed into a solvable boundary value problem. The resulting model is discretized and solved using the finite element method, and the effects of unit-cell geometric parameters and array configurations on acoustic energy transmission are systematically investigated. Numerical simulation results demonstrate that the acoustic energy transmission efficiency of the system is significantly enhanced in the vicinity of the unit-cell resonance frequency, and that the finite element results show good agreement with the predictions of the equivalent medium model. The results validate the effectiveness of the proposed mathematical model and finite element numerical method in analyzing low-frequency cross-media acoustic transmission under millimeter-scale thin-ice and laboratory ideal conditions, providing a feasible numerical analysis tool for the design and optimization of related acoustic structures.

1. 引言

海洋覆盖了地球表面近四分之三，在水下及极地环境中，声波因衰减小、传播距离远而成为主要的信息传输与探测手段。然而，在空气–冰层–水体构成的跨介质体系中，由于介质物性参数差异显著，声阻抗严重失配，声波在界面处易发生强反射，导致声能透射效率显著降低，从而制约相关应用的发展[1]-[6]。

从数学角度看，跨介质声传输可归结为具有不连续系数的声学边值问题。界面参数突变及复杂结构使解析方法难以适用，需借助数值方法进行求解。有限元方法在处理多介质界面与复杂边界方面具有显著优势，已成为研究该类问题的重要工具。同时，基于等效介质理论，可将复杂微结构等效为具有特定声学参数的均匀介质层，从而降低建模复杂度并提高计算效率。

基于此，本文建立了跨介质声传输的等效介质模型与有限元仿真框架，系统分析了结构参数对声能传输特性的影响，并验证了数值结果与理论模型之间的良好一致性。研究结果为跨介质低频声传输问题的数值建模与有限元分析提供了可行方法。

2. 跨介质声传输问题

跨介质声传输是水声研究中的重要基础问题，在水下通信[7]、深海探测[8]及水下目标识别[9]等应用中具有关键意义。实际环境中，声波往往需要在空气、水、冰层及固体结构等多种介质之间传播。然而，由于不同介质在密度、声速及体积模量等物理参数上存在显著差异，其特征声阻抗往往相差数个数量级，导致声波在介质分界面处发生强烈反射，透射声能严重受限，从而显著降低跨介质声传输效率[10]。近年来，超表面通过在亚波长尺度内对界面等效声学参数进行调控，可实现对声阻抗与相位响应的人工设计，为改善跨介质声能传输效率提供了一种新的解决思路。图1给出了一种用于跨介质声调控的超表面结构示意图。

在空气–冰–水等多层介质体系中，冰层不仅具有明显的弹性和各向异性特征，其厚度、结构均匀性及内部裂隙分布还会随环境条件发生变化，使声波在传播过程中伴随复杂的模态转换与能量耗散[11]；在深海高压环境中，水体密度和材料力学性能随深度变化，声–固–流耦合效应显著增强，进一步加剧了声波跨界面传播的不确定性与复杂性。这些因素共同导致传统基于理想边界和均匀介质假设的声学模型在极端环境下适用性受限[12]。

受上述物理机制影响，跨介质声学系统中的声学参数在空间上呈现强不连续特征，使得传统声学材料与结构在宽频范围内难以实现有效的阻抗匹配，进而导致声能传输效率普遍较低。从数学角度来看，该问题对应于一类具有显著参数突变与多场耦合效应的声学边值问题，其解析求解与数值计算均面临较大困难。因此，发展合理的数学模型并引入稳定高效的数值计算方法，以刻画复杂介质界面条件下的声场分布与能量传输特性，已成为当前跨介质声学研究中亟需解决的重要应用数学与计算问题之一[13]-[18]。

Figure 1. An air-water acoustic metasurface across different media [13]

图1. 一种跨介质(空气–水)水声超表面[13]

3. 理论方法

有限元方法(Finite Element Method, FEM)是一种求解复杂物理场边值问题的数值计算方法。其基本思想是将连续的计算区域离散为有限数量的单元，并通过选取合适的插值函数对场变量进行近似表示，从而将原本难以解析求解的偏微分方程转化为代数方程进行数值求解[19]。

3.1. 有限元网格划分

在有限元数值计算中，网格划分是将连续计算区域离散为有限单元的关键步骤，其质量直接影响计算精度、收敛性以及计算效率[20]。合理的网格划分不仅能够准确刻画几何边界和物理场分布特征，还能够在保证计算精度的同时有效降低计算成本。

在有限元分析中，根据计算区域的几何特征及物理问题的复杂程度，可选用不同类型的网格与单元形式[21]。常见的有限元网格类型包括结构化网格与非结构化网格。结构化网格通常由规则的四边形或六面体单元构成，具有拓扑关系清晰、计算效率高等优点，适用于几何形状规则的计算区域；非结构化网格则多由三角形或四面体单元组成，具有较强的几何适应能力，能够灵活刻画复杂边界与局部细节，因此在超表面等复杂结构建模中得到广泛应用。

3.2. 网格质量

网格质量是影响有限元求解稳定性与精度的重要因素。单元形状应尽量规则，避免出现过度扭曲或高纵横比，以降低数值误差并改善矩阵条件数。通常通过网格无关性分析验证网格密度的合理性，即在逐步加密网格的过程中考察关键物理量的收敛性，从而确定满足精度要求的网格划分方案。

在商用有限元软件中，网格偏度常作为评价单元几何质量的重要指标，其数值越小，单元越接近理想形态，计算稳定性与精度越高。一般要求最大偏度不超过0.8，并避免高偏度单元集中分布。图2给出了COMSOL中网格偏度的可视化结果。由图2(a)与图2(b)对比可见，随着网格加密，单元偏度明显降低，计算精度相应提高，但计算成本也随之增加。

Figure 2. Distribution of mesh skewness under different mesh size conditions: (a) COMSOL built-in fine mesh; (b) COMSOL built-in extra-fine mesh

图2. 不同网格尺寸条件下的网格偏度分布：(a) COMSOL内置细化网格；(b) COMSOL内置超细化网格

4. 声传输模型分析

4.1. 声学控制方程

在频域条件下，线性声学问题通常可由声压$p\left(r,\omega \right)$ 的控制方程描述。对于无源、各向同性流体介质，声压场满足Helmholtz方程为[22]

$\nabla \cdot \left(\frac{1}{p}\nabla p\right)+\frac{{\omega }^2}{K}p=0$ (1)

其中，$r=\left(x,y,z\right)$ 表示声压在三维空间中的分布；$\rho$ 为介质密度；$K$ 为体积模量；$\omega =2\pi f$ 表示声波角频率；$f$ 为对应的频率。

为适应复杂几何与边界条件，式(1)通常通过加权残值法转化为弱形式[22]。引入测试函数$q$ ，并在计算域$\Omega$ 上积分，得到：

${\displaystyle {\int }_{\Omega }\frac{1}{p}\cdot \nabla p\cdot \nabla q\text{d}\Omega }-{\displaystyle {\int }_{\Omega }\frac{{\omega }^2}{K}pq\text{d}\Omega }={\displaystyle {\int }_{\Omega }qn\cdot \left(\frac{1}{p}\nabla p\right)\text{d}\Gamma }$ (2)

其中，$\Gamma$ 为边界，$n$ 表示外法向方向。根据有限元的核心思想，对计算域进行离散，并在每个单元内采用形函数$N_i\left(r\right)$ 对声压近似表示：

$p\left(r\right)\approx {\displaystyle \sum _{i=1}^N{N}_i\left(r\right)p_i}$ (3)

其中，$p_i$ 为第$i$ 个节点处的声压自由度，$N$ 为单元节点数。

将式(3)代入(2)，并对整个离散区域进行组装，得到如下线性代数方程组：

$\left(G-{\omega }^{2}M\right)P=F$ (4)

其中，$G$ 和$M$ 分别表示系统的刚度矩阵和质量矩阵；$F$ 为由边界条件或声源引入的激励项；$P={\left(p_1,p_2,p_3,\cdots ,p_N\right)}^{\text{T}}$ 为待求解向量，其表示计算域内声压场的数值近似结果。

4.2. 等效介质数学建模

我们首先建立一种理想化模型，在薄冰层下方引入平坦空气层，如图3(a)所示。假设冰层与水层厚度均远小于入射声波波长，在长波近似条件下，可将冰–水复合层等效为具有均匀参数的介质，如图3(b)所示。此时，等效介质层在声压激励下整体参与运动，而空气层主要表现为体积压缩与膨胀。在此基础上，原多介质声学问题可简化为一维等效模型，其动力学行为可等效为质量–弹簧系统：等效介质层对应集中质量项，空气层对应等效弹性项，如图3(c)所示。该集中参数模型忽略层内传播效应，显著降低了计算复杂度，为后续参数推导与数值分析提供了理论基础。

基于上述假设，该质量–弹簧系统的声阻抗$Z_d$ 可以用$x=-d$ 处的阻抗形式表示，其表达式为

$Z_d=\frac{Z_e}{1+\left(\frac{Z_e^2}{Z_{air}^2}-1\right){\sin }^2\left(k_{a}d\right)}+i\left(\frac{k_{a}h{\rho }_e{c}_a}{S}-\frac{Z_a}{2}\frac{\left(\frac{Z_e^2}{Z_{air}^2}-1\right)\sin \left(2k_{a}d\right)}{1+\left(\frac{Z_e^2}{Z_{air}^2}-1\right){\sin }^2\left(k_{a}d\right)}\right)$ (5)

其中，$S$ 表示气泡型结构的横截面积；$Z_e={\rho }_e{c}_e$ 和$Z_{air}={\rho }_a{c}_a$ 分别为等效介质和空气的声阻抗，${\rho }_e$ 、$c_e$ 以及${\rho }_a$ 、$c_a$ 分别表示等效介质和空气的密度与声速。等效介质的密度与体积模量可由冰层与水层的加权平均形式给出：

${\rho }_e=\frac{1}{h_1+h_2}\left(h_1{\rho }_{ice}+h_2{\rho }_w\right)$ (6)

$K_e={\left(\frac{h_1}{h_1+h_2}\frac{1}{K_{ice}}+\frac{h_2}{h_1+h_2}\frac{1}{K_w}\right)}^{-1}$ (7)

相应的等效声速为

$c_e=\sqrt{\frac{K_e}{\rho {}_e}}$ (8)

其中，${\rho }_{ice}$ 、$K_{ice}$ 以及${\rho }_w$ 、$K_w$ 分别表示冰和水的质量密度与体积模量。

由于冰层和水层的声阻抗均远大于空气的声阻抗[23]，在满足$Z_d=Z_{air}$ 的阻抗匹配条件，并结合约束关系$Z_e\gg Z_{air}$ 和$k_{a}d\ll 1$ 的低频近似下，系统的单位透射特征频率可表示为

$f_{\mu }=\frac{\omega }{2\pi }\sqrt{\frac{K_e}{M_e}}=\frac{c_a}{2\pi }\sqrt{\frac{{\rho }_a}{\left(h_1+h_2\right)d\rho {}_e}}$ (9)

其中$\omega$ 为声波的角频率，$M_e$ 表示等效介质层的等效质量。式(9)成立前提是

$\frac{h_e}{d}=\frac{h_1+h_2}{d}=\frac{c_e}{c_a}$ (10)

该关系表明等效介质层厚度与空气层厚度之间需满足特定的声速比例，以实现系统的有效阻抗匹配。对于一般情况，系统的声能传输系数$T$ 可表示为

$T=\frac{4}{2+\left(\frac{Z_{air}}{Z_{water}}+\frac{Z_{water}}{Z_{air}}\right){\cos }^2\left(k_{a}d\right)+\left(\frac{Z_e^2}{Z_{air}\cdot Z_{water}}+\frac{Z_{air}\cdot Z_{water}}{Z_e^2}\right){\sin }^2\left(k_{a}d\right)}$ (11)

对式(11)关于频率进行一阶求导分析发现，系统最大声能传输系数所对应的频率$f_{\max }$ 与(9)给出的特征频率$f_{\mu }$ 高度接近，表明该质量–弹簧模型能够有效刻画系统在低频区域的透射增强机理。

Figure 3. Schematic of the effective medium theory. (a) Schematic of the idealized bubble structure model. (b) Schematic of the effective medium model, where the effective thickness is denoted by $h_e=h_1+h_2$ . (c) Schematic of the mass-spring equivalent model, in which the bubble layer is modeled as an effective spring, while the thin ice layer and surrounding water are modeled as an equivalent mass block

图3. 等效介质模型示意图。(a) 理想结构模型示意。(b) 等效介质模型示意，其中等效厚度$h_e=h_1+h_2$ 。(c) 质量–弹簧模型示意图，其中气泡等效为弹簧，薄冰–水等效为质量块

在所关注的低频波段内，冰层结构的特征尺寸远小于声波波长，整体厚度处于深亚波长范围，使声场在结构尺度上变化平缓，满足等效介质理论的基本适用条件。同时，低频声波在冰–水复合体系中的传播主要受整体等效声阻抗和等效声速控制，局部结构引起的高阶散射效应可忽略，从而保证了等效介质模型的有效性。对比理论预测与数值仿真结果可发现，二者在透射系数及最大声传输频率方面具有良好一致性，进一步验证了在亚波长冰层条件下采用等效介质理论分析低频声学问题的合理性与可靠性。图4(a)~(c)给出了三个用于验证理论模型适用性的算例，其中冰层厚度分别取为10 mm、5 mm和3 mm。

Figure 4. Verification of the theory. (a) $h_1=10,h_2=5,a=30\text{ mm}$ , and $d=2.4\text{ mm}$ , yielding a theoretically predicted resonance frequency of $f_{\mu }=327.6Hz$ . (b) $h_1=5,h_2=7,a=30\text{ mm}$ and $d=2.25\text{ mm}$ , with a theoretically predicted resonance frequency of $f_{\mu }=373.1Hz$ . (c) $h_1=3,h_2=4,a=30\text{ mm}$ and $d=1.3\text{ mm}$ , corresponding to a theoretically predicted resonance frequency of $f_{\mu }=642.1Hz$

图4. 理论验证图(a) $h_1=10、h_2=5、a=30\text{ mm}$ 及$d=2.4\text{ mm}$ ，理论计算得$f_{\mu }=327.6Hz$ ；(b) $h_1=5、h_2=7、a=30\text{ mm}$ 及$d=2.25\text{ mm}$ ，理论计算得$f_{\mu }=373.1Hz$ ；(c) $h_1=3、h_2=4、a=30\text{ mm}$ 及$d=1.3\text{ mm}$ ，理论计算得$f_{\mu }=642.1Hz$

此外，冰层厚度对结构声学性能的影响主要体现在共振频率的变化上。由式(9)可知，随着薄冰厚度h₁的增大，共振频率逐渐减小，同时最大能量传输系数也随之发生变化。图5(a)~(b)给出了两种不同结构在薄冰厚度变化条件下的共振频率响应结果。可以看出，结构的共振频率随冰层厚度的增加而单调减小，与理论分析结果一致。

Figure 5. Verification of the theory. (a) A parametric sweep of h₁ under the condition $h_2=5,a=30\text{mm}$ and $d=2.4\text{mm}$ . (b) A parametric sweep of h₁ under the condition $h_2=7,a=30\text{mm}$ and $d=2.25\text{mm}$

图5. 理论验证图(a)在$h_2=5、a=30\text{mm}$ 和$d=2.4\text{mm}$ 条件下，参数化扫描h₁；(b) 在$h_2=7、a=30\text{mm}$ 和$d=2.25\text{mm}$ 条件下，参数化扫描h₁

4.3. 有限元仿真结果

在水介质中，气泡结构在低频条件下呈现显著的局域共振行为，其特征频率相较于同尺度的固体结构明显降低，从而在低频声学问题中具有重要的调控作用[23]。然而，从建模与数值计算的角度来看，裸露气泡在水中难以保持几何与力学稳定性，不利于建立可重复、可计算的数学模型。基于上述考虑，构建了用于有限元分析的气泡型结构单元，其几何构型如图6所示。该单元由固体外壳与内部方形空腔两部分组成，外壳边长为a，厚度为d，壁厚为c；方孔边长为b = a − 2c。该结构单元以晶格常数a在平面内进行周期性排列，从而形成规则的二维阵列。模型中，参数e是薄冰的厚度，f表示结构单元在水介质中的浸入深度。上述几何参数的引入使得该结构能够以参数化形式嵌入有限元计算框架，为后续的数值离散、参数分析及收敛性研究提供了统一的数学描述。

Figure 6. Schematic of the bubble-type unit cell

图6. 气泡型单元示意图

数值模拟采用COMSOL Multiphysics有限元软件实现。在空气区域上边界及水域下边界施加平面波辐射边界条件，以保证波动的无反射传播。在空气域内部施加沿−z方向的背景压力场作为入射声场。沿x、y方向设置周期性边界条件，实现单胞的空间重复性。固体–流体界面(空气–固体及水–固体)采用流固耦合边界条件，将流体声压变化与固体结构的弹性形变耦合，以模拟界面相互作用。网格生成采用四面体单元，确保各区域至少五层单元离散，且最大网格尺寸不超过a/25，以保证有限元解的精度和收敛性。

单胞的几何参数分别为：a = 3.5、b = 3.2、c = 0.15、d = 2.2、e = 5 mm以及f = 5 mm。空气和水的密度及声速如下：${\rho }_a=1.21{\text{ kg}/\text{m}}^3$ 和${\rho }_w=1000{\text{kg}/\text{m}}^3$ ，$c_a=343\text{m}/\text{s}$ 和$c_w=1500\text{m}/\text{s}$ 。根据材料力学的基本假设，可将小变形薄冰视为线弹性材料。其中冰和环氧树脂的密度、杨氏模量及泊松比分别为${\rho }_{ice}=900{\text{kg}/\text{m}}^3$ 、$E_{ice}=10\text{ GPa}$ ${\mu }_{ice}=0.31$ ；${\rho }_s=1100{\text{ kg}/\text{m}}^3$ 、$E_s=2\text{ GPa}$ ${\mu }_s=0.35$ 。空气和水的声阻抗分别为$Z_a={\rho }_a{c}_a$ 和$Z_w={\rho }_w{c}_w$ ，声能传输系数的数学模型为

$\tau =\frac{{\left|p_t\right|}^2{z}_a}{{\left|p_i\right|}^2{z}_w}$ (12)

其中，$p_t$ 、$p_i$ 分别为透射声压与入射声压。

通过空气–冰–水界面的传输增强的数学模型表示为

$T_e=20\lg \left(\frac{p_t^b}{p_t^0}\right)$ (13)

其中，$p_t^b$ 、$p_t^0$ 分别是在跨介质附近有或没有超表面的情况下的平均透射声压。

利用气泡在低频下共振特性，将气泡作为基本功能单元周期性嵌入界面结构中，构建气泡型声学超表面，从而实现对水下声能量传输特性的可控调节。图7(a)为气泡型超表面的结构示意图，该结构由气泡单元在x与y方向上周期性阵列构成。图7(b)为单个单元的网格划分示意图。

Figure 7. (a) Schematic diagram of the bubble-type metasurface; (b) Unit cell meshing

图7. (a) 气泡型超表面示意图；(b) 单元网格划分

结构的能量传输系数如图8所示，同时给出了未引入该结构时的参考结果。可以观察到，当在界面上引入含气泡的超表面后，约在461.7 Hz附近出现明显的透射峰。该透射增强主要源于气泡在该频率范围内与固体结构发生耦合振动，从而显著增强声能沿结构的传输效率。

Figure 8. Energy transmission coefficient of the metasurface without energy dissipation

图8. 超表面的能量传输系数(无能量耗散)

为进一步考虑实际环境中不可避免的能量耗散效应，本文在数值模型中引入了热–黏性声学模块。图9展示了该气泡型超表面结构的声能传输系数结果。结果表明，声波在传播过程中存在一定程度的能量耗散与转移，从而导致整体声能传输效率下降。与图8的无耗散情况相比，图9中传输系数峰值位置出现轻微偏移。该现象主要源于热–黏性模块对计算内存和资源的高需求。为提升计算效率，对求解过程中的时间步长进行了适当优化，虽在一定程度上降低了数值精度，但仍能够合理反映声能传输的物理趋势。

Figure 9. Energy transmission coefficient of the metasurface with energy dissipation

图9. 超表面的能量传输系数(有能量耗散)

5. 总结

跨介质声传输问题的本质在于不同介质之间存在显著的声阻抗不匹配，导致多介质界面处声能发生强烈反射，从而限制有效能量的耦合与传输效率。本文提出基于气泡型超结构的建模方案，将气泡作为低频共振单元周期性排列于界面上，并通过调控几何参数实现空气–冰–水体系中声能传输的增强，从而提供一种针对极端环境的可控声学调节策略。

数值分析表明，该方法的性能高度依赖气泡单元的共振特性，对冰层厚度及单元几何参数和工作频段具有较强敏感性。在考虑热–黏性耗散和结构稳定性等现实因素后，有限元数值模拟结果仍显示传输效率存在进一步优化的空间。目前研究主要基于有限元数值求解，实验验证及工程适用性仍需进一步探讨。上述分析明确了气泡型超结构在跨介质声传输建模中的关键参数、性能瓶颈及敏感性特征，为后续有限元优化设计与数值模拟方法提供了理论依据和方向参考。

基金项目

本研究由中国国家自然科学基金委员会(项目编号：12572103)和青岛市级自然科学基金(项目编号：24-4-4-zrjj-129-jch)资助。

NOTES

^*通讯作者。