Page 86 - 《应用声学》2021年第1期
P. 86
82 2021 年 1 月
∂ n+ 1 为 ξ = x、ξ = x, y、ξ = x, y, z,F 和 F 分别为傅里
−
n 2
ρ − ∆tρ 0 u
ξ ∂ξ ξ
ρ n+1 = , 叶变换和傅里叶逆变换,i 为虚数单位,k ξ 表示在 ξ
ξ ∂ n+ 1
1 + 2∆t u 2 方向的波数,∆ξ 表示在 ξ 方向的网格间距,∆t 是时
∂ξ ξ
( ) 2 ) 间步长,κ = sinc (c ref k∆t/2)为k-space算子,c ref 是
p n+1 = c 2 ρ n+1 + B 1 ( ρ n+1 − L d , (3) n+1 ∑ n+1 ,所有公式中的
0
2A ρ 0 参考声速,总密度 ρ = ρ ξ
ξ
密度 ρ n+1 被分解成笛卡尔分量以引入各向异性 上标 n 和n + 1 表示函数在当前和下一时间点的值,
ξ
完全匹配层 (Perfect matched layer, PML) [19] 。公 n − 1/2和n + 1/2表示时间交错点的值,k ξ 和L d 离
式 (3) 中,ξ 在一维空间、二维空间、三维空间分别 散形式如下:
[ ]
N ξ N ξ N ξ 2π
− , − + 1, · · · , − 1 , N ξ 为偶数,
2 2 2 ∆ξN ξ
k ξ = [ ]
(N ξ − 1) (N ξ − 1) (N ξ − 1) 2π
− , − + 1, · · · , , N ξ 为奇数,
2 2 2 ∆ξN ξ
{ { n }}
∂ρ { { }}
L d = τF −1 k y 1 −2 F + ηF −1 k y 1 −1 F ρ n+1 ,
∂t
式中的N ξ 是在ξ 方向的网格点数。 虽然公式 (4) 和公式 (5) 是无黏性流体情况得
方程 (3) 中的离散方程使用时间步长 ∆t = 到的,严格地说只适用于在不存在剪切且剪切模量
CFL∆x/c max 迭代求解,为了能够平衡准确性和 为零的流体中的传播,但它考虑了由黏性衰减和不
计算效率之间的关系,CFL 6 c 0 /c max 。在弱异构介 均匀性可能造成的散射而导致的波动量在封闭面
质,通常 CFL = 0.3 能够为精度和计算速度之间提 内体积沉积,而且在医学超声应用的实际环境中,
供很好的平衡 [19] 。在每个时间步长,可以通过在计 软组织具有相当低的吸收和散射特性,组织中的体
算域内的适当网格点添加源值来引入质量或力源。 模量远大于剪切模量,压缩应力比剪切应力更容易
类似地,模拟的输出可以通过在特定网格点的每个 出现,因此,公式 (4) 和公式 (5) 计算腹壁组织的声
时间步长记录声学变量来获得。 辐射力是合理的,对本文的问题仍然适用。基于公
声辐射力的一般表达式为 [14−15,20−21] 式 (5) 计算的声辐射力所用的声压和质点振速是由
∫∫
1 ⟨ ⟩ 1 ⟨ ⟩ k-Wave 运算过程中记录的声压和非交错网格中的
2
v
F = − 2 p 1 n + ρ 0 2 n
1
S 2ρ 0 c 0 2 质点振速。
⟨ T ⟩
− ρ 0 v 1 v · n dS, (4)
1
式(4)用到扰动分析以及 2 仿真模型的建立
ρ = ρ 0 + ρ 1 + ρ 2 + · · · , 本文的声场是由 100 个阵元组成的线性相控
v = v 1 + v 2 + · · · , 阵换能器, 如图 1(a) 所示; 其相应的物理参数
示意图如图 1(b) 所示。仿真中设置计算区域为
p = p 0 + p 1 + p 2 + · · · ,
(N x × N y × N z ),沿 x 方向的网格数 N x = 88,沿
其中,ρ、p、v 分别表示介质的密度、声压、质点速度,
y 方向的网格数 N y = 108, 沿 z 方向的网格数
n为面S 的法向矢量,⟨ ⟩表示时间平均,下标 0、1、2
N z = 44,阵元宽度为 1 个网格,阵元长度为 12
分别代表静态、一阶、二阶微小量。
个网格,阵元中心距离为 0,换能器的中心频率为
3
对于体积元(N/m )的声辐射力,公式(4)可以
1 MHz,换能器放置在yz 平面中间位置,声波沿x方
表示为
向传播,其聚焦点在 x 方向距离换能器 N x /2 × dx。
( )
∂ 1 ⟨ ⟩ 1 ⟨ ⟩
2
v
F i = − 2 p 1 − ρ 0 2 为了能够保证较高的精度和计算速度取网格大小
1
∂x i 2ρ 0 c 2
0 −4
dx = dy = dz = 1.5 × 10 m 和 CFL = 0.1。生
∂ ⟨v i v j ⟩
− ρ 0 . (5) 物组织样品 (腹壁组织) 放置在声场中,其横截面
∂x j