Page 221 - 《应用声学》2024年第1期

第 43 卷 第 1 期 孔惠元等： 圆柱类构件超声相控阵聚焦模型 217

θ 时折射声束无法进入构件，设此时入射声束与界 fsolve 函数是通过选择不同的参数来选取不同
面两侧交点为 p、q，则有效声束入射点位于p 与q 之 的算法，如牛顿迭代法或 Levenberg-Marquardt 算
间 (这里不考虑波束角大小和波型转换现象，假设 法，再根据初始值和方程组的导数信息逐步接近方
每个阵元发射声束可以覆盖整个构件，以纵波声线 程组的解，本文通过选择 Levenberg-Marquardt 算
路径进行研究，不考虑横波影响)。设第 i 个阵元位 法寻找未知变量的值。Levenberg-Marquardt 算法
置坐标 (x i , 0)，x i 满足式 (1)，设临界点 p 位置坐标 是使用最广泛的非线性最小二乘法，它的关键是用
(x p , y p )，点q 坐标(x q , y q )，由折射定律有 模型函数对待估参数在其领域内做线性近似，忽略
c 1 掉二阶以上的导数项，从而转化为线性最小二乘问
sin θ = . (20)
题，它具有收敛速度快的优点。此外使用fsolve函数
c 2
同时入射角 θ 为向量 po(−x p , h + r − y p ) 和向 时提供一个合适的初始值非常重要，初始值选择不
−→
− →
量 pi(x i − x p , −y p )夹角的补角，即入射角θ 满足 合理可能导致算法无法收敛。初始值的确定与阵元

− → −→ − → −→ 坐标、聚焦点坐标、检测对象参数、坐标系建立方式
po × pi = |po| · pi · sin(π − θ). (21)
等有关，结合实际相控阵换能器各阵元中心间距较
−→
将向量 po、向量 pi 坐标代入，得
−→
小(本文选用中心间距0.6 mm相控阵换能器)，则不
x p y p − (x i − x p )(h + r − y p )
sin θ = √ . (22) 同阵元对于初始值的选取影响较小，且本文以换能
2
r (x i − x p ) + y 2
p 器中心建立坐标系，对于左右两侧的阵元，入射点分
将 sin θ 向量坐标表达式代入式 (20)，并与点 p 别位于构件中心左右两侧，因此本文选择待检测圆
所在曲率圆上方程联立得方程组 柱面上超声换能器阵列中心对应位置 (0,5) 作为初

x p y p − (x i − x p )(h + r − y p ) c 1 始值。

 √ = ,

r (x i − x p ) + y p 2 c 2 (23) 基于费马定律的迭代算法计算延迟时间时，入
2

 射点的坐标精度取决于界面的离散精度，离散精度
x + (y p − h − r) = r , y p 6 h + r.
 2 2 2
p
越小，计算精度越高，但计算时间越长。分析本文的
通过解算式 (23) 即可求得 p(x p , y p )，在 y p 6 h + r
推导过程与方程组的求解过程可知，本文方法避免
范围内临界点 p、q 坐标同时满足式 (23) 方程组，因
了费马定律遍历算法中界面离散精度的影响。然而
此可以解出两组实根，即为临界点p、点q 坐标。
在方程组求解过程中，所选算法不同、初始值选取
2.3 入射点方程组求解 不同也会影响方程组求解的精度与效率。实际检测
在上述延迟时间的计算中，重点在于对非线性 过程中，受阵元坐标测量误差、换能器频率存在微
方程组式 (13) 的求解，目前求解非线性方程组的方 小偏差导致近场区计算误差、构件声速与理论值存
法有很多，如直接降维法、最小二乘法、牛顿迭代法 在偏差等因素影响，也会导致声束聚焦点坐标实际
等，本文使用数据处理软件中相关工具包对方程组 值与理论值存在微小偏差。
进行编程求解，算法思路及核心函数如图3所示。
3 仿真验证与结果分析
Ύၹ
ࠀ˧ԫ᧚nj MatlabFunction
Ԡ஝ԣவሮ Ѧ஝ၷੇѦ஝Բ౻
3.1 仿真条件
为了验证本文方法在计算效率上的优势与有
ࠀ˧Ѻݽག֗ Ύၹfsolve 效性，本文利用数据处理软件进行数据计算，利用
රᝍ٨ᤥᮊ Ѧ஝රᝍ
有限元仿真软件对如图 4 所示的圆柱类构件进行仿
真。模型结构的长 30 mm、高 20 mm，上半部分为
Ύၹoptimoptions ᣥѣරᝍϙnj 耦合介质，下半部分为曲率半径r = 50 mm的半圆，
Ѧ஝ѹथ͖ӑᤥᮊ ѣ԰౎͈
材料为钢，为节省计算时间，只截取 36 圆心角对应
◦
图 3 方程组求解算法流程图 弧长作为耦合介质与钢的交界面。分别对耦合介质
Fig. 3 Flow of equation system solving algorithm 为水和有机玻璃楔块进行验证，材料参数如表 1 所