Page 77 - 《应用声学》2021年第2期

第 40 卷 第 2 期 张玉等： 乳腺超声断层成像系统的有限元仿真分析 245

1.2 边界条件 80
60
波动方程是用来描述声学、光学、电磁学的波
40
动情况的一类重要的偏微分方程，用 COMSOL 中 20
“压力声学，瞬态”接口求解瞬态波动方程 [5] ： Y/mm 0
2 [ ] -20
1 ∂ p t 1
+ ∇ · − (∇p t − q d ) = Q m , (1) -40
2
ρc ∂t 2 ρ
-60
p t = p + p b , (2)
-80
其中，t表示时间，ρ 是流体密度，c是声速，p t 是总声 -80 -40 0 40 80
X/mm
压，q d 是偶极源，Q m 是单极源，p b 是背景压力波。辐
图 3 模型网格剖分图
射边界选择平面波辐射，完美匹配层中的典型波速
Fig. 3 Meshing of model
设为水中的声速即 1500 m/s，声压级参考压力使用
水的参考压力，即P ref,SPL = 1 µPa。
-65.0
固体域和水域的双向自动耦合在水域和换能 -66.0
器边界处实现，边界条件为 -67.0
[ ] -68.0
1
− n · − (∇p t − q d ) = −n · u tt , (3) Y/mm -69.0
ρc
-70.0
F A = p t n, (4)
-71.0
-72.0
其中，n 代表表面法线方向，u tt 是结构加速度，F A
-73.0
是作用在结构上的载荷(每单位面积的力)。 -32 -30 -28 -26 -24 -22 -20
换能器的压电效应在 PZT-5域中实现，通过求 X/mm
解线性方程，使得固体力学方程和静电方程相耦合， 图 4 局部网格剖分图
方程通过将应力和应变与电场和电位移相耦合来 Fig. 4 Local meshing
模拟压电效应。
solver, MUMPS)进行求解计算，该求解器具有完整
1.3 网格划分 的线性系统矩阵，能够实现实数和复数的运算，将求
在有限元仿真中，需要对所建模型进行网格划 解器中解的最大频率设置为f max,sol = 3 MHz。
分，COMSOL Multiphysics中提供交互式网格划分 确定时间步进的方法有向后差分公式、广义 α。
环境，每一个网格划分操作都将添加到网格划分序 在瞬态计算中，通常使用广义 α 的方法。因为向后
列中，当网格划分序列中所有操作完成以后，得到最 差分公式算法会产生散射，波形畸变随计算时间成
终的网格。时间步长以及剖分网格的大小和形状很 正比；广义 α 可以有效规避这种干扰，它用的是前 5
大程度上限制了数值计算的精度。在流体域中，每 个时间步长的解，并且可以对下个时间步长的解进
个波长使用 5∼6 个网格单元，通常设置最大单元尺 行预测 [6] 。本次计算采用广义 α 的方法，手动控制
寸为波长/5，即网格尺寸 = 波长/5。因为此模型中 时间步长，采用默认的时间步长即周期/60，这样的
水、软组织和铁块的波长不同，最大单元网格大小分 计算结果精确高，能达到最佳性能。
别为 0.1 mm、0.101 mm、0.345 mm。完成网格剖分 仿真完成后，对仿真结果进行后处理，得到在
之后网格顶点数为 26738，三角形单元数为 48102， 乳腺超声断层成像中，超声波在乳腺中的声压分布
边单元数为 10242，顶点单元数为1796，整个模型网
情况。图5 是第一个换能器发射后，时间为 27.8 µs、
格划分和局部划分如图3和图4所示。
41.7 µs、69.4 µs、97.2 µs时的分布云图。软组织、铁
块的声速和密度不同，其声阻抗系数也不同，声速、
2 计算结果与分析
密度越大，声阻抗越大 [7] ，超声换能器接收到的信
在网格剖分完成之后，选择并行稀疏直接 号越小。从图 5 中可以看出，超声波在水中逐渐衰
求解器(Multifrontal massive parallel sparse direct 减，当超声波遇到铁块时，能量有明显的衰减。