Comsol相控阵聚焦声压分布仿真：相位补偿到网格剖分全流程解析

做相控阵声场仿真时我踩过最冤枉的一个坑是相位明明算对了但声压分布图上焦点位置就是不对——后来才发现是Comsol默认时谐约定在捣鬼。这篇文章想把Comsol相控阵聚焦声压分布仿真的完整脉络理顺从阵元相位怎么算、边界条件怎么加到网格怎么剖、结果怎么读把那些文档里不会写的细节一次讲清楚。适合刚接触声学模块、想在仿真里验证聚焦性能的人也适合做超声无损检测、医学超声或者声镊的同学拿来当参考。1. 相控阵聚焦的核心逻辑相位补偿与等效球面波前1.1 聚焦的本质是让所有阵元的声波“同时到达”相控阵聚焦听起来很高端本质上就是一句话把每个阵元的发射时间调好让它们发出的声波在目标位置同相叠加。各阵元到焦点的距离不一样声波走的路径就有长有短到达焦点的时间自然有先有后。只要给每个阵元加上一个与声程差对应的相位延迟让所有波前在焦点处对齐你就会看到那个位置的声压明显抬起来其它位置则因为相位不齐而互相抵消。以常见的线性阵列为例阵元沿x轴排布声束沿y轴传播。假设第i个阵元中心坐标是(x_i, 0)焦点坐标是(x_f, y_f)那么该阵元到焦点的距离就是r_i sqrt((x_i - x_f)^2 y_f^2)。为了补偿声程差通常取一个参考距离R0然后给每个阵元加相位phi_i k * (sqrt((x_i - x_f)^2 y_f^2) - R0) k 2 * pi / lambda 2 * pi * f0 / c0这里的R0是任意参考常数常用中心阵元到焦点的距离也可以直接取0。加这个常数只影响绝对相位不影响各阵元之间的相对补偿关系所以完全不用纠结。把相位设计成这个形式之后焦点处实际等效于每个阵元都发射了一个球面波前在焦点处恰好同相形成Constructive Interference。这里有个新手最容易犯的错把相位公式写成phi_i k * R0 - k * r_i符号反了。频域仿真里相位差是正还是负取决于Comsol的时谐约定。也就是下文要专门说的问题。实际判断方法很简单算完之后看焦点处的声压如果是亮斑就对了如果是暗斑把相位表达式里的负号反过来就行。1.2 频域声学模型稳态聚焦问题的最佳选择很多人在做超声相控阵仿真时第一反应是上时域瞬态求解觉得能看到声波传播过程才直观。但如果是纯聚焦问题研究的核心是“某个单频激励下声场长什么样”那Comsol的压力声学频域接口就是最合适的选择。它直接求解Helmholtz方程一次求解就能得到整个计算域内的复声压分布计算量比瞬态小几个数量级。频域接口输出的p是复数实部和虚部合在一起包含了幅值和相位信息。你画声压云图时用abs(p)看幅值需要相位图就atan2(imag(p), real(p))非常方便。瞬态仿真当然也能做但你要么用脉冲激励看波前传播要么研究非线性声场否则没必要给自己加负担。还有一个容易被忽略的点频域求解假设系统是线性的稳态响应。对于医学超声里常见的非线性传播、组织中的谐波成像频域线性模型是不够的那需要压力声学-非线性或者时域求解器。但本文讨论的“聚焦声压分布探索”频域线性模型完全够用而且结果非常干净便于验证相位补偿逻辑是否正确。2. 从几何到激励Comsol建模全流程落地2.1 二维几何搭建阵元排布、计算域与PML我建议先从二维模型开始。相控阵的阵元沿直线排列声场在包含阵列轴线和传播方向的平面上基本是二维的用2D模型能快速验证相位算法等方案成熟了再扩展到3D。下面以16阵元、中心频率1 MHz、水中声速1480 m/s的模型为例说清楚每个参数怎么定。几何尺寸这样设置计算域取x方向从-25 mm到25 mmy方向从0到60 mm。阵元表面位于y0这条边上沿x轴排布阵元宽度w0.8 mm阵元中心间距pitch1 mm共16个阵元分布在x-8 mm到x8 mm左右。焦点放在(0, 40 mm)处也就是声束正前方40 mm。这个配置很常见孔径16 mm左右焦距40 mm算是中等聚焦。计算域四周不能直接设成硬边界否则声波会在边界上来回反射云图里全是干涉条纹完全没法看。标准做法是在左右两侧和远端的下侧这里是y60 mm一侧加完美匹配层PML。PML的作用是把入射波“吃掉”模拟无限大水域。厚度我建议至少取半个波长1 MHz水中波长约1.48 mm所以PML厚度设个2 mm比较稳妥。如果你的计算频率会往下扫比如最低0.8 MHz那波长变成1.85 mmPML厚度就要按最低频率对应的波长来留不然低频成分会在边界处被弹回来。2.2 材料参数与物理场设置域内材料用水的参数即可密度ρ998 kg/m³声速c01480 m/s。打开“压力声学-频域”接口把整个矩形域都设为压力声学域然后在“域设置”里把声速和密度填进去。Comsol内置材料库里也有Water但我习惯直接手填少一层耦合免得材料节点干扰判断。这里要特别说一下理论上声速和密度应该对应到域方程里但如果你将来想把模型改成油、酒精或者其他液体记住声速变了波长也变网格尺寸要跟着改。最容易翻车的就是改了声速忘了重新剖网格结果高频算出来全是数值色散。物理场设置好后需要把PML区域从普通声学域中区分出来。做法是在“定义”里添加“完美匹配层”节点几何选择PML区域坐标类型选直角坐标。添加PML之后Comsol会重写该区域内的控制方程你不需要手动改任何物理场设置。2.3 边界条件一个表达式同时给16个阵元加相位这是整个建模流程里最有技巧的一步。很多人一开始会把16个阵元拆成16个边界条件逐个写相位麻烦不说阵元一多改起来就是灾难。正确做法是把16个阵元表面统一选中做成一个“显式选择”然后施加一个“指定法向加速度”边界条件表达式里直接使用坐标变量x和y让相位随边界上的坐标自动变化。具体操作是在几何里选中y0边上属于阵元的16条线段创建显式选择。然后在压力声学物理场中添加“法向加速度”边界条件选择这个显式选择法向加速度表达式写A0 * exp(i * k0 * (sqrt((x - x_f)^2 (y - y_f)^2) - R_ref))参数设置如下参数名值说明f01 MHz中心频率c01480 m/s水中声速lambdac0 / f0波长k02*pi / lambda波数x_f0 mm焦点x坐标y_f40 mm焦点y坐标R_ref40 mm参考距离可取任意合理值A01 m/s²阵元表面法向加速度幅值这个表达式之所以可行是因为Comsol在求值时会遍历边界上每个高斯点代入该点的x、y坐标。阵元表面的x不同算出来的相位就不同恰好实现“每个阵元按坐标自动补偿相位”。整个过程只需要一个边界条件后续想改成32阵元、64阵元只要几何和显式选择更新表达式一个字符都不用动。有个细节y0边上的阵元之间缝隙区域也要处理。真实换能器的阵元之间是惰性材料基本不振动在声学上视为硬边界。 Comsol里对于外部边界如果不加任何条件默认就是声学硬边界也就是法向速度为0。所以只要把阵元表面的显式选择做对其它区域自动保持硬边界不需要额外设置。2.4 网格划分与求解器选择声学仿真的网格直接决定结果准不准。我自己的经验是压力声学频域问题最大网格尺寸至少要达到λ/6保守一点用λ/8。对于1 MHz水声λ/6约0.25 mm在40 mm x 60 mm的矩形域里三角形网格数量大概几万直接用“物理场控制网格”里的“超细化”档位就能得到不错的结果。但物理场控制网格有个隐藏问题它对PML区域不一定生成理想的分层结构。PML本身需要网格沿吸收方向逐渐变化最好是扫描网格或者较密的规则网格。我的做法是手动剖分先选PML区域用映射网格把厚度方向分成4到6层再对核心计算域用自由三角形网格。这样既保证PML吸收效果又控制核心区域网格密度。求解器方面频域问题用直接求解器就行。我一般选PARDISO鲁棒性比MUMPS好内存占用也能接受。千万别为了省内存去用迭代求解器声学频域问题的矩阵性质不太好迭代很容易不收敛你会在收敛曲线图上浪费大量时间。3. 声压分布结果解读焦点、旁瓣与焦斑尺寸3.1 声压云图怎么读求解完成后二维绘图组里画表面图表达式填abs(p)就能看到声压幅值分布。如果相位补偿方向正确你应该在(0, 40 mm)附近看到一个明显的亮斑这就是焦点。亮斑周围是一圈圈近似椭圆的等高线向外逐渐变暗。焦点前方远场方向会出现一系列平行弯曲的波前焦点后方则是复杂的干涉区亮度较低。拿到云图第一步不是急着赞叹而是用“交叉线定位”核实焦点位置。我通常创建一个“派生值-点计算”输入坐标(0, 40 mm)再在焦点前1 mm和后1 mm分别取点看声压是不是中间最大。如果最大值不在(0, 40 mm)先检查相位表达式里的x_f、y_f是不是写错单位再检查阵元选择是否漏了几条边。有一个非常隐蔽的坑pitch中心间距如果比阵元宽度大很多阵元之间的缝隙也会显著影响聚焦效果但缝隙在云图上不容易看出来最好在几何图上用显式选择高亮检查一遍。云图还能大致看出旁瓣情况。焦点中心越亮、周围越暗说明能量集中度越好。如果你看到焦点周围有一条条明显的放射状亮纹那基本上是旁瓣或者栅瓣。栅瓣通常是因为阵元间距超过半波长引起的后面会专门说。3.2 沿轴与横向声压曲线光看云图只能定性判断定量分析要画沿轴声压分布和横向声压分布曲线。沿轴曲线是把y从0到60 mm画一条直线数据集x固定为0然后绘制abs(p)随y的变化。这条曲线能看到声波从阵元表面出发在焦点处达到峰值之后再缓慢衰减。焦点位置附近你会看到声压从0迅速上升到峰值再下降峰值半高宽就是焦区长度。这类曲线验证聚焦行为非常直观也是判断焦点是否偏的最硬指标。横向曲线更有意思。在焦点所在的y40 mm处创建一条水平线数据集x从-20 mm到20 mm绘制abs(p)。你会看到中心有一个主峰两侧出现一系列次级峰这就是旁瓣。均匀激励线性阵的第一旁瓣理论上约-13.5 dB左右但由于聚焦几何和孔径效应仿真值会略有出入。我第一次跑出来的横向曲线主峰旁边还叠了个小鼓包排查半天发现是左右PML厚度不对称。这个后面在排查表里会再提。3.3 焦点增益与焦斑尺寸的对照分析读曲线的时候有两个数字值得记录焦点声压幅值相对单阵元中心声压的增益以及-6 dB焦斑尺寸。理想情况下N个阵元同相叠加焦点声压理论上接近单阵元贡献的N倍。仿真里由于各阵元到焦点的距离不完全相等传播衰减略有差异实际增益会略低于N。16阵元模型我实测大概在12到14倍之间。如果你算出来只有几倍说明相位补偿有问题或者阵元间距过大产生了栅瓣抵消。焦斑尺寸可以跟理论公式对比。聚焦相控阵横向-6 dB焦斑宽度近似为D_f 0.886 * lambda * R / DR是焦距40 mmD是孔径约16 mmλ1.48 mm代入算出来约3.3 mm。这个量级跟仿真结果非常接近可以作为模型验证的一个标尺。另外你会注意到焦点处的声压不是无限尖焦斑大小取决于阵列孔径和频率。想把焦斑做小要么提高频率要么增大阵元数量撑大孔径——这就是相控阵设计的基本权衡。4. 参数扫描、相位误差与问题排查4.1 参数化扫描焦点从20 mm扫到50 mm静态模型跑通之后最自然的探索方向是参数扫描。我想看焦点距离变化对声压分布的影响于是把y_f设为参数从20 mm到50 mm以5 mm步长扫描。这里的关键是R_ref保持40 mm不变也不会影响结果因为它是个常数偏移。在“研究”里添加“参数化扫描”选择参数y_f填入20, 25, 30, 35, 40, 45, 50 mm运行求解。完成后数据集中会多出每个扫描解。你可以在“一维绘图组”里选择所有解一次性画出不同焦距下的沿轴声压曲线。一个明显趋势是焦距越近焦点峰值越高焦斑越窄焦距越远峰值下降焦斑变宽。原因是传播距离更长几何衰减更大等效孔径对应的聚焦能力变弱。更实用的做法是用“派生值-组计算”在多个扫描解上计算同一个点(0, y_f)的abs(p)会得到一个表格直观看出焦点声压随焦距的变化。我实际跑出来的数据里20 mm焦距的焦点声压比50 mm焦距高了接近一倍这说明短焦距聚焦更“锐”但穿透深度小长焦距穿透深但聚焦增益打折。4.2 相位误差的敏感性实验工程实际中FPGA或模拟电路给阵元加相位时不可能完全精确相位误差会直接影响聚焦性能。这个可以用仿真快速验证在相位表达式里额外加一个随机项delta_phi。做法是定义一个参数error_amp然后相位表达式改为A0 * exp(i * (k0 * (sqrt((x - x_f)^2 (y - y_f)^2) - R_ref) error_amp * sin(100*x)))这里用sin(100*x)当做一个空间变化的伪随机扰动模拟不同阵元的相位偏差。把error_amp从0扫到0.5 rad你会看到焦点峰值逐渐下降旁瓣逐渐抬升。当相位误差超过0.5 rad时横向曲线上的旁瓣可能比主瓣只低6到8 dB聚焦效果已经显著恶化。这个结论对做硬件的人特别有价值它告诉你相位分辨率要做到多少才够。如果每档相位步进太大引入的系统误差会直接映射成旁瓣上升。在工程上相位量化误差一般控制在π/8以内对应的旁瓣恶化是可以接受的。4.3 常见问题排查表我把实际调试过程中遇到最多的问题整理成下表按“现象-原因-解决”三条线列出来现象可能原因解决方式焦点处是暗斑而非亮斑相位符号反了与Comsol时谐约定不匹配将exp(iphi)改为exp(-iphi)焦点位置偏移不在设计坐标x_f/y_f单位写错或部分阵元未选中检查参数单位重现检查显式选择范围云图出现大量横向条纹PML太薄或网格太粗边界反射增加PML厚度到0.5λ以上细化网格焦点增益远低于理论值阵元间距过大产生栅瓣检查pitch是否大于0.5λ适当增加阵元数扫描焦距时结果突变某个扫描点焦点落在PML附近将y_f扫描上限减小留出至少2λ缓冲计算内存不足网格过密或3D模型过大先用2D验证或利用对称性建半模型表格里最后一条的对称性值得单独说当焦点在阵列中心轴上、阵列几何也对称时相位表达式关于x0是对称的整个声场关于x0对称。这时候只建x≥0的半边模型在x0边界上设对称条件计算量直接减半。我后期做参数扫描时一直用半模型速度提升了接近一倍。4.4 关于时谐约定的一个重点提醒这是我踩过最大的坑值得单独强调。Comsol的“压力声学-频域”默认时谐因子是exp(iωt)但不同模块、不同版本可能不一致。当你用复数相位exp(iφ)做激励时如果软件内部用的是exp(-iωt)约定那么你算出来的“焦点”实际上是相消的暗斑。判断方法很简单算完后看焦点声压。如果云图里对应位置是暗的而周围反而亮那基本就是相位符号问题。这种情况不需要改公式只要把指数项里的i换成-i也就是把相位取反重新求解即可。千万别去改全局时谐设置容易影响其它结果理解。5. 一些实操心得与扩展思路5.1 我踩过的几个坑第一个坑是把网格全交给物理场控制。刚开始图省事结果PML区域全是不规则三角形吸收效果很差云图边缘全是反射条纹。后来改成PML用映射网格分5层核心计算域用自由三角形干净很多。第二个坑是焦点坐标单位。在Comsol里几何默认单位是m参数表里填40 mm显示没问题但在相位表达式里的x_f如果忘写单位默认按m理解等于把焦点放到了40米远云图上当然什么都没有。这种错误最坑人因为Comsol不报错你反复检查边界条件都发现不了。第三个坑是参数化扫描时用了数组表达式扫描了20个点结果PML厚度是按最低频率波长设计的高频没问题低频出反射。后来我统一按最低频率波长设置PML厚度多花了一点网格量但所有扫描点一次通过。5.2 还能往哪些方向扩展这个模型验证过关后拓展空间很大。最直接的是改成32阵元或64阵元观察焦点增益和焦斑尺寸的变化趋势也可以研究不同阵元间距对栅瓣的影响这在无损检测里特别重要。更贴近实际应用的做法是把阵元从“法向加速度”边界换成压电换能器结构耦合模型用压电材料外加电压激励模拟真实的机电转换过程。那样模型会变成压电-结构-声场三物理场耦合计算量上一个大台阶但能看到实际阵元振动幅度和相位畸变对聚焦的影响。如果做医学超声方向还可以把焦点附近区域加上生物组织域设置密度、声速、衰减系数观察衰减对聚焦深度的影响。再进一步耦合生物传热方程就能评估HIFU治疗的热剂量分布。这些都需要在此基础上逐步加模块但核心的相控阵聚焦声压分布探索逻辑是不变的。我个人实际操作中的体会是相控阵聚焦仿真最花时间的往往不是求解本身而是相位补偿方向、单位规范和边界条件选择这些“看起来很小”的点。如果你也想做类似探索建议先把单个阵元模型跑通确认声压分布正常再逐步加阵元数量。单阵元都错了16个阵元叠加只会让错误更隐蔽。另外记得多保存几个阶段版本我每次修改相位公式或者网格策略之前都会复制一份模型文件因为调试到后面你根本分不清是哪个参数改出来问题。