摘要:利用部分波展开法求解得到了Gauss声束入射下刚性和非刚性椭圆柱的声散射系数, 推导了一般情况下的声辐射力矩表达式. 在此基础上, 通过一系列数值仿真详细分析了离轴距离、入射角度和束腰半径对声辐射力矩的影响. 结果表明: 正向与负向声辐射力矩均可以在一定条件下存在; 低频情况下刚性椭圆柱比非刚性椭圆柱更容易产生较强的声辐射力矩; 特定频率的入射声场可以激发出非刚性椭圆柱不同阶的共振散射模式, 因而非刚性椭圆柱的声辐射力矩峰值与频率的关系更密切; 增加束腰半径有利于扩大散射截面, 进而增加椭圆柱的声辐射力矩. 该研究结果预期可以为利用声辐射力矩实现粒子的可控旋转和流体黏度的反演提供一定的理论指导.
关键词: 声辐射力矩/
Gauss声束/
椭圆柱/
部分波展开法/
声操控

English Abstract

全文HTML

--> --> -->

当声波入射到物体上时会对其产生力矩的作用, 称为声辐射力矩. 利用声辐射力矩可以实现粒子的可控旋转, 在声学操控^[1-4]、声悬浮^[5,6]、生物工程^[7]和材料科学^[8]等领域具有广泛的应用前景. 声辐射力矩的基本理论最早由Maidanik^[9]于1958年提出. 随后, Fan等^[10]从角动量守恒定律出发, 给出了一般情况下声辐射力矩的计算方法. 在实际应用中, 为了实现对靶向粒子的精确操控, 往往需要对各种情况下的声辐射力矩进行理论与数值计算. 在此基础上, Zhang和Marston^[11,12]系统地总结了声辐射力矩的理论框架, 并对悬浮固体和液滴的声辐射力矩进行了分析. 2012年, Silva和Mitri等^[13-15]成功推导出任意涡旋声束对在轴球形粒子声辐射力矩的表达式, 并以Bessel声束为例进行了仿真. 近年来, Zhang^[16]和Gong等^[17]详细分析了Bessel声束的声辐射力矩特性和力矩方向的决定因素, 其结果对分析涡旋声场中粒子的转动特性具有重要意义. 除了球形粒子外, 柱形粒子在实际应用中也很常见, 如纳米管^[18]、长纤维^[19]和声表面波器件^[20]等. 2007年, Hasheminejad^[21]利用Mathieu函数计算得到了无限长弹性椭圆柱在平面波场中的声辐射力和力矩解析式. 2011年, Wang和Dual^[22]从理论上给出了任意固体圆柱在黏性流体中受到的声辐射力矩大小并得到了很好的数值验证. 2016年, Mitri^[23]利用部分波级数展开法(PWSE)避开了Mathieu函数的烦琐运算, 成功求解得到了平面波以任意角度入射时刚性椭圆柱受到的声辐射力矩. 随后, Mitri^[24]又将研究范围拓展到Gauss聚焦声束, 对黏性流体中圆柱体受到的声辐射力与力矩进行了详细分析, 特别关注了力矩随离轴距离和入射角度的变化规律. 2017年, Mitri^[25]分析了一对黏性柱体的声辐射力矩并考虑了柱体间多次散射对结果的影响, 而后Wang等^[26]考虑了Airy波入射下的同样问题. 为了更精确地模拟实际边界附近的声操控问题, Mitri^[27]在已有模型中引入刚性边界, 利用加法公式推导了此时无限长圆柱的声辐射力矩.
本文在前人工作的基础上研究Gauss声束对理想流体中偏离声束中心的无限长椭圆柱的声辐射力矩, 采用部分波展开法分别对椭圆柱的散射系数进行求解, 进而计算此时的声辐射力矩. 在理论计算的基础上, 详细分析了离轴距离、入射角度和束腰半径等因素对声辐射力矩的影响, 为实现粒子在声场中的可控旋转提供了一定的理论基础. 值得注意的是, 本文模型中Gauss声束的入射方向可以是任意的, 因而具有良好的适用性. 此外, 为了不失一般性, 本文同时考虑了刚性和非刚性椭圆柱两种情况, 并对它们的声辐射力矩特性进行了详细的分析与比较.

2

-->

2.1.理论模型

如图1所示, 一束束腰半径为W₀的Gauss声束入射到理想流体中的无限长椭圆柱上. 以椭圆柱的中心为原点建立空间直角坐标系(x, y, z), 其中z轴垂直于xOy平面且并未在图中画出. 椭圆柱在x方向和y方向的半轴长度分别为a和b. 显然, 当a = b时, 椭圆柱退化为标准的圆柱体. 声束的传播方向与x轴正方向的夹角为α, 声束中心的空间坐标为(x₀, y₀). 为了后面数学讨论的方便, 以O为原点再建立一柱坐标系(r, θ, z), 两坐标系之间满足关系: $x = r\cos \theta,\; y = r\sin \theta,\; z = z$.
图 1 Gauss声束斜入射到无限长椭圆柱上
Figure1. Description for the interaction of a Gauss beam with an infinitely long elliptical cylinder.

2 -->

2.2.声场的级数展开

理想流体中Gauss声束的速度势函数可以表示为级数展开的形式^[28]:

$ {\phi }_{\rm{inc}}={\phi }_{0}{\rm{e}}^{-{\rm{i}}\omega t}{\displaystyle \sum _{n=-\infty }^{+\infty }{b}_{n}{\rm{J}}_{n}\left(kr\right){\rm{e}}^{{\rm{i}}n\left(\theta -\alpha \right)}}, $

其中, ${\phi _0}$表示声波速度势函数的振幅, k = ω/c₀表示声波的波数, ω是声波的角频率, c₀是流体中的纵波声速, J_n表示n阶柱Bessel函数, b_n是Gauss声束的波束形成因子. 当声束中心位于原点O时, 波束形成因子的表达式为^[28]

$ {b}_{n}={{\rm{i}}}^{n}\frac{{\rm{I}}_{n}\left(k{x}_{\rm{R}}\right)}{{\rm{I}}_{0}\left(k{x}_{0}\right)}, $

其中, I_n表示n阶第一类虚宗量柱Bessel函数, ${x_{\rm{R}}} = kW_0^2/2$表示Rayleigh距离. 当kW₀较小时, 声束具有较强的聚焦特性. 随着kW₀的增大, 声束逐渐扩散. 当kW₀趋向于无穷大时, Gauss声束将退化为平面行波. 当声束中心偏离原点O时, 可以通过柱Bessel函数的加法定理求解得到此时的波束形成因子, 具体表达式为

$ {b}_{n}={\displaystyle \sum _{p=-\infty }^{+\infty }{b}_{p}{\rm{J}}_{n-p}\left(kd\right){\rm{e}}^{-{\rm{i}}\left(n-p\right){\varphi }_{d}}}, $

其中, b_p是声束中心位于原点O时的波束形成因子, $d = \sqrt {x_0^2 + y_0^2} $是声束中心到原点O的距离, ${\phi _d} = $$ {\tan ^{ - 1}}\left[ {\left( {y - {y_0}} \right)/\left( {x - {x_0}} \right)} \right]$是偏移角度的大小.
类似地, 椭圆柱体的散射声场也可以表示为级数展开的形式:

$ {\phi }_{\rm{sca}}={\phi }_{0}{\rm{e}}^{-{\rm{i}}\omega t}{\displaystyle \sum _{n=-\infty }^{+\infty }{b}_{n}{s}_{n}{\rm{H}}_{n}^{\left(1\right)}\left(kr\right){\rm{e}}^{{\rm{i}}n\left(\theta -\alpha \right)}}, $

其中, ${\rm{H}}_n^{\left( 1 \right)}$表示n阶第一类柱Hankel函数, s_n是声波散射系数, 将根据椭圆柱表面特定的边界条件来决定.
2 -->

2.3.散射系数的求解

定义椭圆柱的形状函数为^[23]

${S_\theta } = {\left( {{{\left( {\frac{{\cos \theta }}{a}} \right)}^2} + {{\left( {\frac{{\sin \theta }}{b}} \right)}^2}} \right)^{ - \frac{1}{2}}}.$

显然对于圆柱(a = b)而言, 形状函数是一个固定的常数; 对于一般的椭圆柱而言, 形状函数是角度θ的函数.
当椭圆柱的声阻抗远大于周围流体的声阻抗时, 椭圆柱可以近似看作是刚性的, 此时椭圆柱表面需要满足法向速度为零的Neumann边界条件, 即

$ {\nabla \left({\phi }_{\rm{inc}}+{\phi }_{\rm{sca}}\right)\cdot {{n}}|}_{r={S}_{\theta }}=0, $

其中, ${{n}} = {{{e}}_r} - \dfrac{1}{{{S_\theta }}}\dfrac{{{{\rm{d}}}{S_\theta }}}{{{{\rm{d}}}\theta }}{{{e}}_\theta }$表示法向单位矢量, e_r和e_θ分别表示径向和角度方向的单位矢量. 将(1)式和(4)式代入(6)式, 可以得到关于散射系数s_n的方程:

$\sum\limits_{n = - \infty }^{ + \infty } {{b_n}} {{\rm{e}}^{ - {{\rm{i}}}n\alpha }}\left[ {{A_n}\left( \theta \right) + {s_n}{B_n}\left( \theta \right)} \right] = 0,$

其中, A_n和B_n的具体表达式为

$\begin{split}\left\{\begin{array}{l}{A}_{n}\\ {B}_{n}\end{array}\right\}=\;&k{\rm{e}}^{{\rm{i}}n\theta }\left\{\begin{array}{l}{\rm{J}}_{n}'\left(k{S}_{\theta }\right)\\ {{\rm{H}}_{n}^{\left(1\right)}}'\left(k{S}_{\theta }\right)\end{array}\right\}\\&-{\rm{i}}\frac{n}{{S}_{\theta }^{2}}\frac{{\rm{d}}{S}_{\theta }}{{\rm{d}}\theta }{\rm{e}}^{{\rm{i}}n\theta }\left\{\begin{array}{l}{\rm{J}}_{n}\left(k{S}_{\theta }\right)\\ {\rm{H}}_{n}^{\left(1\right)}\left(k{S}_{\theta }\right)\end{array}\right\}, \end{split}$

对于标准的圆柱体而言, (8)式右边的第二项为零, 因此(7)式的左边正是以${{\rm{e}}^{{{\rm{i}}}n\theta }}$为基函数的Fourier级数展开形式, 级数恒为零等价于级数的每一项恒为零, 从而可以直接得到散射系数的表达式$ {s_n} = - {{\rm{J}}_n'}\left( {ka} \right)/{{\rm{H}}_n^{\left( 1 \right)}}'\left( {ka} \right) $. 然而对于一般情况下的椭圆柱, (7)式的左边不再是Fourier级数, 无法通过直接去掉求和号来计算散射系数. 尽管如此, 我们形式上仍然可以将(7)式的左边改写为Fourier级数的形式

$ \begin{split}&{\displaystyle \sum _{n=-\infty }^{+\infty }{b}_{n}}{\rm{e}}^{-{\rm{i}}n\alpha }\left[{A}_{n}\left(\theta \right)+{s}_{n}{B}_{n}\left(\theta \right)\right]\\&={\displaystyle \sum _{n=-\infty }^{+\infty }\left({\overline{A}_{n}}+{s}_{n}{\overline{B}_{n}}\right){\rm{e}}^{{\rm{i}}n\theta }}{=0},\end{split}$

其中, $\overline A_n$和$\overline B_n$是此时的Fourier展开系数, 可以通过Fourier逆变换计算, 即

$\left\{ \begin{array}{l} \overline {{A_l}} \\ \overline {{B_l}} \end{array} \right\} = \frac{{\rm{1}}}{{{\rm{2\pi }}}}\sum\limits_n {{b_n}{{\rm{e}}^{ - {{\rm{i}}}n\alpha }}} \int_0^{{\rm{2\pi }}} {\left\{ \begin{array}{l} {A_n}\left( \theta \right) \\ {B_n}\left( \theta \right) \end{array} \right\}} {{\rm{e}}^{ - {{\rm{i}}}l\theta }}{{\rm{d}}}\theta .$

此时, (9)式正是Fourier级数展开式, 可以直接求解得到散射系数s_n.
当椭圆柱的声阻抗与周围流体相差不大时, 不能再将椭圆柱看作刚性, 即此时必须考虑椭圆柱内部透射声波的存在. 透射声场的速度势函数可以表示为

$ {\phi }_{\rm{l}}={\phi }_{0}{\rm{e}}^{-{\rm{i}}\omega t}{\displaystyle \sum _{n=-\infty }^{+\infty }{b}_{n}{t}_{n}{\rm{H}}_{n}^{\left(1\right)}\left(k_{1}r\right){\rm{e}}^{{\rm{i}}n\left(\theta -\alpha \right)}}, $

其中, k_l = ω/c_l和c_l分别是椭圆柱内部的波数和纵波声速, t_n是声波透射系数. 对于这样的非刚性椭圆柱, 边界条件要求声压和法向速度在表面连续,

$\left\{ \begin{aligned} &{\left. {\left( {{p_{{\rm{inc}}}} + {p_{{\rm{sca}}}}} \right)} \right|_{r = {S_\theta }}} = {\left. {{p_{\rm{l}}}} \right|_{r = {S_\theta }}}, \\ &{\left. {\left( {{{{v}}_{{\rm{inc}}}} + {{{v}}_{{\rm{sca}}}}} \right) \cdot {{n}}} \right|_{r = {S_\theta }}} = {\left. {{{{v}}_{\rm{l}}} \cdot {{n}}} \right|_{r = {S_\theta }}}. \end{aligned} \right.$

假设周围流体和椭圆柱的密度分别为ρ₀和ρ_l, (12)式中的各项声压和速度可以通过它们与速度势的关系求出, 即${p_{{{\rm{inc}}, {\rm{sca}}}}} = - {{\rm{i}}}\omega {\rho _0}{\phi _{{{\rm{inc}}, {\rm{sca}}}}}, \;{p_{{\rm{l}}}} = $$ - {{\rm{i}}}\omega {\rho _{\rm{l}}}{\phi _{{\rm{l}}}}$和${{{v}}_{{{\rm{inc}}, {\rm{sca}}}}} = - \nabla {\phi _{{{\rm{inc}}, {\rm{sca}}}}}, \;{{{v}}_{{\rm{l}}}} = - \nabla {\phi _{{\rm{l}}}}$. 将 (1)式, (4)式和(11)式代入(12)式可以得到如下的方程组:

$\left\{ \begin{aligned} & \quad - {{\rm{i}}}\omega {\rho _0}\sum\limits_{n = - \infty }^{ + \infty } {{b_n}} {{\rm{e}}^{ - {{\rm{i}}}n\alpha }}\left[ {{D_n}\left( \theta \right) + {s_n}{E_n}\left( \theta \right)} \right] \\ &= - {{\rm{i}}}\omega {\rho _{\rm{l}}}\sum\limits_{n = - \infty }^{ + \infty } {{b_n}} {{\rm{e}}^{ - {{\rm{i}}}n\alpha }}{t_n}{F_n}\left( \theta \right), \\ &\quad \sum\limits_{n = - \infty }^{ + \infty } {{b_n}} {{\rm{e}}^{ - {{\rm{i}}}n\alpha }}\left[ {{G_n}\left( \theta \right) + {s_n}{{K}_n}\left( \theta \right)} \right] \\ &= \sum\limits_{n = - \infty }^{ + \infty } {{b_n}} {{\rm{e}}^{ - {{\rm{i}}}n\alpha }}{t_n}{L_n}\left( \theta \right), \end{aligned} \right.$

其中

$ \begin{split}&\left\{\begin{array}{l}{D}_{n}\\ {E}_{n}\\ {F}_{n}\end{array}\right\}={\rm{e}}^{{\rm{i}}n\theta }\left\{\begin{array}{l}{\rm{J}}_{n}\left(k{S}_{\theta }\right)\\ {\rm{H}}_{n}^{\left(1\right)}\left(k{S}_{\theta }\right)\\ {\rm{J}}_{n}\left({k}_{\rm{l}}{S}_{\theta }\right)\end{array}\right\}, \\ &\left\{\begin{array}{l}{G}_{n}\\ {K}_{n}\end{array}\right\}=k{\rm{e}}^{{\rm{i}}n\theta }\left\{\begin{array}{l}{\rm{J}}_{n}'\left(k{S}_{\theta }\right)\\ {{\rm{H}}_{n}^{\left(1\right)}}'\left(k{S}_{\theta }\right)\end{array}\right\}\\&\qquad\quad\qquad-{\rm{i}}\frac{n}{{S}_{\theta }{}^{2}}\frac{{\rm{d}}{S}_{\theta }}{{\rm{d}}\theta }{\rm{e}}^{{\rm{i}}n\theta }\left\{\begin{array}{l}{\rm{J}}_{n}\left(k{S}_{\theta }\right)\\ {\rm{H}}_{n}^{\left(1\right)}\left(k{S}_{\theta }\right)\end{array}\right\},\\ &L_{n}={k}_{\rm{l}}{\rm{e}}^{{\rm{i}}n\theta }{\rm{J}}_{n}'\left({k}_{\rm{l}}{S}_{\theta }\right)-{\rm{i}}\frac{n}{{S}_{\theta }{}^{2}}\frac{{\rm{d}}{S}_{\theta }}{{\rm{d}}\theta }{\rm{e}}^{{\rm{i}}n\theta }{\rm{J}}_{n}\left({k}_{\rm{l}}{S}_{\theta }\right), \end{split}$

仿照前述对于刚性椭圆柱的处理, 同样将(13)式改写为Fourier级数的形式,

$\left\{ \begin{aligned} & \quad - {{\rm{i}}}\omega {\rho _0}\sum\limits_{n = - \infty }^{ + \infty } {\left( {\overline {{D_n}} + {s_n}\overline {{E_n}} } \right){{\rm{e}}^{{{\rm{i}}}n\theta }}}\\ &= - {{\rm{i}}}\omega {\rho _{\rm{l}}}\sum\limits_{n = - \infty }^{ + \infty } {{t_n}\overline {{F_n}} {{\rm{e}}^{{{\rm{i}}}n\theta }}}, \\ &\sum\limits_{n = - \infty }^{ + \infty } {\left[ {\overline {{G_n}} + {s_n}\overline {{K_n}} } \right]{{\rm{e}}^{{{\rm{i}}}n\theta }}} = \sum\limits_{n = - \infty }^{ + \infty } {{t_n}\overline {{L_n}} {{\rm{e}}^{{{\rm{i}}}n\theta }}}, \end{aligned} \right.$

其中

$\left\{ \begin{array}{l} \overline {{D_{\rm{l}}}} \\ \overline {{E_{\rm{l}}}} \\ \overline {{F_{\rm{l}}}} \\ \overline {{G_{\rm{l}}}} \\ \overline {{K_{\rm{l}}}} \\ \overline {{L_{\rm{l}}}} \end{array} \right\} = \frac{{\rm{1}}}{{{\rm{2\pi }}}}\sum\limits_n {{b_n}} {{\rm{e}}^{ - {{\rm{i}}}n\alpha }}\int_0^{2{\rm{\pi }}} {\left\{ \begin{array}{l} {D_n}\left( \theta \right) \\ {E_n}\left( \theta \right) \\ {F_n}\left( \theta \right) \\ {G_n}\left( \theta \right) \\ {{{K}}_n}\left( \theta \right) \\ {{{L}}_n}\left( \theta \right) \end{array} \right\}} {{\rm{e}}^{ - {{\rm{i}}}l\theta }}{{\rm{d}}}\theta .$

此时, 通过求解(15)式可以得到散射系数s_n.
2 -->

2.4.声辐射力矩的计算

声辐射力矩的基本计算公式最早在文献[9]中给出, 其表达式为

$ {\langle {{N}}\rangle }_{kr\to \infty }=-{\rho }_{0}{\displaystyle {\iint }_{S}{{v}}\langle {{r}}\times {{v}}\rangle {\rm{d}}{{S}}}, $

其中, 符号$\left\langle {} \right\rangle $表示对时间平均, r是以O为起点的径矢, ${{\rm{d}}}{{S}} = {\rm{d}}S{{{e}}_r}$, dS = Lrdθ是将椭圆柱包含在内的一个半径为r的大积分圆柱面上长度为L的面积微元, ${{{e}}_r} = \cos \theta {{{e}}_x} + \sin \theta {{{e}}_y}$是积分表面的外法向单位矢量, ${{{e}}_x}$和${{{e}}_y}$是xOy坐标系沿坐标轴方向的单位矢量. 考虑到周围的流体是理想的, 不存在对声波动量或角动量的吸收, 因此可以用这一大的积分柱面代替原来的椭圆柱面从而利用柱函数的远场近似来简化我们的运算.
在本文的模型中, 只有z方向的声辐射力矩不为零, 利用速度和速度势函数之间的关系及两个复数量乘积的运算性质, 可以将z方向的力矩表示为

$\left\langle {{N_z}} \right\rangle = \left\langle {{{{N}}_z}} \right\rangle \cdot {{{e}}_z} = \frac{{{\rho _0}}}{2}{\rm{Im}} \left\{ {\iint_S {\frac{{\partial {\phi^* _{{\rm{total}}}}}}{{\partial r}} {{\hat L_z}} {\phi _{{\rm{total}}}}{{\rm{d}}}{{S}}}} \right\},$

其中, ${\phi _{{\rm{total}}}} = {\phi _{{\rm{inc}}}} + {\phi _{{\rm{sca}}}}$表示流体中的总声场, ${{\hat L_z}} = - {{\rm{i}}}\dfrac{\partial }{{\partial \theta }}$表示z方向的角动量算符, Im表示对复数量取虚部. 在远场近似下, (1)式和(4)式所表示的入射和散射声场可以表示为

$ \begin{split}&{{\phi }_{\rm{inc}}|}_{kr\to \infty }\\=\;&{\phi }_{0}{\rm{e}}^{-{\rm{i}}\omega t}\!\sqrt{\frac{2}{{\rm{\pi}}kr}}{\displaystyle\! \sum _{n=-\infty }^{+\infty }\!{b}_{n}{\rm{cos}}\!\left(kr\!-\!\frac{2n+1}{4}{\rm{\pi}}\!\right)\!{\rm{e}}^{{\rm{i}}n\left(\theta -\alpha \right)}}, \end{split}$

$ {{\phi }_{\rm{sca}}|}_{kr\to \infty }\!={\phi }_{0}{\rm{e}}^{{\rm{i}}\left(kr-\omega t\right)}\!\sqrt{\frac{2}{{\rm{i \pi}}kr}\!}{\displaystyle \sum _{n=-\infty }^{+\infty }\!{b}_{n}{s}_{n}{{\rm{i}}}^{-n}{\rm{e}}^{{\rm{i}}n\left(\theta -\alpha \right)}}, $

将(19)式和(20)式代入(18)式中, 并利用如下的关系

$ {\displaystyle {\int }_{\rm{0}}^{\rm{2 \pi}}{\rm{e}}^{{\rm{i}}\left(n'-n\right)\theta }{\rm{d}}\theta }=2{\rm{\pi}}{\delta }_{n,n'}, $

其中, $ {\delta _i}_{, j}$是Kronecker符号. 经过一系列数学运算, 最终可以得到长度为L的椭圆柱受到的声辐射力矩为^[24]

$ \langle {N}_{z}\rangle ={\rm{\pi}}{b}^{2}LE{\tau }_{z}, $

其中, $E = {\rho _{\rm{0}}}{k^2}\phi _{\rm{0}}^{\rm{2}}/{\rm{2}}$是声波的平均空间能量密度, τ_z称为声辐射力矩函数, 其具体表达式为

${\tau _z} = - \frac{4}{{{\rm{\pi }}{{\left( {kb} \right)}^2}}}\sum\limits_{n = - \infty }^{ + \infty } {n{{\left| {{b_n}} \right|}^2}\left[ {{\rm{Re}} \left( {{s_n}} \right) + {{\left| {{s_n}} \right|}^2}} \right]},$

其中, Re表示对复数量取实部.
进一步考察声辐射力矩函数(23)式, 不难发现单极辐射项(n = 0)对声辐射力矩没有贡献. 在不考虑柱体对声波黏滞吸收的情况下, 对于标准的圆柱体而言${\rm{Re}} \left( {{s_n}} \right) + {\left| {{s_n}} \right|^2}$恒为零, 即此时不存在声辐射力矩. 这一现象在关于声辐射力矩的其他文献^[23,24]中也多次提及. 但需要注意的是, 对于椭圆柱而言, 即使不考虑声波的黏滞吸收, ${\rm{Re}} \left( {{s_n}} \right) + {\left| {{s_n}} \right|^2}$也不一定为零, 即声辐射力矩是可以存在的.

根据(23)式可以对椭圆柱的声辐射力矩进行数值仿真. 在计算过程中, 以上的所有无穷级数必须进行截断. 若截断的阶数太低, 则级数收敛的精度不够, 进而导致结果不准确; 若截断的阶数过高, 则会增加不必要的计算成本. 因此, 实际中必须予以统筹兼顾, 在合适的阶数将无穷级数进行截断处理并确保结果收敛. 文献[23]中对这一问题已经进行了详细讨论, 证明了在$\left| {{s_{n + {n_{\max }}}}/{s_0}} \right| \sim {10^{ - {\rm{16}}}}$处截断时截断误差可以忽略, 因此本文所有的计算中也采用了这一截断准则. 根据文献[23]中的分析, 随着长短轴之比的增加, 椭圆柱偏离圆柱越明显, 所需要的截断阶数也越高. 在本文的所有仿真实例中, 长短轴之比最大为2, 此时需要的截断阶数仍小于30, 可以提供足够的运算效率. 有必要指出, 在整个计算过程中, (10)式和(16)式涉及的对角度θ的定积分需要采用数值积分方法进行求解, 这里我们选用Simpson方法. 在本文的所有仿真实例中, 均假设椭圆柱位于水中, 水的密度为1000 kg/m³, 纵波声速为1480 m/s. 此外, 对于非刚性椭圆柱的情形, 我们将生物医学超声中常见的两种材料聚二甲基硅氧烷(PDMS)和四溴乙烷(TBE)混合成密度刚好为1000 kg/m³的液体(PDMS-TBE)并以此为例进行仿真, 在常温下, 这两种物质均以液体形式存在, 合成得到的PDMS-TBE材料的纵波声速为930 m/s^[29].
2

-->

3.1.椭圆柱的位置对声辐射力矩的影响

图2给出了声束中心沿y轴上下移动(kx₀ = 0)时, a/b取不同值的情况下刚性椭圆柱和PDMS-TBE椭圆柱的声辐射力矩函数τ_z随kb和ky₀的变化关系, Gauss声束的束腰半径满足kW₀ = 3, 入射角度α = π/4. 可以看出, 由于声束是倾斜入射的, 所以位于声束中心的椭圆柱(ky₀ = 0)也会受到声辐射力矩的作用. 此外, 声辐射力矩可以为正或为负. 根据力矩方向的定义法则, 当力矩为正值时, 椭圆柱将绕z轴作逆时针转动, 反之将作顺时针转动. 特别地, 对于标准的圆柱体(a/b = 1), 无论是刚性还是非刚性的情况下声辐射力矩均恒为零, 这与前面的理论分析是相符的. 对于刚性椭圆柱而言, 在0 < kb < 1的低频范围内声辐射力矩出现峰值, 并集中在ky₀ = 0附近. 当a/b < 1时峰值为正, 且随着kb的增加峰值逐渐向ky₀的负半轴移动, 而当a/b > 1时峰值为负, 且随着kb的增加逐渐向ky₀的正半轴移动. 无论是何种形状的刚性椭圆柱, 随着频率的增加, 声辐射力矩的峰值均出现明显的衰减. 对于PDMS-TBE椭圆柱而言, 声辐射力矩函数图像关于ky₀ = 0基本呈现奇对称的特征, 即当声束中心在关于y轴对称的位置力矩恰好大小相等而方向相反. 从图2中还可以看出, PDMS-TBE椭圆柱的声辐射力矩函数峰值形成了若干平行于ky₀轴的条带, 随着a/b的增加, 这些条带逐渐变细, 即力矩峰值对频率的依赖型更强. 刚性椭圆柱的形状偏离圆柱越明显, 力矩的峰值越大, 但PDMS-TBE椭圆柱却不满足这一规律.
图 2 椭圆柱的声辐射力矩函数随kb和ky₀的变化关系(kx₀ = 0, kW₀ = 3)　(a) a/b = 1/2, 刚性; (b) a/b = 1/2, PDMS-TBE; (c) a/b = 2/3, 刚性; (d) a/b = 2/3, PDMS-TBE; (e) a/b = 1, 刚性; (f) a/b = 1, PDMS-TBE; (g) a/b = 3/2, 刚性; (h) a/b = 3/2, PDMS-TBE; (i) a/b = 2, 刚性; (j) a/b = 2, PDMS-TBE
Figure2. The acoustic radiation torque function plots for an elliptical cylinder versus kb and ky₀ (kx₀ = 0, kW₀ = 3): (a) a/b = 1/2, rigid; (b) a/b = 1/2, PDMS-TBE; (c) a/b = 2/3, rigid; (d) a/b = 2/3, PDMS-TBE; (e) a/b = 1, rigid; (f) a/b = 1, PDMS-TBE; (g) a/b = 3/2, rigid; (h) a/b = 3/2, PDMS-TBE; (i) a/b = 2, rigid; (j) a/b = 2, PDMS-TBE.

为了进一步理解非刚性椭圆柱声辐射力矩函数仿真图中出现的条带, 对其远场散射声场进行进一步分析. 将散射声场的远场近似式(20)式改写为

$ {{\phi }_{\rm{sca}}|}_{kr\to \infty }={\phi }_{0}\sqrt{\frac{b}{2r}}{f}^{\rm{total}}, $

其中, ${f^{{\rm{total}}}}$是总的远场散射函数, 其具体表达式为

${f^{{\rm{total}}}} = {{\rm{e}}^{{{\rm{i}}}\left( {kr - \omega t} \right)}}\frac{2}{{\sqrt {{\rm{\pi }}kb} }}\sum\limits_{n = - \infty }^{ + \infty } {{b_n}{s_n}{{{\rm{i}}}^{ - n}}{{\rm{e}}^{{{\rm{i}}}n\left( {\theta - \alpha } \right)}}} .$

根据文献[30,31]的讨论, 为了得到单纯的远场共振散射函数, 需要从(25)式中减去刚性椭圆柱的散射函数, 即

${f^{{\rm{res}}}} \!=\! {{\rm{e}}^{{{\rm{i}}}\left( {kr - \omega t} \right)}}\!\frac{2}{{\sqrt {{\rm{\pi }}kb} }}\!\sum\limits_{n = - \infty }^{ + \infty }\! {{b_n}\!\left(\! {{s_n} \!-\! {s_n}^{{\rm{rigid}}}} \!\right)\!{{{\rm{i}}}^{ - n}}{{\rm{e}}^{{{\rm{i}}}n\left( {\theta - \alpha } \right)}}} .$

以a/b = 3/2的PDMS-TBE椭圆柱为例, 利用(26)式绘出当kb = 5.5与kb = 6.1时远场共振散射函数的幅值随角度的分布图像, 分别如图3(a)和图3(b)所示. 这两个kb值分别是图2(h)中的两条带对应的频率. 椭圆柱位于kx₀ = 0, ky₀ = 6处, Gauss声束入射角度α = π/4, 束腰半径kW₀ = 3. 从图3(a)可以看出, kb = 5.5时共振散射函数的指向性图案在两个相反的方向出现了极大值, 这恰好是偶极辐射(n = 1)的典型特征. 此时偶极辐射项占据主导地位, 而其余项可以忽略不计. 同样地, 从图3(b)可以看出, kb = 6.1时共振散射函数的指向性图案出现了4个极大值, 此时四极辐射项(n = 2)占据主导地位. 因此, 对于非刚性椭圆柱而言, 其声辐射力矩函数图像的峰值形成的条带正好对应着不同阶的共振散射模式. 由于单极模式(n = 0)对声辐射力矩没有贡献, 因而在低频区域PDMS-TBE椭圆柱受到的声辐射力矩很微弱, 这与刚性情况形成了鲜明的对比. 在操控非刚性椭圆柱的旋转时, 需要根据实际需要选取合适的频率激发出其自身的共振散射模式, 从而得到明显的声辐射力矩.
图 3 PDMS椭圆柱的共振散射函数幅值$\left| {{f^{{\rm{res}}}}} \right|$随角度θ的变化关系$(kx_0 \!=\! 0, ky_0 \!=\! 6, \alpha \!=\! \pi/4, kW_0 \!=\! 3) ~~{\rm(a)}~ kb \!=\! 5.5; ~{\rm (b)}~ kb \!=\! 6.1$
Figure3. The resonance scattering function modulus $\left| {{f^{{\rm{res}}}}} \right|$ for a PDMS-TBE elliptical cylinder versus the angle θ (kx₀ = 0, ky₀ = 6, α = π/4, kW₀ = 3): (a) kb = 5.5; (b) kb = 6.1.

图4给出了椭圆柱在沿平行于x轴的方向移动时声辐射力矩函数τ_z随kb和kx₀的变化关系, 其中ky₀ = –3, kW₀ = 3, α = π/4. 与图2的结果类似, 刚性椭圆柱的声辐射力矩在a/b < 1时主要为正, 而在a/b > 1时主要为负. 值得注意的是, 在低频范围内, 当声束中心位于y轴附近时力矩较小, 随着kx₀的增大, 力矩也随之增强. 然而在中高频范围内, 声辐射力矩却随着声束中心偏离y轴而逐渐衰减. 对于非刚性椭圆柱而言, 其声辐射力矩函数图像依然保持了关于kx₀ = 0奇对称的特性, 且其峰值形成了一系列平行于kx₀轴的条带. 对比图2和图4可以发现, 对于特定形状的PDMS-TBE椭圆柱, 这些条带所对应的kb值是完全相同的. 如前所述, 这些条带反映了不同阶散射项的共振频率, 因而只与椭圆柱的材料和形状有关. 同样地, a/b的增加会使仿真图中的条带变细, 即力矩峰值对频率的依赖性更强.
图 4 椭圆柱的声辐射力矩函数随kb和kx₀的变化关系(ky₀ = –3, kW₀ = 3)　(a) a/b = 1/2, 刚性; (b) a/b = 1/2, PDMS-TBE; (c) a/b = 2/3, 刚性; (d) a/b = 2/3, PDMS-TBE; (e) a/b = 1, 刚性; (f) a/b = 1, PDMS-TBE; (g) a/b = 3/2, 刚性; (h) a/b = 3/2, PDMS-TBE; (i) a/b = 2, 刚性; (j) a/b = 2, PDMS-TBE
Figure4. The acoustic radiation torque function plots for an elliptical cylinder versus kb and kx₀ (ky₀ = –3, kW₀ = 3): (a) a/b = 1/2, rigid; (b) a/b = 1/2, PDMS-TBE; (c) a/b = 2/3, rigid; (d) a/b = 2/3, PDMS-TBE; (e) a/b = 1, rigid; (f) a/b = 1, PDMS-TBE; (g) a/b = 3/2, rigid; (h) a/b = 3/2, PDMS-TBE; (i) a/b = 2, rigid; (j) a/b = 2, PDMS-TBE.

2 -->

3.2.声束入射角度对声辐射力矩的影响

声束的入射角度也是影响声辐射力矩的重要因素. 图5给出了当声束中心固定于kx₀ = ky₀ = –3的位置时在不同角度入射下的声辐射力矩函数仿真示意图, 束腰半径仍满足kW₀ = 3. 对于a/b < 1的刚性椭圆柱而言, 随着入射角度的增大, 在kb < 1的低频范围内声辐射力矩由正转负. 在kb > 1的中高频范围内, 小角度入射更容易产生较强的负向声辐射力矩. 随着角度的增加, 这一负向力矩逐渐减小并过渡到正向力矩. 对于a/b > 1的刚性椭圆柱, 结论恰好与a/b < 1的情况相反. 需要注意的是, 在α = 0°或α = 90°的情况下, 由于椭圆柱此时偏离Gauss声束的声轴, 声辐射力矩仍然不为零, 这与文献[23]中平面波入射下的结果很不相同. 当平面波垂直入射到椭圆柱上时, 声辐射力矩恒为零. 对于PDMS-TBE椭圆柱而言, 其声辐射力矩的一系列峰值源于不同阶的共振散射, 因而其声辐射力矩函数几乎与入射角度无关, 这有利于我们在一般情况下利用Gauss声束操控非刚性粒子的转动.
图 5 椭圆柱的声辐射力矩函数随kb和α的变化关系(kx₀ = –3, ky₀ = –3, kW₀ = 3)　(a) a/b = 1/2, 刚性; (b) a/b = 1/2, PDMS-TBE; (c) a/b = 2/3, 刚性; (d) a/b = 2/3, PDMS-TBE; (e) a/b = 1, 刚性; (f) a/b = 1, PDMS-TBE; (g) a/b = 3/2, 刚性; (h) a/b = 3/2, PDMS-TBE; (i) a/b = 2, 刚性; (j) a/b = 2, PDMS-TBE
Figure5. The acoustic radiation torque function plots for an elliptical cylinder versus kb and α (kx₀ = –3, ky₀ = –3, kW₀ = 3): (a) a/b = 1/2, rigid; (b) a/b = 1/2, PDMS-TBE; (c) a/b = 2/3, rigid; (d) a/b = 2/3, PDMS-TBE; (e) a/b = 1, rigid; (f) a/b = 1, PDMS-TBE; (g) a/b = 3/2, rigid; (h) a/b = 3/2, PDMS-TBE; (i) a/b = 2, rigid; (j) a/b = 2, PDMS-TBE.

2 -->

3.3.声束束腰半径对声辐射力矩的影响

束腰半径是描述Gauss声束的重要参量, 反映了声束聚焦特性的强弱. 图6显示了声辐射力矩函数τ_z随kb和kW₀变化的图像, 其中kW₀的变化范围为1 < kW₀ < 10. 仿真结果显示, 随着kW₀的增大, 刚性椭圆柱的声辐射力矩也随之增强. 这是由于束腰半径的增加使得椭圆柱散射截面增大, 从而提升了声辐射力矩大小. 在不同的频率范围内, 声辐射力矩增大的速率也有所不同. 在kb < 1的低频范围内, 椭圆柱声辐射力矩的增强最为明显, 而在中高频范围内束腰半径对声辐射力矩的影响较为微弱. 因此, 随着束腰半径的增加, 刚性椭圆柱的声辐射力矩随kb变化而出现的振荡也更加显著. 对PDMS-TBE椭圆柱而言, 随着波束变宽, 声辐射力矩的峰值也逐渐增强, 但峰值所对应的kb值并不改变. 值得一提的是, 随着kW₀的不断增大, 声束的聚焦特性越来越弱, 最终将退化为平面波入射的情形.
图 6 椭圆柱的声辐射力矩函数随kb和kW₀的变化关系(kx₀ = –3, ky₀ = –3, α = π/4)　(a) a/b = 1/2, 刚性; (b) a/b = 1/2, PDMS-TBE; (c) a/b = 2/3, 刚性; (d) a/b = 2/3, PDMS-TBE; (e) a/b = 1, 刚性; (f) a/b = 1, PDMS-TBE; (g) a/b = 3/2, 刚性; (h) a/b = 3/2, PDMS-TBE; (i) a/b = 2, 刚性; (j) a/b = 2, PDMS-TBE
Figure6. The acoustic radiation torque function plots for an elliptical cylinder versus kb and kW₀ (kx₀ = –3, ky₀ = –3, α = π/4): (a) a/b = 1/2, rigid; (b) a/b = 1/2, PDMS-TBE; (c) a/b = 2/3, rigid; (d) a/b = 2/3, PDMS-TBE; (e) a/b = 1, rigid; (f) a/b = 1, PDMS-TBE; (g) a/b = 3/2, rigid; (h) a/b = 3/2, PDMS-TBE; (i) a/b = 2, rigid; (j) a/b = 2, PDMS-TBE.

2 -->

3.4.声吸收对声辐射力矩的影响

在以上的所有计算中, 我们均忽略了椭圆柱的声吸收效应. 然而, 对于PDMS和TBE这样的高分子化合物而言, 声吸收效应往往比较显著, 可能对其声辐射力矩特性产生一定的影响. 为了分析声吸收对声辐射力矩的影响, 在具体的计算过程中只需要给波数增加一项虚部即可, 即此时椭圆柱中的声波波数为$ {{\tilde k_{\rm{l}}}} = {k_{\rm{l}}}\left( {1 - {{\rm{i}}}{\gamma _{\rm{l}}}} \right)$, 这里γ_l是椭圆柱的归一化纵波吸声系数. 图7同时给出了在考虑声吸收和不考虑声吸收的情况下a/b = 1/2的PDMS-TBE椭圆柱的声辐射力矩函数随kb的变化曲线, 其中kx₀ = 0, ky₀ = 6, kW₀ = 3, α = π/4, γ_l = 0.02. 不难发现, 两条曲线的变化趋势几乎完全一致, 声吸收对声辐射力矩的影响主要体现在kb > 6的高频范围. 当考虑声吸收时, PDMS-TBE椭圆柱的声辐射力矩在高频会略微减小, 但差异并不显著. 因此在考虑非刚性椭圆柱的声辐射力矩时, 忽略声吸收效应的影响是一种可以接受的合理近似.
图 7 PDMS-TBE椭圆柱的声辐射力矩函数随kb的变化关系(kx₀ = 0, ky₀ = 6, kW₀ = 3, α = π/4)
Figure7. The acoustic radiation torque function plots for a PDMS-TBE elliptical cylinder versus kb (kx₀ = 0, ky₀ = 6, kW₀ = 3, α = π/4).

2 -->

3.5.利用声辐射力矩进行流体黏度的反演

至此, 我们对理想流体中的椭圆柱在Gauss声束作用下的声辐射力矩问题进行了详细的参数化分析, 这些结果有助于在实际的操控中选取合适的声场和材料参数从而达到预期的操控目标. 对于Gauss声场中的椭圆柱而言, 声辐射力矩将作为驱动力矩使其产生绕z轴的加速转动. 此外, 在实际情况下, 周围流体往往是非理想的, 会对椭圆柱产生反向的黏滞力矩. 综上, 根据动量矩定理不难写出此时椭圆柱的转动方程:

${\rm{\pi }}{b^2}{E_0}{\tau _z} - {N_{{\rm{d}}}} = J\frac{{{{\rm{d}}}{\omega _z}}}{{{{\rm{d}}}t}}, $

其中${N_{{\rm{d}}}} = f\eta {\omega _z}$是单位长度椭圆柱受到的黏滞力矩^[32], 其与角速度大小ω_z和周围流体的黏度η成正比; f为比例系数; J为单位长度的椭圆柱绕自身轴线的转动惯量, 其具体数值可以根据椭圆柱的长短轴大小与密度计算得到. 根据微分方程(27)式, 可以得到转动角速度随时间的变化关系. 特别地, 当时间足够长时, 椭圆柱受到的声辐射力矩和黏滞力矩恰好相互抵消, 从而达到了稳定转动的状态, 这一稳定转动的角速度ω_zst可以通过令(27)式右边为零得到, 其计算结果为

${\omega _{z{\rm{st}}}} = \frac{{{\rm{\pi }}{b^2}{E_0}{\tau _z}}}{{f\eta }}.$

注意到(28)式是周围流体黏度的显函数, 因而在已知声辐射力矩和最终稳定转动角速度的情况下可以对流体的黏度进行反演.

本文从声波的散射理论出发, 利用部分波展开法将Gauss声束进行级数展开, 根据刚性椭圆柱和非刚性椭圆柱表面不同的边界条件计算得到声散射系数, 进而给出了以任意角度入射的Gauss声束对无限长离轴椭圆柱的声辐射力矩解析式. 在此基础上, 对水中的刚性椭圆柱和PDMS-TBE非刚性椭圆柱进行了大量的数值仿真, 得到如下结论.
1)无论是刚性还是非刚性的情况, τ_z均会出现正值或负值, 这依赖于kb大小和椭圆柱的具体位置, 因而椭圆柱的逆时针与顺时针转动均可以实现. 在倾斜入射的情况下, 即使椭圆柱位于声束中心也存在不为零的声辐射力矩.
2)对于刚性椭圆柱而言, a/b < 1和a/b > 1时的声辐射力矩的变化规律类似但方向相反. 对于PDMS-TBE椭圆柱而言, τ_z关于椭圆柱的长短轴大致呈奇对称, 并且其峰值与不同阶的共振散射模式一一对应. 在低频范围内, 刚性椭圆柱的τ_z远大于PDMS椭圆柱.
3) Gauss声束的入射角度α会对刚性椭圆柱的声辐射力矩产生显著影响, 但PDMS-TBE椭圆柱的声辐射力矩与声场频率有更密切的依赖关系, 而与α则几乎无关. 随着束腰半径kW₀的增加, 椭圆柱的散射截面随之增大, 从而增加了τ_z大小.
4)椭圆柱本身的声吸收效应会使高频范围内的声辐射力矩略微减小. 对于实际情况下非理想流体中的椭圆柱而言, 在声辐射力矩和周围流体黏滞力矩的共同作用下, 可以根据椭圆柱最终稳定转动的角速度ω_zst反演出流体的黏度η.
本研究结果预期可以为粒子的声操控特别是实现粒子的可控旋转、无容器测量技术提供理论指导, 在医学超声、材料科学等领域具有潜在的应用价值.