1.   引　言

作为一种有效且无创的非接触式操作细胞和微粒的方法, 声操控有望成为无创治疗、提高药物疗效的有效工具. 在过去的十几年间, 声操控在医学领域中得到广泛应用^[1,2]. 除了用作医生病理分析诊断的辅助工具外^[3], 许多研究人员发现了声操控技术在医学治疗的广阔前景^[4]. Wu等^[5]通过实验对比了有无超声条件下药物的输送率, 结果显示利用超声确实能提高将药物输送到指定位置的传输效率. Wang等^[6]基于声流诱导碰撞建立了一种细胞裂解装置, 在20 s工作时间内实现了95%的裂解率. Mishra等^[7]通过实验证明辐射力也会使红细胞变形且与光镊产生的变形程度相当. Silva等^[8]使用多种波长的超声波驻波装置将囊泡变形, 结合形变图像与薄壳弹性理论获得囊泡的杨氏模量等力学性质. Zhang等^[9]通过研发一种装置将平面声波转换成涡流束来旋转虾卵, 并提出根据细胞旋转行为确定细胞参数的可能. 细胞的声操控有着广阔的应用前景, 是未来值得探索的领域.

细胞的声操控是实现细胞筛选、分离的关键, 声操控应用的基础是声辐射力(ARF)的控制. 悬浮在声场中的粒子受到由声波散射而产生的时间平均力称为声辐射力^[10]. 目前已经建立了许多模型用于计算不同声场中的声辐射力^[11-14]. King^[15]首先提出了刚性球面上的声辐射力公式. Rajabi等^[16]研究了可平移振荡的刚性振动球上出现的负向声辐射力, 为实现平面波操控粒子平移运动提供了理论支持. 考虑到细胞是带有膜结构的, 结合薄壳弹性理论来建立模型, 而关于在无界空间中作用在弹性球壳的声辐射力已有很多研究^[17-19]. Hasegawa等^[20]从理论上研究了浸没在流体中的球壳和柱壳上所产生的声辐射力, 指出随壳的厚度增大球壳所受声辐射力幅值减小. Junger^[21]建立了平面波声场中薄弹性壳的声散射理论, 发现外部流体介质的阻尼效应限制球壳的共振响应. 随后, Mitri^[22]推导了平面驻波场中的弹性球壳的声辐射力函数表达式, 指出弹性球壳的散射特性会因外部液体介质改变而发生变化. 随着对细胞操控和靶向治疗研究的深入, 陆续有研究者提出采用两层和三层球模型来模拟细胞和药粒, 以探究该情况下细胞和药粒的声辐射力. Wang等^[23]研究了高斯驻波场中多层球面上的声辐射力理论, 求解了双层球以及三层球模型的声辐射力表达式. Wang等^[24]提出用三层球模型来模拟细胞结构, 随着内芯半径或外壳厚度的增大模型所受声辐射力明显增大.

在实际中细胞核在细胞中的位置会发生改变, 因此需要考虑细胞核偏心对细胞所受声辐射力的影响. Thompson^[25] 首先提出了偏心球模型并计算了偏心球源的声辐射力, 随后Roumeliotis等^[26]发现由于引入偏心内球而引起的硬壁球形腔声共振频率偏移. Hasheminejad和Azarpeyvand^[27]对内置的偏心辐射源散射场进行研究, 指出密封外壳的存在导致模态阻力和惯性振幅的明显增强. Zang等^[28]建立理论模型探究了位于零阶贝塞尔声束中双层偏心液滴中偏心球的位置对负向声辐射力产生的影响, 指出当偏心球位于液体球心右侧时会产生负向声辐射力. 通过对不同类型的细胞进行操控和分选以获得某些组织的病变信息是医学诊断的重要手段, 而现有的理论模型缺少对细胞核位置及大小变化的准确描述, 这就使得细胞核在发生异常的增大或偏移时细胞周围声散射变化被忽略. 为此, 需要发展一种更接近实际细胞结构的模型来估算细胞所受到的声辐射力, 从而提高诊断的可信度. 所以本文在Zang等^[28]的模型的基础上, 发展了一个内含偏心液滴的弹性充液球壳的三层球模型来模拟药粒和细胞, 并重点讨论了偏心液滴的不同参量对充液球壳所受的声辐射力的影响. 本工作扩展了声辐射力的理论, 有助于不同类型细胞的操控和分选, 病变组织的区分检测及靶向治疗的发展.

2.   理论模型

本文的研究对象是内部含有偏心液滴的球形密闭液体腔, 液体腔被厚度为h的弹性薄球壳包裹构成的复合粒子, 如图1所示, 偏心液滴和弹性球壳半径分别为$a$和$b$, 球壳外部为无界理想液体. 为方便分析, 建立两套球坐标系分析球形密闭液体腔在平面波声场中的声散射和声辐射力. 以球形偏心液滴中心为坐标原点建立球坐标系O₁, 坐标为$\left( {{r_1}, {\theta _1}} \right)$. 设球形液滴中心与液体腔中心间的偏心距离d, 以液体腔中心为原点建立球坐标系O₂, 坐标为$\left( {{r_2}, {\theta _2}} \right)$. 根据液体腔的几何形貌, 将液体腔分为3个区域, 偏心液滴内部为区域Ⅰ, 液滴外的液体球腔为区域Ⅱ, 薄球壳外部为区域Ⅲ, 其中液滴、腔内液体、外壳和外部流体密度及声速分别为(ρ₁, c₁), (ρ₂, c₂), (ρ_s, c_s), (ρ₃, c₃).

假定流体介质是可压缩无黏性的理想流体. 利用球坐标系下的分波级数法, 为简化分析, 忽略时谐因子${{\text{e}}^{ - {\text{i}}\omega t}}$, 则薄球壳外部入射平面波及其散射声场可表示为

$${p_{{\text{int}}}} = {p_0}{\text{ }}\sum\limits_{n{\text{ = 0}}}^\infty {{\gamma _n}{{\text{j}}_n}({k_3}{r_2}){P_n}(\cos {\theta _2})} ,$$

(1)

$${p_{{\text{sca}}}} = {p_0}{\text{ }}\sum\limits_{n{\text{ = 0}}}^\infty {{\gamma _n}{C_n}{\text{h}}_n^{(1)}({k_3}{r_2}){P_n}(\cos {\theta _2})} ,$$

(2)

其中γ_n = (2n+1)iⁿ, ${{\text{j}}_n}\left( \cdot \right)$表示n阶球Bessel函数, ${\text{h}}_n^{(1)}\left( \cdot \right)$表示n阶第一类球Hankel函数, P_n表示勒让德多项式, C_n为声散射系数. 因此, 弹性球壳外部的总声压为

$$\begin{split} {p_3} =\;& {p_0}{\text{ }}\sum\limits_{n{\text{ = 0}}}^\infty {\gamma _n}\Big[ {{\text{j}}_n}({k_3}{r_2}) + {C_n}{\text{h}}_n^{(1)}({k_3}{r_2}) \Big] \\& \times {P_n}(\cos {\theta _2}) , \end{split}$$

(3)

式中k_i = ω/c_i (i = 1, 2, 3)为声波在对应区域介质中的波数. 弹性球壳内偏心液滴外部液体中的声压可表示为

$$\begin{split} {p_2} =\;& {p_0}{\text{ }}\sum\limits_{n{\text{ = 0}}}^\infty {\gamma _n}\Big[ {E_n}{{\text{j}}_n}({k_2}{r_2}){P_n}(\cos {\theta _2}) \\& +{B_n}{\text{h}}_n^{(1)}({k_2}{r_1}){P_n}(\cos {\theta _1}) \Big] , \end{split}$$

(4)

同理可得偏心液滴内部声压为

$${p_1} = {p_0}{\text{ }}\sum\limits_{n{\text{ = 0}}}^\infty {{\gamma _n}{D_n}{{\text{j}}_n}({k_1}{r_1}){P_n}(\cos {\theta _1})} .$$

(5)

球坐标系下弹性壳的位移用球壳的中间层的偏移量表示, 故基于球壳的广义径向和切向位移分量的傅里叶级数式可表示为

$$W(\omega ,{\theta _2}) = \sum\limits_{n{\text{ = 0}}}^\infty {{W_n}{P_n}(\cos {\theta _2})} ,$$

(6)

$$V(\omega ,{\theta _2}) = \sum\limits_{n{\text{ = 0}}}^\infty {{V_n}P_n^1(\cos {\theta _2})} ,$$

(7)

式中$P_n^1\left( x \right) = {\left( {1 - {x^2}} \right)^{1/2}}{\text{d}}{P_n}\left( x \right)/{\text{d}}x$, W_n和V_n分别是球壳中间层的径向和切向位移分布, 在压力场中满足的方程为^[27]

$$\begin{split}& \left[{\varOmega }^{2}-(1+{\beta }^{2})(\upsilon +{\lambda }_{n}-1)\right]{V}_{n}\\&-\left[{\beta }^{2}(\upsilon +{\lambda }_{n}-1)+(1+\upsilon )\right]{W}_{n}=0\text{, }~~~~ \end{split}$$

(8)

$$\begin{split} & -{\lambda }_{n}\left[{\beta }^{2}(\nu +{\lambda }_{n}-1)+(1+\nu )\right]{V}_{n}\\& +\left[{\varOmega }^{2}-2(1+\nu )-{\beta }^{2}{\lambda }_{n}(\nu +{\lambda }_{n}-1)\right]{W}_{n}\\&= -\frac{{a}^{2}{(1-\nu )}^{2}}{Eh}\Delta {P}_{n}(\omega )\text{, } \end{split}$$

(9)

式中, $\beta^2=h^2/(12b^2)$, ${\lambda _n} = n(n + 1)$, $\varOmega = \omega b/{c_{\text{p}}}$为无量纲参数, $c_{\text{p}}^2 = E/(1 - {\upsilon ^2}){\rho _{\text{s}}}$是弹性外壳的纵波相速度, ${\rho _{\text{s}}}$为弹性壳密度, h为弹性壳厚度, E是弹性壳杨氏模量, $\nu$为泊松比, ΔP_n(ω)是球壳内外声压差. 联立(8)式和(9)式, 可得:

$${W_n}{\text{ = }} - \frac{{{a^2}{{(1 - \nu )}^2}}}{{Eh{Z_n}}}\Delta {P_n}(\omega ),~~~~~~~~~~~~~~~~~~~$$

(10)

$${V_n}{\text{ = }}\frac{{{\beta ^2}(\nu + {\lambda _n} - 1) + (1 + \nu )}}{{{\varOmega ^2} - (1 + {\beta ^2})(\nu + {\lambda _n} - 1)}}{W_n},~~~~~~~~$$

(11)

$$\begin{split} {Z_n} =\;& \left[ {{\varOmega ^2} - 2(1 + \nu ) - {\beta ^2}{\lambda _n}(\nu + {\lambda _n} - 1)} \right] \\& -{\lambda _n}[ {\beta ^2}(\nu + {\lambda _n} - 1) \\&+(1 + \nu ) ]\frac{{{\beta ^2}(\nu + {\lambda _n} - 1) + (1 + \nu )}}{{{\varOmega ^2} - (1 + {\beta ^2})(\nu + {\lambda _n} - 1)}}, \end{split}$$

(12)

$$\begin{split} \Delta {P_n}(\omega ) =\;& {\gamma _n}\bigg[ {{\text{j}}_n}({k_3}{a_ + }) + {C_n}{\text{h}}_n^{(1)}({k_3}{a_ + }) - {E_n}{{\text{j}}_n}({k_2}{a_ - }) \\&- \sum\limits_{m = 0}^\infty {{B_m}} {\text{h}}_n^{(1)}({k_2}{a_ - }){Q_{mn}}({k_2}d) \bigg], \\[-18pt] \end{split}$$

(13)

式中, ${a_ - } = a - h/2$和${a_ + } = a + h/2$分别为表示弹性球壳内壁面和外壁面径向坐标. 在应用边界条件过程中需要将各声参量利用平移加法定理^[29]在球坐标O₁和O₂中进行转换, 因此在(13)式中引入了球坐标转换系数Q_mn, 其具体表达式见附录A.

待定散射系数C_n, B_n, E_n, D_n可由Ⅰ, Ⅱ流体界面和Ⅱ, Ⅲ弹性球壳界面的边界条件确定. 在液滴表面处(r₁=b)声压和速度连续, 有

$$\begin{split} {D_n}{{\text{j}}_n}({k_1}b) =\;& \sum\limits_{m = 0}^\infty {{E_m}} {{\text{j}}_n}({k_2}b){Q_{mn}}({k_2}d) \\& +{B_n}{\text{h}}_n^{(1)}({k_2}b), \end{split}$$

(14)

$$\begin{split} \frac{{{\rho _2}{c_2}}}{{{\rho _1}{c_1}}}{D_n}{{\text{j}}_n'}({k_1}b) =\;& \sum\limits_{m = 0}^\infty {{E_m}} {{\text{j}}_n'}({k_2}b){Q_{mn}}({k_2}d) \\&+ {B_n}{\text{h}}_{n}^{(1)'}({k_2}b), \end{split}$$

(15)

联立 (14)式和(15)式, 系数B_n可表示为

$${B_n} = {M_n}\sum\limits_{m = 0}^\infty {{E_m}} {Q_{mn}}({k_2}d),$$

(16)

即${B_n}$可用${E_m}(m = 0, 1, 2, \cdots )$的求和表示, 且有

$${M_n} = \frac{{{{\text{j}}_n}({k_1}b){{{\text{j}}}_n'}({k_2}b) - \dfrac{{{\rho _2}{c_2}}}{{{\rho _1}{c_1}}}{{{\text{j}}}_n'}({k_1}b){{\text{j}}_n}({k_2}b)}}{{\dfrac{{{\rho _2}{c_2}}}{{{\rho _1}{c_1}}}{\text{h}}_n^{(1)}({k_2}b){{{\text{j}}}_n'}({k_1}b) - {\text{h}}_{n}^{(1)'} ({k_2}b){{\text{j}}_n}({k_1}b)}}.$$

(17)

在r₂ = a处, 由球壳内外表面的径向速度连续性, 得

$$\begin{split} & - {\text{i}}\omega {W_n}\left( {\omega ,\theta {}_2} \right) \\ =\;& \frac{1}{{{\text{i}}\omega {\rho _3}}}{\left. {\frac{{\partial {p_3}}}{{\partial {r_2}}}} \right|_{{r_2} = {a_ + }}} = \frac{1}{{{\text{i}}\omega {\rho _2}}}{\left. {\frac{{\partial {p_2}}}{{\partial {r_2}}}} \right|_{{r_2} = {a_ - }}}, \end{split}$$

(18)

将 (3)式和(4)式代入(18)式, 有

$$\omega {W_n} - \frac{{{\gamma _n}}}{{{\rho _3}{c_3}}}\left[ {{{{\text{j}}}_n'}({k_3}{a_ + }) + {C_n}{\text{h}}{{_n^{\left( 1 \right)'}} }\left( {{k_3}{a_ + }} \right)} \right] = 0,$$

(19)

$$\begin{split} & \omega {W_n} - \frac{{{\gamma _n}}}{{{\rho _2}{c_2}}}\Bigg[ {E_n}{{{\text{j}}}_n'}({k_2}{a_ - }) \\&+ \sum\limits_{m = 0}^\infty {B{}_m} {\text{h}}{{_n^{\left( 1 \right)'}} }\left( {{k_2}{a_ - }} \right)Q{}_{mn}\left( {{k_2}d} \right) \Bigg] = 0. \end{split}$$

(20)

将(16)式代入(19)式和(20)式, 得

$${S_n} + {R_n}{C_n} - {g_n}{\text{h}}_n^{\left( 1 \right)}\left( {{k_2}{a_ - }} \right)\sum\limits_{m = 0}^\infty {{M_m}\left( {\sum\limits_{j = 0}^\infty {\frac{{{\eta _j} + {\varsigma _j}{C_j}}}{{{d_j}}}{Q_{jm}}\left( {{k_2}d} \right)} } \right){Q_{mn}}\left( {{k_2}d} \right)} = 0,$$

(21)

即为C_n (n = 0, 1, 2, ···)的求解式, 方程数量与n的取值相关, 且有:

$${S_n} = \frac{{{\eta _n}}}{{{d_n}}} + \frac{{{{{\text{j}}}_n'}({k_3}{a_ + })}}{{{\rho _3}{c_3}}} + {g_n}{{\text{j}}_n}({k_3}{a_ + }),~~~~{R_n} = \frac{{{\varsigma _n}}}{{{d_n}}} + \frac{{{\text{h}}{{_n^{\left( 1 \right)'}} }\left( {{k_3}{a_ + }} \right)}}{{{\rho _3}{c_3}}} + {g_n}{\text{h}}_n^{\left( 1 \right)}({k_3}{a_ + }),$$

$$\begin{split} {\eta _n} = \frac{1}{{{\rho _2}{c_2}}}{{\text{j}}_n'} ({k_3}{a_ + }){\text{h}}_n^{(1)'} \left( {{k_2}{a_ -}} \right) + \frac{{{\rho _3}{c_3}}}{{{\rho _2}{c_2}}}{g_n}{{\text{j}}_n}({k_3}{a_ + }){\text{h}}_{n}^{\left( 1 \right)'} \left( {{k_2}{a_ - }} \right) - {g_n}{{\text{j}}_n'}({k_3}{a_ + }){\text{h}}_n^{\left( 1 \right)}\left( {{k_2}{a_ - }} \right), \end{split}$$

$$\begin{split} {\varsigma _n} =\;& \frac{1}{{{\rho _2}{c_2}}}{\text{h}}_{n}^{\left( 1 \right)'} \left( {{k_3}{a_ + }} \right){\text{h}}{_n^{\left( 1 \right)'}}\left( {{k_2}{a_ - }} \right) + \frac{{{\rho _3}{c_3}}}{{{\rho _2}{c_2}}}{g_n}{\text{h}}_n^{\left( 1 \right)}\left( {{k_3}{a_ + }} \right){\text{h}}_{n}^{\left( 1 \right)'} \left( {{k_2}{a_ - }} \right) - {g_n}{\text{h}}_{n}^{\left( 1 \right)'} \left( {{k_3}{a_ + }} \right){\text{h}}_n^{\left( 1 \right)}\left( {{k_2}{a_ - }} \right), \end{split}$$

$$\begin{split} {d_n} =\;& \frac{{{\rho _3}{c_3}}}{{{\rho _2}{c_2}}}{g_n}\left[ {{\text{j}}_n}({k_2}{a_ - }){\text{h}}{{_n^{\left( 1 \right)'}} }\left( {{k_2}{a_ - }} \right) - {{{\text{j}}}_n'}({k_2}{a_ - }){\text{h}}_n^{\left( 1 \right)}\left( {{k_2}{a_ - }} \right) \right], ~~~~ {g_n} = \frac{{\omega {a^2}}}{{c_{\text{p}}^2h{\rho _{\text{s}}}{Z_n}}}. \end{split}$$

平面波声场中复合粒子受到的声辐射力取决于粒子的散射声场. 基于方程组(21)可得到弹性球壳的n阶声散射系数C_n, 若n趋于无穷, 则为一无限维方程组. 在求解过程中, 通常取一截断n值后通过求解线性方程组的方法确定系数C_n.

3.   声辐射力

在无界理想液体中, 平面声波场中球状粒子所受声辐射力常用附加压力的时间平均值的球面积分表示, 与声波在介质中形成的非线性扰动相关, 即弹性球壳上的声辐射力函数为

$$\left\langle {\boldsymbol{F}} \right\rangle = - \iint_{S(t)} {\left\langle {{p_3} - {p_0}} \right\rangle }{{\boldsymbol{e}}_r}{\text{d}}S,$$

(22)

忽略三阶小量后用速度势表示声压为

$${p_3} - {p_0} = {\rho _3}\frac{{\partial \psi }}{{\partial t}} - \frac{1}{2}{\rho _3}{\nu ^2} + \frac{1}{2}\frac{{{\rho _3}}}{{{c_3}}}{\left( {\frac{{\partial \psi }}{{\partial t}}} \right)^2},$$

(23)

基于声场中声压扰动与速度势以及速度场分布相关的二阶场效应, 将(23)式代入(22)式力函数为^[30]

$$\begin{split} \left\langle {\boldsymbol{F}} \right\rangle =\;& - \left\langle {\iint_{{S_0}} {{\rho _3}\left( {{\upsilon _r}{{\boldsymbol{e}}_r} + {\upsilon _\theta }{{\boldsymbol{e}}_{\boldsymbol{\theta }}}} \right){\upsilon _r}{\text{d}}S}} \right\rangle \\&+ \left\langle {\iint_{{S_0}} {\frac{{{\rho _3}}}{2}{{\left| \upsilon \right|}^2}{{\boldsymbol{e}}_r}{\text{d}}S}} \right\rangle \\&- \left\langle {\iint_{{S_0}} {\frac{{{\rho _3}}}{{2{c_3}^2}}{{\left( {\frac{{\partial \psi }}{{\partial t}}} \right)}^2}{{\boldsymbol{e}}_r}{\text{d}}S}} \right\rangle , \end{split}$$

(24)

其中, e_r, e_θ分别是单位矢量的径向分量和切向分量, S₀为散射体边界表面. 在边界(${r_2} = {a_ + }$)处粒子速度的径向分量和切向分量分别为

$${\upsilon _r} = - \frac{{\partial \psi }}{{\partial r}},$$

(25)

$${\upsilon _\theta } = - \frac{1}{r}\frac{{\partial \psi }}{{\partial \theta }},$$

(26)

式中$\psi$是充液球壳外部液体环境中一阶散射场速度势$\varphi$的实部, 即

$$\psi = {{\rm{Re}}} \left[ \varphi \right] = \left| {{\varphi _0}} \right|\sum\limits_{n = 0}^\infty {\left( {2n + 1} \right){R_n}{P_n}\left( {\cos {\theta _2}} \right)} ,$$

(27)

式中${R_n} = {{\rm{Re}}} \left[ {{{\left( { - {\text{i}}} \right)}^n}\left( {{U_n} + {\text{i}}{V_n}} \right){{\text{e}}^{{\text{i}}\omega t}}} \right]$, 已知速度势与声压之间有转换关系$\varphi = \dfrac{{{p_3}}}{{{\text{i}}\omega {\rho _3}}}$, ${\varphi _0} = \dfrac{{{p_0}}}{{{\text{i}}\omega {\rho _3}}}$, 在进行无量纲化的处理过程中${\varphi _0}$会消掉, 故无论是用声压还是速度势表示散射场, 声辐射力Y_p的值都相等. 在散射声场中有:

$${U_n} = \left( {1 + {\alpha _n}} \right){{\text{j}}_n}\left( x \right) - {\beta _n}{{\text{y}}_n}\left( x \right),$$

(28)

$${V_n} = {\beta _n}{{\text{j}}_n}\left( x \right) + {\alpha _n}{{\text{y}}_n}\left( x \right),$$

(29)

式中${\alpha _n} = {{\rm{Re}}} \left[ {{C_n}} \right],\; {\beta _n} = {{\rm{Im}}} \left[ {{C_n}} \right],\; x = {k_3}{r_2}$$,\; {{\text{y}}_n}\left( \cdot \right)$表示n阶第二类球Bessel函数.

因此, 在平面波传播方向(x方向)上充液球壳所受的声辐射力分量可表示为

$$\left\langle {{F_x}} \right\rangle = \left\langle {{F_r}} \right\rangle + \left\langle {{F_\theta }} \right\rangle + \left\langle {{F_{r\theta }}} \right\rangle + \left\langle {{F_t}} \right\rangle ,$$

(30)

其中:

$$\begin{split} & \left\langle {{F_r}} \right\rangle = - 2{\text{π}}\rho {x^2}{\varphi _0}^2\sum\limits_{n = 0}^\infty \left( {n + 1} \right)\left( {{V}_n'}{{U}_{n + 1}'} - {{U}_n'}{{V}_{n + 1}'} \right) , \\& \left\langle {{F_\theta }} \right\rangle = 2{\text{π}}\rho {\varphi _0}^2\sum\limits_{n = 0}^\infty n\left( {n + 1} \right)\left( {n + 2} \right)\left( {V_n}{U_{n + 1}} - {U_n}{V_{n + 1}} \right) , \\& \left\langle {{F_{r\theta }}} \right\rangle = 2{\text{π}}\rho x{\varphi _0}^2\sum\limits_{n = 0}^\infty \big[ n\left( {n + 1} \right)\left( {V_n}{{U'}_{n + 1}} - {U_n}{{V'}_{n + 1}} \right) \\ & \qquad \qquad - \left( {n + 1} \right)\left( {n + 2} \right)\left( {{V'}_n}{U_{n + 1}} - {{U'}_n}{V_{n + 1}} \right) \big] , \\& \left\langle {{F_t}} \right\rangle = - 2{\text{π}}\rho {x^2}{\varphi _0}^2\sum\limits_{n = 0}^\infty \left( {n + 1} \right)\left( {V_n}{U_{n + 1}} - {U_n}{V_{n + 1}} \right). \end{split}$$

(31)

引入无量纲声辐射力函数分析粒子在声场中受到的力效应, 即

$${Y}_{\text{p}}= \dfrac{\langle {F}_{x}\rangle}{{S}_{\text{c}}\overline{E} },$$

(32)

式中$\overline E = \dfrac{1}{2}{\rho _3}k_3^2{\left| {{\varphi _0}} \right|^2}$为入射平面波的平均能量密度, ${S_{\text{c}}} = {\text{π}}r_2^2$是球壳的横截面积. 将(30)式代入(32)式, 无量纲声辐射力函数Y_p可表示为

$${Y}_{\text{p}}=\frac{4}{{x}^{2}}{\displaystyle \sum _{n=0}^{\infty }\left({A}_{n}+{B}_{n}+{C}_{n}+{D}_{n}\right)}\text{, }$$

(33)

其中:

$$\begin{split}& {A_n} = \left( {n + 1} \right)\left( {{{V}_n'}{{U}_{n + 1}'} - {{U}_n'}{{V}_{n + 1}'}} \right){x^2}, \\& {B_n} = - n\left( {n + 1} \right)\left( {n + 2} \right)\left( {{V_n}{U_{n + 1}} - {U_n}{V_{n + 1}}} \right), \\& {C_n} = \left[ n\left( {n + 1} \right)\left( {{U_n}{{V}_{n + 1}'} - {V_n}{{U}_{n + 1}'}} \right)\right. \\&\qquad \left.- \left( {n + 1} \right)\left( {n + 2} \right)\left( {{{U}_n'}{V_{n + 1}} - {{V}_n'}{U_{n + 1}}} \right) \right]x, \\& {D_n} = \left( {n + 1} \right)\left( {{V_n}{U_{n + 1}} - {U_n}{V_{n + 1}}} \right){x^2}. \end{split}$$

将(28)式和(29)式代入(33)式即可计算复合粒子的声辐射力函数. 粒子受到的声辐射力与粒子属性密切相关, 对于充液弹性球壳内含液性微粒, 其物质分布不均匀, 在声波作用下的响应同时取决于球壳介质、壳内液体以及液滴声学特性, 同时还与结构特征有关. 为分析声波对非均匀球形液体粒子的操控影响, Zang等 ^[28]研究构建了贝塞尔声场中偏心粒子的辐射力模型, 分析了小球的偏心对声辐射力的影响, 发现偏心小球位于液体球心右侧时在高频范围内粒子会受到负向声辐射力, 而偏心小球在液体球心左侧时粒子会受到正向的声辐射力. 偏心球粒子的材料种类改变时, 粒子所受声辐射力有巨大差异.

4.   数值分析

为拓展薄壳包裹的非均匀液性粒子在声操控应用中的基础理论, 分析粒子结构以及物质特性对平面波场中复合粒子所受声辐射力, 本文选取材料为钢的弹性薄壳包裹的复合液滴以及基于细胞结构模型化后的复合微粒结构开展数值分析, 相关参数见附录B.

4.1   充液不锈钢弹性薄壳复合粒子所受声辐射力

为讨论薄壳复合粒子在平面波声场中受到的声辐射力, 基于(33)式分析了液体特性的影响. 薄球壳材料采用不锈钢, 外部流体为水, 区域Ⅱ内充满汞, 区域Ⅰ为偏心放置的水滴(d/a = 0.3), 水滴与球壳半径比为b/a = 0.3. 图2给出了壳层厚度分别为h/a = 0.01, 0.05和0.10三种情形下的声辐射力随无量纲粒子半径k₃a的变化曲线. 为表述方便, 后文中将k₃a简写为ka. 通过比较可以看出, 在 0.7 < ka < 1.6范围内, 3种壳层厚度的粒子对应的辐射力函数曲线均出现明显共振峰, 但峰值分布特征不同. 当ka = 0.7时, 壳层厚度0.01的声辐射力函数曲线出现第一个峰值, ka = 1.1时出现第二个峰值; 壳层厚度0.05的声辐射力函数曲线在ka = 1.3时达到峰值; 壳层厚度0.1时声辐射力函数曲线在ka = 1.3附近存在负向声辐射力. 随着壳层厚度的增大, 辐射力函数曲线峰值逐渐变小. 峰值位置不同反映了声场中复合粒子声响应特性的变化, 由于考虑了粒子薄壳层的振动响应的影响, 随着弹性球壳厚度增大, 壳层的径向振动响应减弱, 对入射声波能量的共振散射能力减弱, 进而导致辐射力共振峰峰值下降. 比较基于声波扰动下薄球壳应力应变响应理论相关的复合粒子所受辐射力曲线发现, 壳层厚度变化将影响粒子辐射力的共振响应, 主要表现为壳层越薄, 粒子在低频区的共振峰越多且越显著, 这可能与复合结构粒子各个组成部分的物质结构特性相关. 球壳厚度增加导致与径向振动相关的形变减小, 可在一定程度上抑制特定频率范围内的背向散射, 进而导致在球壳厚度为h/a = 0.1时出现的负向辐射力分布. 随着无量纲参量ka的增大, 薄壳厚度变化对辐射力的影响在一定的频率范围内呈现出减弱的现象, 如在2.4 < ka < 3.8范围内, 3种壳层厚度的粒子所受声辐射力差异较小, 这可能与该频段内频率增大而波长减小, 背向散射主要取决于粒子尺度大小有关. 随着ka的继续增大, 薄壳厚度影响又逐渐增大, 这可能与声波引起的结构响应的耦合效应有关, 即复合结构的变化会改变粒子的一阶散射场分布以及各级共振散射特性, 进而导致辐射力函数曲线共振峰位置分布不同. 事实上, 随着壳层厚度增大, 壳层声响应模式也不同, 如弹性厚壳内需考虑横波影响, 因此, 相比较而言, 薄壳振动理论更适用于h/a较小的情形, 因此, 在后面的分析中, 选择球壳厚度为h/a = 0.01.

复合粒子所处的液体环境不同, 辐射力响应曲线特征不同, 分析了外部流体介质Ⅲ分别为甘油和汞时(其他参数与图2一致)的声辐射力函数曲线, 如图3所示. 对比发现, 辐射力函数曲线均呈现出不同程度的峰值响应, 且主峰值的位置随着液体不同而变化, 甘油中主峰值的位置在水银的右侧, 且甘油中能够得到较大的峰值. 随着ka的增大, 辐射力函数起伏变化, 说明辐射力对粒子尺寸和频率敏感. 从辐射力函数数值看, 同等条件下介质Ⅲ为甘油时声辐射力峰值远大于介质Ⅲ为汞, 特别是在ka较小的范围内, 且外部介质为水银时在ka = 1.1附近区域内出现负向辐射力. 因此, 改变外部流体介质阻抗时, 会改变球壳的共振模式, 且周围介质阻抗与壳内介质阻抗差异越大, 球壳的散射能力越强, 高阻抗梯度的介质环境可以产生较大的辐射力. 在内部介质为水银, 外部介质为甘油时, 内部与外部介质阻抗梯度极大, 内部可视为硬边界, 这意味着有更多的入射波能量被散射. 从共振峰的位置看, 主要取决于入射声波频率、粒子内外介质特性等, 因此, 在图3(a)出现明显的辐射力函数的主共振、二分频共振和三分频共振现象, 主共振峰出现在ka = 5.3附近, 二分频共振共振峰在ka = 3.1附近, 三分频共振共振峰在ka = 1.6附近, 从数值看, 峰值位置会发生一定程度的漂移, 这可能与辐射力函数的非线性特性相关. 两种液体环境内看辐射力函数曲线在低频区均表现出比较显著的起伏特征, 因此, 从复合粒子的筛选, 可以考虑利用小于主共振频率的低频声波实现.

4.2   理想流体中类细胞结构的声辐射力

因细胞内部有细胞质和细胞核等, 可视作一充液壳形结构^[8], 当利用声波对细胞进行声操控时, 需考虑细胞结构对操控的潜在影响. 基于细胞特征模型化后的复合粒子对应介质设置为: 用聚糖壳模拟细胞膜, 区域Ⅲ介质为水, 区域Ⅱ, Ⅰ介质分别为细胞质和细胞核. 为探究壳内微粒偏心距对细胞所受声辐射力的影响, 将细胞核与细胞半径比b/a固定为0.6, 比较不同偏心距下声辐射力曲线, 如图4所示, 其中壳内微粒偏心距d/a = 0, 0.05, 0.10和0.20. 对比发现, 当壳内微粒处于偏心状态时可能影响辐射力函数在ka < 1范围内的响应状态, 即出现多个峰值位置, 即低频区也可能激励出较强的辐射力. 在ka > 2的范围内观察到辐射力函数较为明显的起伏变化, 腹点位置的分布特征表现为随着偏心距的增大, 在2 < ka < 3范围内的腹点位置有细微的左移趋势, 而在5 < ka < 6之间的腹点位置有较为明显右移趋势, 且偏心程度越大, 辐射力函数值越大, 意味着壳内微粒的偏心效应对辐射力幅值以及分布状态均有一定影响. 细胞核偏心程度的提高改变了细胞内质量分布, 细胞质心位置向相应位置的细胞膜偏移, 声波引起的细胞膜振荡的不稳定性受到细胞核的影响也增强, 同时, 细胞结构的变化导致共振频率的变化, 进而影响细胞在不同频率范围内的声散射特征, 声辐射力响应也随之变化. 因此, 在利用超声技术辅助病情诊断时, 可以根据细胞散射特性随细胞核的位置发生变化这一特征, 借助声操控分离出正常细胞和病变细胞从而达到诊断目的^[31].

组织发生病变时细胞往往出现异常的增殖, 为了适应细胞增殖速度, 癌变细胞的细胞核会出现异常的增大^[32]. 为了探究这一现象, 了解细胞核大小对细胞周围散射场的影响, 在保持细胞核位于细胞中心且其他参数不变的情形下分析了细胞核与细胞半径比的影响, 分别设置半径比b/a为0.5, 0.6和0.7时的辐射力函数曲线如图5所示. 对比发现, 随着细胞核半径的增大, 辐射力曲线出现明显差异, 主要表现为当b/a = 0.7时在ka < 3范围内的局域共振峰值点增多的现象, 表明对确定尺寸的壳形细胞结构而言, 内部因异常增生形成的核结构将影响其在低频区受到的辐射力的共振响应状态, 在局域峰值点附近的声波频率的驱动下可能更有利于细胞的筛选. 在ka>3的范围内, 细胞受到的声辐射力强度随ka变化起伏变化且腹点位置有左移的趋势.

为进一步探究细胞核偏心与否和细胞内核尺寸的影响, 对比分析偏心距分别为d/a = 0和d/a = 0.1时细胞结构模型的声辐射力曲线, 如图6所示. 声辐射力函数曲线在ka < 3的范围内出现多个共振峰, 表明有局域共振响应发生. 当壳层结构内核向入射声场一侧偏移时, 偏移量将影响辐射力的起伏变化特征, 即可能出现局域略微增强的情形. 在ka < 1的范围内出现了小范围的负辐射力函数值分布区, 说明随着结构的变化在低频区的响应有利于不同结构特征的复合壳层粒子的分离. 当细胞核偏离细胞中心时, 声辐射力变大, 相比较而言, 细胞核偏心距对细胞的声辐射力幅度变化影响更加显著. 细胞核同时发生偏离和半径变化时二者对细胞所受声辐射力的影响可形成叠加效应.

为了比较细胞质和细胞核对细胞所受到声辐射力强度的调控作用, 基于细胞核和细胞质特性阻抗展开分析. 在分析过程中将特性阻抗进行了无量纲化, 具体的细胞核阻抗Z_n1和细胞质阻抗Z_n2的取值参看附录B表B1. 图7给出的数值分析结果表明, 由于细胞周围的散射场同细胞核、细胞质的特性阻抗密切相关, 因此, 辐射力函数受细胞核以及细胞质阻抗变换的影响显著. 图7(a)中给出细胞核相对特性阻抗分别为0.8, 0.9, 1, 1.1和1.2 时的辐射力函数随ka变化曲线, 相比较而言, 其特性阻抗的变化主要影响不同ka值时辐射力的大小, 对变化趋势影响不明显. 随着细胞核阻抗的增大, 在ka = 5附近的起伏幅度逐步增大, 且腹点位置有右移的趋势. 图7(b)给出了细胞质相对特性阻抗分别为0.8, 0.9, 1.0, 1.1和1.2 时的辐射力函数曲线, 对比发现, 当细胞质特性阻抗比小于1.0时, 在1.1 < ka < 2.6范围内, 细胞受到负向的声辐射力, 阻抗比为0.8时起伏变化更加显著, 说明特性阻抗比的微小变化可对声辐射力形成影响. 对比细胞质阻抗比大于1.0的3种情况发现, 除了没有负向辐射力出现外, 三者的起伏变化特征也存在明显差异, 因此, 辐射力函数随着ka的变化对细胞质特性阻抗非常敏感. 相比较而言, 细胞质的特性阻抗对声辐射力强度影响更为显著. 当细胞质发生变化时, 壳内介质与外部介质阻抗梯度变大, 散射增强, 声辐射力也增强. Jo和Guldiken ^[33]在实验中就曾发现密度较大的细胞受力较大, 伴随着较大的横向位移从而运动到通道两侧, 密度较小的细胞横向位移可以忽略, 从而实现细胞的分选. 细胞质在细胞总体积中占比更大, 当其密度增大时, 细胞特性阻抗改变更大. 由此产生更强的声辐射力, 可以利用这一特性实现对异质细胞的筛选.

5.   结　论

本文从有核细胞的声操控出发, 提出用内置偏心液滴的充液弹性球壳模拟有核细胞, 以探究充液球壳各部分参数对其所受声辐射力的影响. 基于薄球壳理论推导了内置偏心液滴的充液弹性球壳在平面波入射声场中的声辐射力公式, 讨论了弹性球壳的厚度、内置液滴偏心距、半径大小以及外部液体环境对充液球壳所受声辐射力的影响. 弹性球壳的厚度越大, 充液球壳所受声辐射力越小. 随内置液滴偏心距的增大, 充液球壳所受声辐射力增大. 偏心液滴半径大小变化时, 充液球壳所受声辐射力曲线在ka < 3范围内共振峰增多, 在ka > 3范围内极值点有左移的趋势. 以上结果表明, 内置液滴位置和大小是影响充液球壳所受声辐射力的重要因素, 即可以根据充液球壳散射特性对腔内液滴的大小和位置进行判断. 液体腔内介质及偏心液滴的特性阻抗的改变会影响声辐射力强度, 二者对比发现液体腔内介质的特性阻抗对于充液球壳所受声辐射力的变化起着主要调配作用, 其特性阻抗减小时在特定频率范围内会产生负向声辐射力. 这意味着可以通过充液球壳受到的声辐射力来确定内部液体的密度和声速等物理参数. 本文理论模型有助于分析基于声波对细胞核形态和位置异常细胞的筛选的物理机制, 对声操控技术实现病变组织进行监测和判断提供理论支持.

附录A

球坐标转换系数${Q_{mn}}$具体表达式如下^[27]:

$$\begin{split}& {Q_{mn}}({k_2}d) \\ =& \; {{\text{i}}^{m - n}}\sum\limits_{\mu = \left| {m - n} \right|}^{m + n} {{\varepsilon _{^\mu }}{{( - {\text{i}})}^\mu }(2\mu + 1)} b_m^{n\mu }{{\text{j}}_\mu }({k_2}d), \end{split}\tag{A1}$$

其中${\varepsilon _{^\mu }} = 1({\varepsilon _{^\mu }} = - {1^\mu })$当O₁在O₂的左(右)边, $b_m^{n\mu } = {(\left. {n\mu 00} \right|m0)^2}$, Clebsch-Gordan系数为$q = \left( {\mu + n + m} \right)/2$, 当2q是偶数时表达式为

$$\begin{split} & (\left. {n\mu 00} \right|m0) = \frac{{{{( - 1)}^{m + q}}q!}}{{\left( {q - n} \right)!\left( {q - \mu } \right)!\left( {q - m} \right)!}} \\&\times \sqrt {\frac{{(2m + 1)}}{{\left( {2q + 1} \right)!}}(2q - 2n)!(2q - 2\mu )!(2q - 2m)!} , \end{split} \tag{A2}$$

当2q是奇数时$(\left. {n\mu 00} \right|m0){\text{ = }}0$.

附录B

数值计算中使用介质的声速、密度及阻抗等参数见表B1, 弹性球壳的密度、杨氏模量和泊松比见表B2^[24].