COMSOL超声相控阵仿真:压力声学与固体力学模型选型与实现 做超声相控阵仿真这几年我见过太多人在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 是脉冲中包含的周期数一般取 35。比如 5 MHz、5 周期脉冲的持续时间是 1 µs。在 COMSOL 的解析函数里我会写成带延时参数的表达式方便每个阵元共用同一个函数、只传不同的延时进去A0*sin(2*pi*f0*(t-td))*((1-cos(2*pi*f0*(t-td)/N))/2)*(ttd ttdN/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 mmd 只要小于 0.59 mm 即可。很多阵列间距标的是钢中的半波长。如果你模拟的是水浸场景这点要特别小心水里看起来合理的间距可能已经低于半波长要求反而带来不必要的栅瓣风险。孔径大小则决定了波束的指向性和聚焦能力。孔径越大焦点附近的波束越窄近场区也越长。在仿真里孔径增大意味着阵元数量增加计算量直线上升所以通常先用 816 阵元做参数验证再扩展到全阵列。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 厚度设置在 12 个波长以上太薄会反射。边界层如果不加 PML反射波回来会和各种信号混在一起你在云图里看到的均匀背景条纹基本都是边界反射很容易被误判成散射信号。求解器方面压力声学的瞬态问题用 BDF 格式就挺稳。BDF 阶数建议 23超过 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)水10001480-碳钢785059003200铝270063203130PMMA119027301430钛4500607031254.2 缺陷加入与接收信号提取固体力学模型里缺陷的引入方式决定了仿真目标。如果模拟一个圆孔缺陷直接在几何里画一个圆设置圆边界为自由边界即可因为超声波入射到孔壁会产生反射、衍射。如果是模拟裂纹建议用一条细长椭圆或矩形槽来代替尖角裂纹更接近真实裂纹的回波特征网格也更好处理。接收信号的提取通常是在阵元所在边界上做积分平均。比如我想看第 8 个阵元收到的时间-位移曲线就在该边上添加探针或者全局计算输出速度场均方根或者位移的 y 分量随时间变化。然后在后处理里把接收信号按同一套延时表做延时叠加就能得到 A 扫和聚焦信号。这个过程中最容易犯的错是接收端的延时叠加方向和发射端是相反的。发射时远阵元要提前接收时要让远端到的信号延后对齐所以很多人在后处理里直接把发射延时表拿来用结果聚焦效果反而更差。我自己是习惯分开写两套延时数组发射用一套、接收用一套。4.3 裂纹/孔洞仿真的网格局部细化策略固体力学的网格要求比压力声学更苛刻因为它包含横波波长比同频率纵波短得多。同样 5 MHz钢中纵波波长 1.18 mm横波波长只有 0.64 mm。如果你按纵波波长划分网格横波会严重失真。网格划分时我用的策略是全局按横波波长控制缺陷附近再加密。全局最大单元取横波波长的 1/10大约 0.06 mm缺陷边界处加一个局部尺寸控制取到 0.020.03 mm。这样既保证波传播准确又不会让整个模型网格数爆炸。有一个技巧很实用在 COMSOL 里可以用自适应网格细化先跑一遍粗网格根据波场的梯度分布自动细化但超声瞬态问题我不太推荐全程自适应因为波前一直在动每步都重新划分网格会大幅增加耗时。更稳妥的做法是先用均匀网格跑通再针对缺陷区域手动局部加密。5. 两个模型的结果对比谁更适合你做的问题5.1 波场形态的差异纵波只给声压固体有横波把两个模型放在同一套阵列参数下跑完观察波场快照你会发现非常明显的差别。压力声学模型里只有标量声压的明暗条纹波前是干净的圆弧/倾斜直线看不到任何横波影子。这是因为流体里根本不存在剪切刚度这是物理本质决定的。如果你用压力声学去模拟钢块内部的斜入射检测会丢失波型转换这一关键信息可能会导致漏检。固体力学模型里位移矢量的分量云图会让你看到波型转换的完整过程纵波入射到自由表面或缺陷边界时会同时产生反射纵波和反射横波在斜入射时还会产生透射横波。这些波在云图里是不同方向、不同速度的多组波前叠加虽然看起来乱但恰恰是真实检测信号里会出现的回波。把焦点处的声压和位移幅度曲线对比会发现压力声学得到的聚焦信号更干净旁瓣水平更低固体力学的信号含有更多杂散波实际是横波和边界反射之间的相互作用。这个差异不一定是坏事它只是更接近真实。5.2 聚焦效果和缺陷响应的对比检验阵列参数是否正确最直观的方法是看焦点处有没有显著的能量汇聚。以偏转 30° 的线性阵列为例压力声学模型在焦点附近会形成一个椭圆形的声压极大区横向宽度大约 12 个波长固体力学模型在同样位置会看到位移幅值的聚集。实际对比后我建议把两套模型的声轴能量分布曲线画在一起在焦点附近垂直声轴方向截一条线提取声压幅值或位移幅值。两条曲线的主瓣宽度如果吻合说明延时表和阵列几何设置跨模型正确如果主瓣位置偏移先回头查延时表的声速是否统一。缺陷响应也有区别。同样一个 2 mm 圆孔压力声学模型里你看到的是孔表面的二次辐射源回波信号在 A 扫里是一个典型的双极性波形固体力学模型里孔洞产生的反射波里既有纵波也有沿表面爬行的表面波A 扫里会出现多个脉冲簇。用实验中真实探头信号对比固体力学的回波结构往往更贴近实测。5.3 计算代价与收敛情况这一点想提醒所有刚接触仿真的人固体力学的代价通常远超压力声学。压力声学每个节点只有 1 个自由度二维模型网格 30 万节点求解规模只有 30 万自由度固体力学每个节点 23 个自由度同样网格直接翻两三倍。再加上横波波长更短需要更细网格综合起来固体力学模型的自由度总数经常是压力声学的 510 倍。从我跑过的模型数据来看压力声学二维瞬态 20 µs 物理时间60 万自由度单台工作站大约 46 小时固体力学同样几何、同样时长自由度到 150 万直接要 12 天。所以如果你只是验证阵列聚焦性能优先用压力声学做缺陷定量分析再用固体力学仔细算。收敛性方面压力声学基本无脑收敛唯一需要留意的是 PML 层里网格太粗导致的高频毛刺。固体力学则要小心加载瞬间的冲击载荷激励开始的前几个时间步位移变化剧烈普通 BDF 容易不收敛。我一般给激励信号加一个平滑的斜坡起始区让初始冲击幅度从零缓慢上升能显著改善收敛窗口。6. 建模踩过的那些坑数值参数与常见误区6.1 网格太粗波形碎掉声速都算不准这是我反复提醒自己的一句话网格尺寸每增加一倍数值色散误差大概增加到四倍。网格太粗时5 MHz 的脉冲在传播 50 mm 后波形会从干净的汉宁窗正弦变成一串高低不平的锯齿波极大值和零点的位置全都偏移。检测方法很简单在波传播路径上放两个相距 20 mm 的探针点看两点之间波峰到达的时间差反算仿真声速。如果算出来声速和理论值偏差超过 1%网格就必须加密。这一步建议在正式跑全模型之前就做我用一个单独的网格收敛测试模型一分钟就能得到合适的网格尺寸。6.2 边界反射污染PML用少了的后果PML 太少会怎样我举个例子模型宽度 40 mmPML 厚度只有 0.2 mm结果在 15 µs 后波场里出现了一大片从左右边界卷回来的弧形波前正好出现在感兴趣区域把缺陷回波完全淹没了。PML 厚度原则上要超过一个中心波长实际做的时候我用 2λ 甚至 3λ同时保证 PML 内部网格沿外法向逐步拉伸。还有一种替代方案是用低反射边界条件但不适合大角度入射波我用过几次吸收效果都不理想。最终我个人还是倾向于 PML宁可多占一点计算域。6.3 延时单位与时钟错位的检查技巧COMSOL 里的时间单位问题直接能让仿真结果变得让你怀疑人生。有一次我在解析函数里写了t-tdtd 数值是 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 万。这种规模本地工作站基本跑不动就需要考虑集群或降维处理。这两套模型我到现在还会定期打开——毕竟每次换材料、换探头频率、换缺陷类型都要回到基础模型做验证。超声相控阵仿真的门槛并不高物理场选择、延时计算、网格和时间步控制这几关过了模型基本就稳了。剩下的无非是耐心调参数以及不断告诉自己看到诡异波形时先查网格再查延时最后查单位。这三板斧下来九成问题上都能定位。