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
   81   82   83   84   85   86   87   88   89   90   91