做超声相控阵仿真这几年,我见过太多人在COMSOL里一上来就选物理场接口,结果算出来的波场一团糟,或者干脆算不动。其实相控阵仿真最关键的岔路口就是开头那一步:压力声学还是固体力学。这两个模型分别对应完全不同的物理假设、不同的计算代价、不同的适用场景。我手上这套COMSOL超声相控阵仿真模型,就是刻意把两条路线各做了一个版本放在一起,方便对比和选型。
这篇文章我会把两套模型的搭建思路、物理场设置、网格与求解器参数全部展开讲,包括我自己踩过的坑。无论你是做无损检测、医疗超声,还是声学超材料方向的仿真,看完应该都能少走不少弯路。
1. 为什么同样做相控阵,要拆成两个模型来做
1.1 压力声学与固体力学:一个管流体,一个管固体
先说本质区别。COMSOL里的压力声学(Pressure Acoustics)接口,求解的核心变量是声压 p,控制方程是声波方程,它默认介质是理想流体,只存在纵波,也就是压缩波。在这个接口里,你能看到的云图通常是声压分布,单位是帕斯卡(Pa)。
固体力学(Solid Mechanics)接口求解的核心变量是位移场 u,控制方程是弹性波方程,也就是纳维方程。它描述的是固体介质里的完整弹性波传播,包含纵波(P波)、横波(S波)、表面波、板波等等。云图通常是位移、速度或应力分布。
这就决定了它们的适用范围完全不同。我经常用一个类比:压力声学像是"在水面扔石子看涟漪",你只关心水表面波的能量怎么传;固体力学则像是"摇晃一整块果冻",内部的剪切形变、压缩形变都要管。果冻显然比水面复杂得多,计算量也大得多。
表格整理一下两者的核心差异:
| 对比项 | 压力声学 | 固体力学 |
|---|---|---|
| 求解变量 | 声压 p(标量) | 位移场 u(矢量) |
| 控制方程 | 亥姆霍兹/声波方程 | 弹性波纳维方程 |
| 介质假设 | 理想流体 | 各向同性/异性弹性固体 |
| 波类型 | 纵波(压缩波) | 纵波+横波+表面波 |
| 自由度/节点 | 1个 | 2D中2个,3D中3个 |
| 典型应用 | 水浸检测、医疗超声、空气声 | 焊缝检测、固体内部缺陷、岩石 |
1.2 实际检测场景里怎么选
实际工作中接到一个相控阵仿真任务,我一般先问一个问题:超声波的传播路径上,主要介质是流体还是固体?
如果是模拟水浸法检测,探头和工件之间有水层,超声从水进入钢再反射回来,那么从探头辐射出来的那一程,压力声学是最合适的建模方式。流体域里的声压云图能直接对应实验中水听器测到的信号。
如果是模拟接触法检测,比如压电晶片贴着钢块表面直接激发超声波,波绝大部分时间都在固体里传播,那必须用固体力学。因为在钢材中,纵波、横波、表面波可能会同时存在,波型转换也是重要的物理现象,压力声学完全描述不了这些。
回到这套模型本身:链接里给的两个模型,一个用压力声学,一个用固体力学,本质上就是让你在同一套阵列参数下分别看"流体中声场"和"固体中弹性波场"各是怎么回事。这种做法我认为很聪明,因为很多人做课题时根本不确定自己该用哪个接口,先把两条路都跑通,再回到自己的实际工况去选,效率反而最高。
2. 先算好延时再谈建模:阵列参数与激励信号准备
2.1 偏转延时的几何关系
相控阵的原理,说穿了就是"用时间换方向"。通过控制每一路阵元的激发时间,让波前在特定方向上同相叠加,实现声束偏转或聚焦。这个时间差的计算是所有后续模型搭建的前提。
假设阵元间距为 d,介质声速为 c,目标偏转角为 θ,那么第 n 个阵元相对于基准阵元的延时为:
Δt_n = n · d · sinθ / c
如果是聚焦,设焦点坐标为 (x_f, z_f),第 n 个阵元坐标为 (x_n, 0),那么它在焦点处的到达时间为:
t_n = √((x_f - x_n)² + z_f²) / c
每个阵元的延时就是 max(t_n) - t_n(让最远阵元先发,保证所有波前同时到达焦点)。这里注意一个反直觉的点:延时最大的是离焦点最近的阵元,因为它声程最短,需要延迟发出才能和其他阵元同时到达。很多新手第一次都会在这里翻车。
实际算延时,我通常直接用 Python 脚本生成一个延时表。以线性阵列 16 阵元、中心频率 5 MHz、阵元间距 0.6 mm、水中偏转 30° 为例:
import numpy as np c = 1480.0 # 水中的声速,m/s d = 0.6e-3 # 阵元间距,m theta = 30.0 # 偏转角,度 n_elem = 16 theta_rad = np.deg2rad(theta) delays = np.array([n * d * np.sin(theta_rad) / c for n in range(n_elem)]) delays -= delays.min() # 归一化到非负延时 for i, td in enumerate(delays): print(f"阵元 {i:2d}: {td*1e6:.4f} µs")这段代码跑出来,你会看到靠近基准端的阵元延时为 0,沿阵列方向线性增大。把这个表填到 COMSOL 的激励信号表达式里就行。
2.2 窗口函数和阵元激励信号
有了延时表,接下来是激励信号。做超声仿真,激励信号一般都用汉宁窗正弦脉冲。裸正弦波会在频域产生很宽的旁瓣,还会让波包在传播过程中严重拖尾;加窗之后信号频带收窄,波包更干净,仿真结果也更接近真实探头。
汉宁窗正弦脉冲的表达式写成:
A(t) = A₀ · sin(2πf₀t) · [1 - cos(2πf₀t / N)]
其中 N 是脉冲中包含的周期数,一般取 3~5。比如 5 MHz、5 周期脉冲的持续时间是 1 µs。
在 COMSOL 的"解析函数"里,我会写成带延时参数的表达式,方便每个阵元共用同一个函数、只传不同的延时进去:
A0*sin(2*pi*f0*(t-td))*((1-cos(2*pi*f0*(t-td)/N))/2)*(t>=td && t<=td+N/f0)这里 td 就是该阵元的延时。注意 COMSOL 里时间变量默认带单位,t 的单位是秒,所以如果延时量是微秒,需要写td*1e-6,这个单位问题非常容易出鬼,后面我会单独讲。
2.3 阵元间距、孔径与栅瓣控制
延时到位了,阵列几何设计也不能马虎。相控阵里有个必须满足的经验准则:阵元间距 d ≤ λ/2,否则会出现栅瓣,就是主瓣之外会莫名其妙多出几个大能量方向,检测时这些栅瓣会产生伪缺陷回波。
λ 是介质中的最小波长。以 5 MHz 探头为例,水里声速 1480 m/s,波长约 0.296 mm,所以 d 要小于 0.148 mm?实际工业相控阵探头很少这么密,那是因为在钢中纵波声速 5900 m/s,波长约 1.18 mm,d 只要小于 0.59 mm 即可。很多阵列间距标的是钢中的半波长。如果你模拟的是水浸场景,这点要特别小心:水里看起来合理的间距,可能已经低于半波长要求反而带来不必要的栅瓣风险。
孔径大小则决定了波束的指向性和聚焦能力。孔径越大,焦点附近的波束越窄,近场区也越长。在仿真里,孔径增大意味着阵元数量增加,计算量直线上升,所以通常先用 8~16 阵元做参数验证,再扩展到全阵列。
3. 模型一:基于压力声学的相控阵实现细节
3.1 控制方程与边界条件怎么落位
压力声学接口的瞬态控制方程,可以理解为以声压为变量的波动方程:
(1/c²) ∂²p/∂t² + ∇ · (-(1/ρ) ∇p) = 0
这个方程默认流体无粘、无流动、小振幅,绝大多数超声水浸场景都满足。你在 COMSOL 里新建模型时,选择"压力声学,瞬态"即可。
几何方面,我习惯建三个域:水浸区域(流体域)、试块区域(固体域,如果你关心声波进入工件后的传播)、外加最外层的完美匹配层(PML)。如果只关注探头辐射的声场本身,可以暂时不画试块,减少模型规模。
材料参数直接在材料库里选水就行,关键项是密度 1000 kg/m³、声速 1480 m/s。如果用水浸纵波到钢里的场景,钢的密度 7800 kg/m³、纵波速度 5900 m/s 也要在固体域设置好,然后在多物理场节点里耦合"声-固边界"。
3.2 每个阵元独立延时的COMSOL实现
压力声学接口里,激励源最常用的做法是给每个阵元所在的边界设置法向加速度边界条件。边界上一给法向加速度,就相当于在该处有一个活塞式振动源向外辐射声波。
具体操作:在模型树里选中某一个阵元对应的边界,添加"法向加速度"特征,然后把加速度表达式填成上面那段解析函数。COMSOL 里一个边界只能有一个激励条件,因此每个阵元都得单独添加一个特征,16 个阵元就是 16 个特征。
如果你用的是较新版本 COMSOL,可以用"阵列"功能或者把几何阵列与边界选择一起处理,减少重复操作。不过我个人仍然倾向逐个阵元手动添加特征,因为这样延时检错最直观——哪个阵元延时写错了,在模型树里一眼就能定位。
要注意激励信号导入方式。工程上常用任意波形发生器导出 CSV 波形,再在 COMSOL 里用插值函数导入。如果波形数据点不够密(一个周期少于 100 个采样点),仿真出来的波包会带明显高频抖动。我通常会在 Python 里把波形插值到至少每个周期 200 个点再导入。
3.3 网格与求解器设置:时间步、PML、稳定性
声学仿真的网格准则非常硬核:每个波长内至少要有 6 个二阶单元,追求稳定的话取 10 个以上。对应 5 MHz 水声(波长 0.296 mm),最大网格尺寸应设在 0.03 mm 左右。如果你按波长占比设了 6 个单元,波形相位误差还能接受;少于 5 个,你会发现声速都变了——波峰位置明显偏后,这就是数值色散。
时间步长的经验公式是:
Δt < Δx / c_max
实际操作中我会取更保守的值:Δt = Δx / (20 · c_max),也就是一个网格单元内波要走 20 个时间步。这样虽然计算时间长一点,但能避免很多莫名其妙的振荡。
PML 厚度设置在 1~2 个波长以上,太薄会反射。边界层如果不加 PML,反射波回来会和各种信号混在一起,你在云图里看到的"均匀背景条纹"基本都是边界反射,很容易被误判成散射信号。
求解器方面,压力声学的瞬态问题用 BDF 格式就挺稳。BDF 阶数建议 2~3,超过 4 容易引入数值高频振荡。直接求解器选 MUMPS,预处理方式为嵌套剖分,在二维模型里内存压力不大。
4. 模型二:基于固体力学的相控阵实现细节
4.1 弹性波方程与材料参数换算
固体力学里的控制方程写成张量形式是:
ρ ∂²u/∂t² = ∇ · σ + F
其中 σ 是应力张量,u 是位移矢量。在各向同性弹性固体里,弹性波有两个独立的传播速度,纵波速度 c_L 和横波速度 c_T 分别由杨氏模量 E、泊松比 ν 决定:
c_L = √(E(1-ν) / (ρ(1+ν)(1-2ν)))
c_T = √(E / (2ρ(1+ν)))
仿真时材料参数如果直接从手册里抄 E、ν、ρ,算出来的波速和实测值经常对不上。我的经验是先把目标声速确定,再反推 E 和 ν。比如要做钢中纵波 5900 m/s、横波 3200 m/s 的标准碳钢模型,用上面的公式反推,比直接抄一个"弹性模量 210 GPa"更可靠。
常见材料的声学参数参考如下:
| 材料 | 密度(kg/m³) | 纵波速度(m/s) | 横波速度(m/s) |
|---|---|---|---|
| 水 | 1000 | 1480 | - |
| 碳钢 | 7850 | 5900 | 3200 |
| 铝 | 2700 | 6320 | 3130 |
| PMMA | 1190 | 2730 | 1430 |
| 钛 | 4500 | 6070 | 3125 |
4.2 缺陷加入与接收信号提取
固体力学模型里,缺陷的引入方式决定了仿真目标。如果模拟一个圆孔缺陷,直接在几何里画一个圆,设置圆边界为自由边界即可,因为超声波入射到孔壁会产生反射、衍射。如果是模拟裂纹,建议用一条细长椭圆或矩形槽来代替尖角裂纹,更接近真实裂纹的回波特征,网格也更好处理。
接收信号的提取通常是在阵元所在边界上做积分平均。比如我想看第 8 个阵元收到的时间-位移曲线,就在该边上添加"探针"或者全局计算,输出速度场均方根或者位移的 y 分量随时间变化。然后在后处理里把接收信号按同一套延时表做延时叠加,就能得到 A 扫和聚焦信号。
这个过程中最容易犯的错是:接收端的延时叠加方向和发射端是相反的。发射时,远阵元要提前;接收时,要让远端到的信号延后对齐,所以很多人在后处理里直接把发射延时表拿来用,结果聚焦效果反而更差。我自己是习惯分开写两套延时数组,发射用一套、接收用一套。
4.3 裂纹/孔洞仿真的网格局部细化策略
固体力学的网格要求比压力声学更苛刻,因为它包含横波,波长比同频率纵波短得多。同样 5 MHz,钢中纵波波长 1.18 mm,横波波长只有 0.64 mm。如果你按纵波波长划分网格,横波会严重失真。
网格划分时,我用的策略是"全局按横波波长控制,缺陷附近再加密"。全局最大单元取横波波长的 1/10,大约 0.06 mm;缺陷边界处加一个局部尺寸控制,取到 0.02~0.03 mm。这样既保证波传播准确,又不会让整个模型网格数爆炸。
有一个技巧很实用:在 COMSOL 里可以用"自适应网格细化"先跑一遍粗网格,根据波场的梯度分布自动细化,但超声瞬态问题我不太推荐全程自适应,因为波前一直在动,每步都重新划分网格会大幅增加耗时。更稳妥的做法是先用均匀网格跑通,再针对缺陷区域手动局部加密。
5. 两个模型的结果对比:谁更适合你做的问题
5.1 波场形态的差异:纵波只给声压,固体有横波
把两个模型放在同一套阵列参数下跑完,观察波场快照,你会发现非常明显的差别。
压力声学模型里只有标量声压的明暗条纹,波前是干净的圆弧/倾斜直线,看不到任何横波影子。这是因为流体里根本不存在剪切刚度,这是物理本质决定的。如果你用压力声学去模拟钢块内部的斜入射检测,会丢失波型转换这一关键信息,可能会导致漏检。
固体力学模型里,位移矢量的分量云图会让你看到波型转换的完整过程:纵波入射到自由表面或缺陷边界时,会同时产生反射纵波和反射横波;在斜入射时还会产生透射横波。这些波在云图里是不同方向、不同速度的多组波前叠加,虽然看起来乱,但恰恰是真实检测信号里会出现的回波。
把焦点处的声压和位移幅度曲线对比会发现:压力声学得到的聚焦信号更干净,旁瓣水平更低;固体力学的信号含有更多杂散波,实际是横波和边界反射之间的相互作用。这个差异不一定是坏事,它只是更接近真实。
5.2 聚焦效果和缺陷响应的对比
检验阵列参数是否正确,最直观的方法是看焦点处有没有显著的能量汇聚。以偏转 30° 的线性阵列为例,压力声学模型在焦点附近会形成一个椭圆形的声压极大区,横向宽度大约 1~2 个波长;固体力学模型在同样位置会看到位移幅值的聚集。
实际对比后,我建议把两套模型的"声轴能量分布"曲线画在一起:在焦点附近垂直声轴方向截一条线,提取声压幅值或位移幅值。两条曲线的主瓣宽度如果吻合,说明延时表和阵列几何设置跨模型正确;如果主瓣位置偏移,先回头查延时表的声速是否统一。
缺陷响应也有区别。同样一个 2 mm 圆孔,压力声学模型里你看到的是孔表面的二次辐射源,回波信号在 A 扫里是一个典型的双极性波形;固体力学模型里孔洞产生的反射波里既有纵波也有沿表面爬行的表面波,A 扫里会出现多个脉冲簇。用实验中真实探头信号对比,固体力学的回波结构往往更贴近实测。
5.3 计算代价与收敛情况
这一点想提醒所有刚接触仿真的人:固体力学的代价通常远超压力声学。压力声学每个节点只有 1 个自由度,二维模型网格 30 万节点,求解规模只有 30 万自由度;固体力学每个节点 2~3 个自由度,同样网格直接翻两三倍。再加上横波波长更短需要更细网格,综合起来固体力学模型的自由度总数经常是压力声学的 5~10 倍。
从我跑过的模型数据来看,压力声学二维瞬态 20 µs 物理时间,60 万自由度,单台工作站大约 4~6 小时;固体力学同样几何、同样时长,自由度到 150 万,直接要 1~2 天。所以如果你只是验证阵列聚焦性能,优先用压力声学;做缺陷定量分析,再用固体力学仔细算。
收敛性方面,压力声学基本无脑收敛,唯一需要留意的是 PML 层里网格太粗导致的高频毛刺。固体力学则要小心加载瞬间的冲击载荷,激励开始的前几个时间步位移变化剧烈,普通 BDF 容易不收敛。我一般给激励信号加一个平滑的斜坡起始区,让初始冲击幅度从零缓慢上升,能显著改善收敛窗口。
6. 建模踩过的那些坑:数值参数与常见误区
6.1 网格太粗:波形碎掉,声速都算不准
这是我反复提醒自己的一句话:网格尺寸每增加一倍,数值色散误差大概增加到四倍。网格太粗时,5 MHz 的脉冲在传播 50 mm 后,波形会从干净的汉宁窗正弦变成一串高低不平的锯齿波,极大值和零点的位置全都偏移。
检测方法很简单:在波传播路径上放两个相距 20 mm 的探针点,看两点之间波峰到达的时间差,反算仿真声速。如果算出来声速和理论值偏差超过 1%,网格就必须加密。这一步建议在正式跑全模型之前就做,我用一个单独的网格收敛测试模型,一分钟就能得到合适的网格尺寸。
6.2 边界反射污染:PML用少了的后果
PML 太少会怎样?我举个例子:模型宽度 40 mm,PML 厚度只有 0.2 mm,结果在 15 µs 后,波场里出现了一大片从左右边界"卷"回来的弧形波前,正好出现在感兴趣区域,把缺陷回波完全淹没了。
PML 厚度原则上要超过一个中心波长,实际做的时候我用 2λ 甚至 3λ,同时保证 PML 内部网格沿外法向逐步拉伸。还有一种替代方案是用"低反射边界"条件,但不适合大角度入射波,我用过几次吸收效果都不理想。最终我个人还是倾向于 PML,宁可多占一点计算域。
6.3 延时单位与时钟错位的检查技巧
COMSOL 里的时间单位问题,直接能让仿真结果变得让你怀疑人生。有一次我在解析函数里写了t-td,td 数值是 0.1016(我本来以为在 COMSOL 里数值就用秒),实际上我想表达的是 0.1016 µs,结果所有阵元的波形全都挤在一起,波束聚焦效果完全消失。
后来我养成一个习惯:所有延时变量统一用微秒标注,在表达式里显式乘以1e-6,并在文件名里写清单位。比如延时表列头就写td_us,值填 0.1016,在 COMSOL 里引用时写td_us*1e-6。这样即使隔了几个月回来再打开模型,也不会因为单位问题重新踩坑。
6.4 二维与三维的边界条件差异
很多人把二维模型的结果直接对应三维实验,这是个大误区。二维模型在 COMSOL 里默认是"面外无限延伸"的线源,声波会以柱面波形式扩散;三维模型才是真实的点源聚焦,几何衰减规律完全不同(柱面波 1/√r 衰减,球面波 1/r 衰减)。
如果最终目标是三维仿真,我建议分两步走:先用二维把阵列延时、网格参数、缺陷尺寸等基本参数验证对,再映射到三维。三维模型的计算量不是线性增加,是立体爆发的。以 16 阵元阵列、5 MHz 钢块 20 mm×20 mm×20 mm 为例,网格要 600 万节点以上,自由度轻松破 2000 万。这种规模本地工作站基本跑不动,就需要考虑集群或降维处理。
这两套模型我到现在还会定期打开——毕竟每次换材料、换探头频率、换缺陷类型,都要回到基础模型做验证。超声相控阵仿真的门槛并不高,物理场选择、延时计算、网格和时间步控制这几关过了,模型基本就稳了。剩下的无非是耐心调参数,以及不断告诉自己:看到诡异波形时,先查网格,再查延时,最后查单位。这三板斧下来,九成问题上都能定位。