臧雨宸 蘇暢3)? 吳鵬飛3) 林偉軍3)

1) (中國科學院聲學研究所,北京 100190)

2) (中國科學院大學,北京 100049)

3) (北京市海洋深部鉆探測量工程技術研究中心,北京 100190)

聲輻射力和聲輻射力矩的計算是實現粒子精準操控的重要基礎.基于經典聲散射理論的偏波級數展開法較難直接用于復雜模型的研究,而純數值的方法則不利于進行系統(tǒng)的參數化分析.基于Born 近似的基本原理,推導了低頻情況下零階Bessel 駐波場中心任意粒子的聲輻射力和力矩表達式.在此基礎上,以球形粒子、橢球形粒子和柱形粒子為例進行詳細地計算,并考慮聲參數的非均勻性對聲輻射力和力矩的影響.仿真結果表明,在低頻范圍內Born 近似具有很高的精度,隨著頻率的增加和粒子與流體的阻抗匹配變差,Born 近似的精度逐漸下降.對于傾斜放置于零階Bessel 駐波場中的橢球形粒子和柱形粒子,非對稱性會導致其受到聲輻射力矩的作用.在粒子尺寸遠小于波長的情況下,聲輻射力特性與粒子的具體形狀幾乎無關,但聲輻射力矩不然.最后,引入周圍流體的黏滯效應并對聲輻射力的表達式進行了修正.該研究預期可以為生物醫(yī)學、材料科學等領域利用駐波場聲鑷子實現微小粒子的精準操控提供一定的理論指導.

1 引言

當聲波在障礙物表面發(fā)生反射、散射、折射等物理效應時,會對物體產生力或力矩的作用,分別稱為聲輻射力和聲輻射力矩.基于聲輻射力和力矩的粒子精準操控技術在聲鑷子[1?6]、聲懸浮[7,8]等領域展現出巨大的應用前景,同時也對各種情況下聲輻射力和力矩的理論計算提出了更高的要求.傳統(tǒng)的聲輻射力和力矩研究主要基于經典的聲散射理論,利用偏波級數展開法求解散射系數進而計算聲輻射力和力矩,這一方法可以給出無窮級數形式的解析解,并在剛性球、剛性柱、彈性球、彈性柱等簡單模型中得到了成功運用[9?18].對于同心球殼、柱殼等簡單的非均勻粒子,偏波級數展開法依然適用,Hasegawa 等[19]和Mitri[20,21]最早給出了這一類問題的解析解,并詳細討論了球殼厚度和材料對聲輻射力的影響.Wang 等[22]計算了三層同心可壓縮球的聲輻射力從而來模擬對生物細胞的操控,隨后Peng 等[23,24]將其拓寬到聚焦聲束和聲表面波入射的情況.Gauss 駐波場、Bessel 駐波場等新型聲場對多層彈性球的聲輻射力特性也陸續(xù)為學者所研究[25?27].近年來,Mitri[28,29]還將該方法與柱函數的加法公式結合,成功研究了偏心柱的聲輻射力和力矩特性.在實驗探究方面,Mitri[30]通過實驗驗證了彈性柱在平面行波場中的解析解.Aglyamov等[31]使用1.5 MHz 的單陣元換能器對剛性小球施加聲輻射力并觀測小球的運動,其實驗測量結果和理論預測的結果符合得很好.Nikolaeva 等[32]通過實驗探究了傾斜入射的超聲束對聲吸收散射體的輻射力.Garbin 等[33]對盤狀粒子在聲輻射力和力矩作用下產生的聲電泳現象進行了詳細的數值和實驗研究,其結果為紅細胞等各種盤狀粒子的聲操控具有指導意義.Johnson 等[34]針對黏滯效應對平面行波場中聲輻射力的影響進行了實驗測量,并與已有的理論進行對比,結果顯示即使對于弱黏性流體,黏滯效應對聲輻射力的影響也十分顯著.Qiao等[35]對自由空間中球形粒子在黏性流體中的聲輻射力進行了理論和實驗研究,詳細討論了聲壓振幅、聲場頻率和流體黏度對輻射力的影響.

盡管如此,當粒子的形狀或結構更加復雜時,散射系數的解析計算變得較為困難甚至根本無法實現,此時有必要借助一些純數值的方法來進行求解.Wijaya 和Kim[36]利用邊界元方法成功研究了非球形粒子在駐波場中的聲輻射力和力矩,并以橢球體等形狀為例進行了計算.Glynne-Jones 等[37]則借助有限元方法分析了任意形狀的彈性和流體微粒的聲輻射力.誠然,數值方法的運用大大豐富了聲輻射力的研究模型,然而與解析方法相比,數值方法的計算量急劇增加,并且很難根據得到的結果進行模型的參數化分析.幸運的是,在某些情況下,尋求復雜模型的解析解或近似解是可能的.Wei 等[38]給出了聲場垂直入射時無限長橢圓柱在長波近似下的聲輻射力表達式,Hasheminejad 等[39]利用橢圓坐標系中的馬丟函數給出了適用于所有頻率范圍的結果.Marston[40]給出了低頻近似下長橢球與扁橢球在駐波場中的漸近表達式,其中要求橢球體的對稱軸與聲軸完全重合.Mitri[41?45]則沿用圓柱坐標系和球坐標系下的偏波級數展開法,借助廣義Fourier 級數展開理論和數值積分方法,利用半解析的方法成功計算了橢圓柱和橢球體的聲輻射力和力矩,其中要求橢圓柱和橢球的長短軸之比不能過大.Silva 等[46]研究了低頻近似下任意聲場中的橢球形粒子受到的聲輻射力矩,并以平面波為例進行了詳細分析.Jerome 等[47,48]則直接借助于橢球坐標系計算了可壓縮橢球體的聲輻射力和力矩.

Bessel 波束作為一類典型的非衍射波束,在光學和聲學領域都受到了廣泛關注.Marston[49,50]和Mitri[51?55]最早對Bessel 聲束的輻射力特性進行了詳細研究,特別關注了其產生的負向聲輻射力現象.對于高階Bessel 聲束而言,其攜帶的軌道角動量還可以對物體施加聲輻射力矩作用,Zhang 等[56?59],Gong 等[60]系統(tǒng)總結了這一方面的理論框架,并對Bessel 聲束聲輻射力矩的反轉現象進行了理論分析.近年來,Jerome 等[61?63]借鑒量子力學中的Born 近似方法避開了求解散射系數的繁復過程,直接利用曲面積分推導了平面駐波場作用下任意粒子的聲輻射力和力矩,并與精確解相比較,結果顯示在低頻范圍內兩者符合得很好.本文將Born近似方法拓展到Bessel 駐波場入射的情形,推導得到了零階Bessel 駐波場中任意粒子聲輻射力和力矩的一般積分表達式.在此基礎上,以球形粒子、橢球形粒子和柱形粒子為例進行了詳細分析,并考慮了粒子密度與聲速分布的非均勻性.此外,我們還引入背景流體的黏度,分析其對聲輻射力的影響.有必要指出,Born 近似方法的使用存在兩個前提條件:1)粒子的密度和聲速與周圍流體相差不能太大,2)是粒子處在低頻駐波場中.

2 聲輻射力和力矩的Born 近似

2.1 聲輻射力的一般表達式

考慮位于理想流體中的一半徑為R的球形粒子,以球心為原點建立球坐標系(r,θ,φ).根據波動方程和級數展開理論,任意的單頻入射聲場都可以展開為無窮級數的形式:

式中,p0是入射聲壓的復振幅,k=ω/c0是聲波在流體中的波數,ω和c0分別是聲波的角頻率和流體中的縱波聲速,jn和Ynm分別是n階球Bessel函數和n階球諧函數,anm是入射聲波的波束展開因子.為了簡便,我們略去了時間簡諧因子 e?iωt.當入射聲場與周向角φ無關時,(1)式可以進一步簡化為單重求和的級數:

式中,Pn是n階Legendre 函數,an是入射聲波的波束展開因子.散射聲場同樣也可以表示為級數展開的形式:

式中,ρ0是流體的密度,φinc和φsca分別是入射聲場和散射聲場對應的速度勢函數,Re 表示對物理量取實部,S是將粒子包圍在內的一個大封閉球面.將(2)式和(3)式帶入(4)式中,可以得到沿z方向的聲輻射力表達式:

有必要指出,這一表達式是聲輻射力的精確解,且適用于任意頻率的入射聲波.盡管如此,散射系數的求解通常是較為繁復的過程,因此有必要尋求進一步的化簡方法.

2.2 低頻駐波場聲輻射力的Born 近似

低頻駐波場是實際聲操控中的常見聲場,所謂低頻即滿足kR<1.低頻近似下,只需要考慮單極(n=0)和偶極(n=1)散射項即可,其余散射項均可以忽略.單極和偶極散射項對應的散射系數分別可以近似為[64?66]

式 中,ρm(r) 和cm(r) 分 別是粒 子的密 度與縱 波聲速.對于均勻粒子而言,其密度與聲速處處相等,因而f1(r) 和f2(r) 是與空間位置無關的常數;對于非均勻介質而言,f1(r) 和f2(r) 則是與空間位置有關的函數.將(6)式代入(5)式可以得到聲輻射力的低頻近似表達式為

考慮到球形粒子的體積為 4 πR3/3,(8)式可以進一步改寫為

當粒子的半徑很小時,球形粒子將退化為體積微元,(9)式變?yōu)樵擉w積微元受到的聲輻射力:

該式反映了聲輻射力和體積之間存在的線性關系,因而任意形狀粒子的聲輻射力都可以直接通過(10)式對作體積分得到,這也正是Born 近似方法的核心思想.需要注意的是,單個體積微元的散射聲場同樣會產生對其他體積微元的輻射力,因此嚴格而言通過直接對(10)式作體積分來計算聲輻射力是不準確的.因此,Born 近似方法僅適用于粒子的密度和聲速與周圍流體比較接近的情況,此時散射聲場的能量遠小于入射聲場,從而可以忽略不計.根據(7)式可以看出,此時 |f1(r)| 和 |f2(r)| 均遠小于1.對于生物體內的細胞與組織,這一條件是容易滿足的.

還有必要指出,Born 近似僅適用于駐波場中聲輻射力和聲輻射力矩的計算,對于行波場是不能運用的,這是由于行波場中體積微元的聲輻射力正比于dV的平方,從而不滿足線性疊加原理.事實上,駐波場與行波場產生聲輻射力的機理并不相同.前者不攜帶聲能流,其輻射力效應是由聲場動能和勢能的空間梯度引起;后者則攜帶聲能流,其輻射力效應則反映了聲波與物體之間的動量和能量傳遞.

以上主要針對聲輻射力進行了討論,聲輻射力矩的計算是類似的,只需將徑矢與其作矢量積運算即可.

2.3 零階Bessel 駐波場作用下的Born 近似

考慮位于零階Bessel 駐波場中的任意粒子.為了討論的方便,我們假定零階Bessel 駐波場的聲軸沿z方向且粒子恰好固定于聲束中心(坐標原點),此時聲壓可以表示為[67]

式 中,Jn是n階 柱Bessel 函 數,kr=ksinβ和kz=kcosβ分別是徑向和軸向的波矢分量,β是Bessel 聲束的半錐角,h是粒子中心到距離其最近的聲壓波節(jié)的距離.容易看出,當β=0 時,(11)式將退化為平面駐波場的表達式.(11)式與周向角無關,當然可以展開為(2)式的形式,其中的波束展開因子為[49]

將(12)式代入(10)式可以得到零階Bessel駐波場對體積微元的聲輻射力:

有必要指出,盡管我們只給出了軸向聲輻射力的結果,這并不意味著橫向聲輻射力 (Fx,Fy) 為零.但根據文獻[54]的結論,當滿足Born 近似的條件時,橫向聲輻射力比軸向聲輻射力小至少兩個數量級,因此只考慮軸向分量是合理的.

根據式(15),很容易通過矢量積運算得到粒子所受聲輻射力矩的表達式:

式中ez是z方向的單位矢量.

至此,在已知粒子形狀和密度與聲速空間分布的前提下,原則上可以根據(15)式和(16)式求得任意粒子聲輻射力和力矩的Born 近似解.對于密度與聲速均勻分布的粒子而言,(15)式和(16)式中的f1(r) 和f2(r)是與徑矢無關的常數,因而與之有關的項可以提到積分號外,此時聲輻射力和力矩的表達式簡化為:

對于密度與聲速均軸對稱分布的粒子,(15)式和(16)式亦可以進一步簡化.如圖1 所示,假定粒子的對稱軸位于xOz平面內,其對稱軸為z’軸,且與z軸所夾的角度為θs,以O為原點再建立一柱坐標系(ρ′,φ′,z′),兩坐標系之間存在以下關系:

圖1 傾斜放置于零階Bessel 駐波場中心的任意軸對稱粒子Fig.1.An arbitrary object with axisymmetric geometry obliquely positioned in a zero-order standing Bessel beam.

(22)式和(23)式中對于角度的積分可以直接利用柱Bessel 函數的性質得到解析解,于是聲輻射力和力矩的表達式分別簡化為:

式中J1表示一階柱Bessel 函數.進一步地,如果粒子的非均勻性僅僅在z’方向體現,即f1(r) 和f2(r) 僅僅是z′的函數,則(24)式和(25)式中關于極徑的積分可以直接進行計算,最終兩式變?yōu)槿缦碌膯沃胤e分:

當粒子尺寸遠小于聲波波長時(kR(z′)?1,kL ?1),還可以利用柱Bessel 函數的小宗量近似對(26)式和(27)式再次進行化簡,此時聲輻射力和力矩的近似表達式為:

3 球形粒子的Born 近似

3.1 均勻球形粒子

式中V=4πR3/3是球形粒子的體積.(31)式給出了聲輻射力矩為零的結果,對于均勻分布的球形粒子,這是可以預料的.

(30)式是通過Born 近似得到的結果,其精確性可以通過聲輻射力的精確表達式(5)式進行檢驗.為了方便比較,我們對(5)式和(30)式均用進行歸一化處理,這里S0=πR2是球形粒子的散射截面.圖2 同時給出了利用(30)式和(5)式對零階Bessel 駐波場作用下水中均勻球形粒子受到的歸一化聲輻射力的計算結果,水的密度ρ0=1000 kg/m3,聲速c0=1480 m/s ,β=π/6,kzh=π/4.假定粒子的密度與水相同,圖2(a),(b)和(c)分別對應著cm/c0=1.01 ,cm/c0=1.05和cm/c0=1.1 的情況.結果顯示,粒子受到的聲輻射力均可正可負,在kR<1 的低頻范圍內,Born近似的結果與精確解符合得很好,隨著kR的增加,Born 近似解的誤差開始增大.這是由于在中高頻范圍內,四極散射項(n=2)、八極散射項(n=3)等的影響不能忽略,同時這也再次驗證了Born 近似僅適用于kR<1 的低頻范圍.盡管如此,Born近似解和精確解在中高頻范圍內的變化趨勢基本保持一致,且聲輻射力和力矩極值所對應的kR基本不變.因此,在精度要求不是很高的情況下,Born近似方法可以為實際的聲操控提供一定的理論指導.此外,粒子與周圍流體的聲速相差越大,Born近似解和精確解的差異也越大.隨著兩者聲參數差異的增加,散射聲場對聲輻射力的貢獻不能再忽略,從而導致Born 近似方法誤差的增大.

圖2 零階Bessel 駐波場中心均勻球形粒子受到的歸一化聲輻射力隨kR 的變化(β=π/6,kzh=π/4,ρm/ρ0=1)(a) cm/c0=1.01;(b) cm/c0=1.05;(c) cm/c0=1.1Fig.2.The dimensionless acoustic radiation force plots for a homogeneous sphere versus kR in a zero-order standing Bessel beam (β=π/6,kzh=π/4,ρm/ρ0=1):(a) cm/c0=1.01;(b) cm/c0=1.05;(c) cm/c0=1.1.

圖3 給出了kR=0.5 時均勻球形粒子的歸一化聲輻射力隨聲束半錐角β的變化關系,粒子聲參數的設置與圖2 完全相同.容易看出,隨著粒子與周圍流體聲參數差異的增加,兩者的阻抗匹配變差,聲輻射力也隨之出現明顯的增強.此外,半錐角的增加均會導致聲輻射力明顯減小,當β=0 時,結果將退化為平面駐波場入射的情形.

圖3 零階Bessel 駐波場中心均勻球形粒子受到的歸一化聲輻射力隨β 的變化(kR=0.5,kzh=π/4,ρm/ρ0=1)Fig.3.The dimensionless acoustic radiation force plots for a homogeneous sphere versus β in a zero-order standing Bessel beam (kR=0.5,kzh=π/4,ρm/ρ0=1).

進一步考慮kR ?1 的情況.我們當然可以利用(28)式計算此時的聲輻射力,但直接計算(30)式在kR ?1 時的極限是更為簡便的方法.此時聲輻射力的公式簡化為:

3.2 非均勻球形粒子

實際的聲操控中粒子往往是非均勻的,即f1(r)和f2(r) 是空間位置的函數.考慮粒子密度和聲速僅僅隨徑向坐標r變化的情形,并滿足如下的關系:

容易看出,f1(r)+f2(r)/2 和f2(r) 在球心處分別為fA和fC,在球形粒子表面分別為fA+fB和fC+fD,并且隨徑向坐標線性變化.將(33)式和(34)式代入(15)式和(16)式,通過在球坐標中作體積分可以得到聲輻射力和力矩的表達式:

顯然,當fB=fD=0 時,(35)式將退化為均勻球形粒子的結果(30)式.此外,盡管此時粒子存在非均勻性,但對稱性使其仍然不受到聲輻射力矩的作用.

利用(35)式對水中非均勻球形粒子的聲輻射力進行分析.假設球心處的密度和聲速與水完全相同,此時顯然有fA=fC=0.圖3 計算了歸一化聲輻射力隨kR的變化關系,其中β=π/6 ,kzh=π/4.粒子表面的聲速和密度分別設置為四種情形:1420 kg/m3和940 m/s,1460 kg/m3和980 m/s,1500 kg/m3和1020 m/s,1540 kg/m3和1060 m/s.一般而言,密度越大的介質聲速越大,因此這樣的設置具有一定的合理性.經計算,在這4 種情況下分別有:fB=?0.250,fD=?0.042 ;fB=?0.077,fD=?0.014 ;fB=0.071,fD=0.013 ;fB=0.197,fD=0.039.結果顯示,粒子的聲參數與周圍流體相差越大,聲輻射力的峰值也越大,但聲輻射力峰值所對應的kR值幾乎保持不變.有必要指出,這是因為在Born 近似的過程中忽略了散射波對聲輻射力的貢獻.事實上,聲參數的改變必然會改變粒子的本征振動模式,進而影響聲輻射力峰值的位置.另一方面,如前所述,Born 近似解僅在kR<1時保持較高的精度,中高頻范圍內的結果參考意義有限.基于這一考慮,圖5 給出了kR=0.5 時歸一化聲輻射力隨半錐角的變化曲線,粒子聲參數的設置與圖4 完全相同.不難看出,隨著聲束半錐角的增加,聲輻射力的幅值明顯減小并最終趨于零.當β=0時聲輻射力的幅值最大,此時正對應著平面駐波場入射的情況.

圖4 零階Bessel 駐波場中心非均勻球形粒子受到的歸一化聲輻射力隨kR 的變化(fA=fC=0,β=π/6,kzh=π/4)Fig.4.The dimensionless acoustic radiation force plots for an inhomogeneous sphere versus kR in a zero-order standing Bessel beam (fA=fC=0,β=π/6,kzh=π/4).

圖5 零階Bessel 駐波場中心非均勻球形粒子受到的歸一化聲輻射力隨β 的變化(fA=fC=0,kR=0.5,kzh=π/4)Fig.5.The dimensionless acoustic radiation force plots for an inhomogeneous sphere versus β in a zero-order standing Bessel beam (fA=fC=0,kR=0.5,kzh=π/4).

同樣地,(35)式在kR ?1 的情況下還可以作進一步近似,其取極限的結果為

與(32)式類似,(37)式歸一化后的結果同樣與kR成正比(固定kzh).注意到當fB=fD=0 時,粒子將不存在非均勻性,(37)式將退化為(32)式.

從以上過程可以看出,對于由(33)式和(34)式表示的聲參數分布,聲輻射力的積分式是可以解析求解的,但對于更復雜分布的非均勻球形粒子,可能需要通過數值積分方法進行求解.

4 橢球形粒子的Born 近似

4.1 均勻橢球形粒子

橢球形粒子在實際的聲操控中也十分常見,如駐波場中的懸浮液滴在聲輻射力、重力和表面張力的耦合作用下通常呈現為橢球的形狀[8].考慮位于零階Bessel 駐波場中心的均勻橢球形粒子,其在對稱軸z’方向的半軸長度為b,垂直于對稱軸方向的橫截圓面的最大半徑為a,顯然a>b對應著扁橢球的情形,a

對于均勻橢球形粒子而言,其聲輻射力和力矩的積分表達式(17)式和(18)式都沒有顯式的解析結果,只能通過數值積分方法來進行計算.圖5 給出了零階Bessel 駐波場中心均勻橢球形粒子的歸一化聲輻射力和力矩隨kb的變化曲線,其中聲束的半錐角β=π/6 ,對稱軸的取向角θs=π/6,粒子的縱橫比分別為1/4,1/2,1,2 和4.粒子與水的密度和聲速之比分別為ρm/ρ0=1,cm/c0=1.05.與聲輻射力相對應,聲輻射力矩的歸一化因子為τ0=bF0.需要注意的是,在計算聲輻射力時設置kzh=π/4,計算聲輻射力矩時則設置kzh=0,這是為了讓粒子的聲輻射力和力矩效應均最為顯著.從圖6(a)不難看出,當增加粒子的縱橫比時,聲輻射力的峰值也隨之增加,且峰值所對應的kb值也隨之改變.盡管如此,所有曲線在kb ?1 的低頻范圍內幾乎完全重合,這說明在粒子尺寸遠小于波長時,聲輻射力特性與粒子的具體形狀幾乎無關.從圖6(b)可以發(fā)現,當粒子的對稱軸與聲軸既不平行也不垂直時,橢球形粒子會受到不為零的聲輻射力矩作用,且其方向會隨著kb的變化而變化,這預示著從y′軸正向看去,粒子既可以逆時針轉動也可以順時針轉動.當縱橫比為1 時,粒子將退化為球形,因而不存在聲輻射力矩.同時,橢球形粒子偏離球形的程度越大,聲輻射力矩的峰值也越大.注意到即使是在kb ?1 的范圍內,粒子縱橫比對聲輻射力矩特性的影響也十分顯著,這與聲輻射力的情形有所不同.

圖6 零階Bessel 駐波場中心均勻橢球形粒子受到的歸一化聲輻射力和力矩隨kb 的變化(β=π/6,θs=π/6,ρm/ρ0=1,cm/c0=1.05) (a)歸一化聲輻射力(kzh=π/4);(b)歸一化聲輻射力矩(kzh=0)Fig.6.The dimensionless acoustic radiation force and torque plots for a homogeneous spheroid versus kb in a zeroorder standing Bessel beam (β=π/6,θs=π/6,ρm/ρ0=1,cm/c0=1.05):(a) Dimensionless acoustic radiation force(kzh=π/4);(b) dimensionless acoustic radiation torque(kzh=0).

圖7(a)和圖7(b)分別給出了kb=0.5 時歸一化聲輻射力和力矩隨粒子對稱軸取向角的變化關系,其余參數設置與圖6 完全相同.從仿真結果可以看出,橢球形粒子的聲輻射力關于θs=π/2 呈偶對稱,這是由模型的對稱性決定的.在θs=π/2 處,長橢球形粒子的聲輻射力達到最大值,而扁橢球形粒子的聲輻射力達到最小值.對于球形粒子而言,聲輻射力曲線退化為與θs無關的一條水平直線.與聲輻射力的結果不同,橢球形粒子的聲輻射力矩關于θs=π/2 呈奇對稱,在θs=0,π/2,π 處力矩均恒為零,這也是模型對稱性所要求的必然結果.在0<θs<π/2的范圍內,長橢球形粒子的聲輻射力矩取負值,扁橢球形粒子的聲輻射力取正值,在π/2<θs<π的范圍內則恰好相反.此外,球形粒子的聲輻射力矩恒為零,這是與預期相符的.

圖7 零階Bessel 駐波場中心均勻橢球形粒子受到的歸一化聲輻射力和力矩隨θs 的變化(kb=0.5,β=π/6,ρm/ρ0=1,cm/c0=1.05) (a)歸一化聲輻射力(kzh=π/4);(b)歸一化聲輻射力矩(kzh=0)Fig.7.The dimensionless acoustic radiation force and torque plots for a homogeneous spheroid versus θs in a zeroorder standing Bessel beam (kb=0.5,β=π/6,ρm/ρ0=1,cm/c0=1.05):(a) Dimensionless acoustic radiation force(kzh=π/4);(b) dimensionless acoustic radiation torque(kzh=0).

在kb ?1 的條件下,可以利用(28)式和(29)式來計算此時聲輻射力和力矩的近似表達式:

式中V=4πa2b/3 是橢球形粒子的體積.(38)式形式上與(32)式完全相同,這再次驗證了在粒子尺寸遠小于波長的情況下,聲輻射力特性與粒子的具體形狀和粒子對稱軸的取向角θs均無關.(39)式在一般情況下是非零的,且依賴于粒子對稱軸的取向角θs,當粒子的對稱軸與聲軸平行或垂直時力矩為零.容易看出,當a=b時聲輻射力矩也將消失,這是因為此時粒子退化為球形.有必要指出,在固定kzh的條件下,歸一化后的聲輻射力矩并非與kb成正比,而是與成正比,這與聲輻射力的情形有所不同.

4.2 非均勻橢球形粒子

嘗試將非均勻性引入橢球形粒子.考慮粒子的密度和聲速僅僅隨z′坐標變化的情形,并滿足如下的關系:

由于此時的非均勻性僅僅體現在z′坐標上,因此可以通過(26)式和(27)式對橢球形粒子的聲輻射力和力矩進行計算.與均勻橢球形粒子的情況類似,此時仍然無法得到積分的解析結果,需要進行數值積分運算.圖8 計算了零階Bessel 駐波場中心非均勻橢球形粒子的歸一化聲輻射力隨kb的變化關系,圖8(a)中的橢球形粒子在z′=?L/2 處ρm=1000 kg/m3,cm=1480 m/s,在z′=L/2 處ρm=1080 kg/m3,cm=1560 m/s,圖8(b)中的橢球形粒子在z′=?L/2 處m=1000 kg/m3,cm=1480 m/s,在z′=L/2 處m=920 kg/m3,cm=1400 m/s,經 計算在這兩種情況下分別有:fA=0.137,fB=0.254,fC=0.026,fD=0.051 ;fA=?0.160,fB=?0.349,fC=?0.027,fD=?0.056,其余參數設置與圖6完全相同.結果顯示,圖8(a)的變化規(guī)律與圖6(a)基本一致,但聲輻射力的幅值更大.圖8(b)的仿真結果則與圖8(a)異號,這是由于此時粒子的聲阻抗小于周圍流體.圖9 給出了與圖8 相同情況下歸一化聲輻射力矩的計算結果,其中圖9(a)和圖6(b)的變化規(guī)律幾乎一致,而圖9(b)則與圖9(a)異號.

圖8 零階Bessel 駐波場中心非均勻橢球形粒子受到的歸一化聲輻射力隨kb 的變化(β=π/6,θs=π/6,kzh=π/4):(a) fA=0.137,fB=0.254,fC=0.026,fD=0.051;(b) fA=–0.160,fB=–0.349,fC=–0.027,fD=–0.056Fig.8.The dimensionless acoustic radiation force plots for an inhomogeneous spheroid versus kb in a zero-order standing Bessel beam (β=π/6,θs=π/6,kzh=π/4):(a) fA=0.137,fB=0.254,fC=0.026,fD=0.051;(b) fA=–0.160,fB=–0.349,fC=–0.027,fD=–0.056.

圖9 零階Bessel 駐波場中心非均勻橢球形粒子受到的歸一化聲輻射力矩隨kb 的變化(β=π/6,θs=π/6,kzh=0) (a) fA=0.137,fB=0.254,fC=0.026,fD=0.051;(b) fA=–0.160,fB=–0.349,fC=–0.027,fD=–0.056Fig.9.The dimensionless acoustic radiation torque plots for an inhomogeneous spheroid versus kb in a zero-order standing Bessel beam (β=π/6,θs=π/6,kzh=0):(a) fA=0.137,fB=0.254,fC=0.026,fD=0.051;(b) fA=–0.160,fB=–0.349,fC=–0.027,fD=–0.056.

圖10 和圖11 給出了kb=0.5 時歸一化聲輻射力和力矩隨粒子對稱軸取向角的變化關系,其余參數的設置與圖8 和圖9 完全相同.同樣地,圖10(a)和圖11(a)的計算結果與圖7(a)和圖7(b) 完全類似,圖10(b)和圖11(b)的結果則分別與圖10(a)和圖11(a)異號.

圖10 零階Bessel 駐波場中心非均勻橢球形粒子受到的歸一化聲輻射力隨θs 的變化(kb=0.5,β=π/6,kzh=π/4)(a) fA=0.137,fB=0.254,fC=0.026,fD=0.051;(b) fA=–0.160,fB=–0.349,fC=–0.027,fD=–0.056Fig.10.The dimensionless acoustic radiation force plots for an inhomogeneous spheroid versus θs in a zero-order standing Bessel beam (kb=0.5,β=π/6,kzh=π/4):(a) fA=0.137,fB=0.254,fC=0.026,fD=0.051;(b) fA=–0.160,fB=–0.349,fC=–0.027,fD=–0.056.

圖11 零階Bessel 駐波場中心非均勻橢球形粒子受到的歸一化聲輻射力矩隨θs 的變化(kb=0.5,β=π/6,kzh=0) (a) fA=0.137,fB=0.254,fC=0.026,fD=0.051;(b) fA=–0.160,fB=–0.349,fC=–0.027,fD=–0.056Fig.11.The dimensionless acoustic radiation torque plots for an inhomogeneous spheroid versus θs in a zero-order standing Bessel beam (kb=0.5,β=π/6,kzh=0):(a) fA=0.137,fB=0.254,fC=0.026,fD=0.051;(b) fA=–0.160,fB=–0.349,fC=–0.027,fD=–0.056.

考察kb ?1 且縱橫比b/a ?1 的情況.利用(28)式和(29)式得到此時聲輻射力和力矩的近似表達式:

可以看出,聲輻射力和力矩均與粒子對稱軸的取向角有關.若fB和fD均取正值,則當θs=0 時粒子的聲輻射力最大,當θs=π/2 時粒子的聲輻射力最小,但這兩種情況下聲輻射力矩均由于對稱性而消失.考慮到kb ?1 ,表征粒子聲參數非均勻性的fB和fD對聲輻射力的貢獻遠小于粒子聲參數的平均值fA和fC,而聲輻射力矩的情況則恰好相反.值得注意的是,在固定kzh的前提下,此時的歸一化聲輻射力和力矩與頻率的依賴關系較為復雜,兩者不再存在簡單的冪次關系.對于均勻橢球形粒子而言,只需置fB=fD=0 即可,此時(42)式退化為(38)式,但(43)式并不退化為(39)式.若給(39)式也加上縱橫比b/a ?1 的前提,則兩者完全等價,這也是預期的結果.

5 柱形粒子的Born 近似

5.1 均勻柱形粒子

Born 近似同樣也適用于柱形粒子.考慮截面直徑為D、長度為L的一均勻柱形粒子,根據幾何關系顯然有R(z)=D/2,此時(26)式和(27)式存在解析的積分結果,聲輻射力和力矩分別為:這里V=(π/4)D2L是柱形粒子的體積.從(45)式可以看出,當柱形粒子的對稱軸平行或垂直于聲軸時,聲輻射力矩將由于對稱性消失,這是預料之中的結果.

圖12 給出了零階Bessel 駐波場中心均勻柱形粒子的歸一化聲輻射力和力矩隨kL的變化曲線,其中聲束的半錐角β=π/6,對稱軸的取向角θs=π/6,D/L分別為1/2,1 和2.粒子與水的密度和聲速之比仍然設置為ρm/ρ0=1,cm/c0=1.05.對于柱形粒子,此時聲輻射力和聲輻射力矩的歸一化因子分別為,其中S0=(π/4)D2是柱形粒子的橫截面積.同樣地,在計算聲輻射力時設置kzh=π/4,計算聲輻射力矩時則設置kzh=0.結果顯示,不同的聲輻射力曲線在kL ?1 的低頻范圍內幾乎完全重合,這同樣反映了粒子尺寸遠小于波長時聲輻射力特性和粒子的具體形狀無關,而聲輻射力矩則無此性質.聲輻射力和力矩的方向均同時受到kL和D/L的影響,可取正值或負值.

圖12 零階Bessel 駐波場中心均勻柱形粒子受到的歸一化聲輻射力和力矩隨kL 的變化(β=π/6,θs=π/6,ρm/ρ0=1,cm/c0=1.05) (a)歸一化聲輻射力(kzh=π/4);(b)歸一化聲輻射力矩(kzh=0)Fig.12.The dimensionless acoustic radiation force and torque plots for a homogeneous cylinder versus kL in a zero-order standing Bessel beam (β=π/6,θs=π/6,ρm/ρ0=1,cm/c0=1.05):(a) Dimensionless acoustic radiation force (kzh=π/4);(b) dimensionless acoustic radiation torque (kzh=0).

圖13 給出了kL=0.5 時歸一化聲輻射力和力矩隨粒子對稱軸取向角的變化關系,其余參數設置與圖12 完全相同.與橢球形粒子的結果類似,柱形 粒子的聲輻射力 關于θs=π/2 偶對 稱,且 在θs=π/2處取得極大值或極小值.聲輻射力矩則關于θs=π/2 奇對稱,且在θs=0,π/2,π 力矩均由 于對稱性消失.

圖13 零階Bessel 駐波場中心均勻柱形粒子受到的歸一化聲輻射力和力矩隨θs 的變化(β=π/6,kL=0.5,ρm/ρ0=1,cm/c0=1.05) (a)歸一化聲輻射力(kzh=π/4);(b)歸一化聲輻射力矩(kzh=0)Fig.13.The dimensionless acoustic radiation force and torque plots for a homogeneous cylinder versus θs in a zeroorder standing Bessel beam (β=π/6,kL=0.5,ρm/ρ0=1,cm/c0=1.05):(a) Dimensionless acoustic radiation force (kzh=π/4);(b) dimensionless acoustic radiation torque (kzh=0).

考慮橫截面半徑和粒子長度均遠小于波長的情形(kD ?1,kL ?1),此時的聲輻射力和力矩可以近似表示為更簡潔的形式:

(46)式和(32)式形式上完全相同,再次驗證了粒子尺寸遠小于波長時聲輻射力特性和粒子的具體形狀無關.(47)式則與D和L均有關.有趣的是,當時,無論對稱軸的取向角如何,粒子的聲輻射力矩均恒為零,這對于實際聲操控中如何保持粒子的力矩平衡有重要的理論指導意義.

5.2 非均勻柱形粒子

嘗試在柱形粒子模型中引入非均勻性.同樣考慮粒子的密度和聲速僅僅隨z’坐標變化的情形,并滿足(40)式和(41)式所示的線性關系.顯然,粒子的聲參數仍然隨徑向坐標呈線性變化,f1(z′)+f2(z′)/2和f2(z′) 在柱形粒子的一端z′=?L/2 分別取值為fA?fB和fC?fD,在另一端z′=L/2 分別取值為fA+fB和fC+fD.將(40)式和(41)式代入(26)式和(27)式可以得到此時聲輻射力和力矩的表達式:

根據(48)式和(49)式,可以直接對非均勻柱形粒子的聲輻射力和力矩進行仿真,這里不再詳細展開.

進一步考察kD ?1,kL ?1 且D ?L的情況,(48)式和(49)式可以近似表示為:

同樣地,考慮到kL ?1,表征粒子非均勻性的參數fB和fD對聲輻射力的貢獻遠小于粒子聲參數的平均值fA和fC,而聲輻射力矩的情況則恰好相反.對于均勻柱形粒子(fB=fD=0),(50)式將退化為(46)式,(51)但式并不退化為(47)式.類似地,若給(51)式也加上D ?L的前提條件,則(51)式和(47)式完全一致.

6 周圍流體黏滯效應的影響

在實際的聲操控中,周圍流體往往存在一定的黏滯效應,該效應不可避免地會對聲輻射力產生一定的影響,因而需要對理想流體中的結果作修正.根據流體力學的基本理論,黏性流體中聲波的邊界層厚度定義為[68,69]

式中,ν=η/ρ0是流體的動力學黏滯系數,η是流體本身的黏滯系數.對于室溫下水中1 kHz 的聲波而言,聲波的邊界層厚度約為0.6 μm.對于半徑為R的球形粒子,定義歸一化邊界層厚度為

根據文獻[70]中的討論,在考慮流體黏滯吸收的情況下,單極散射系數f1(r) 保持不變,偶極散射系數f2(r) 則會變?yōu)槿缦碌膹碗s形式:

一般而言,(54)式是復數量.這里我們考慮兩種特殊的情況:對于密度與周圍流體相等的粒子,有f1(r)≡0,此時的偶極散射系數與流體的黏滯效應無關.對于粒子尺寸較小且流體黏度較大的情況(?1),此時的偶極散射系數近似為一個實數量f2(r)≈2[ρm(r)?ρ0(r)]/(3ρ0(r)).

對于均勻球形粒子,(54)式與空間位置無關.圖14(a)和圖14(b) 分別給出了此時偶極散射系數f2的實部和虛部隨歸一化邊界層厚度的變化曲線,其中粒子與周圍流體的密度比分別設置為ρm/ρ0=0.9, 0.95, 1, 1.05, 1.1.為了便于比較,圖14(a)中各曲線均作了歸一化處理,歸一化因子為理想流體(=0)中偶極散射系數的實部(記為 R e(f20)).從圖中可以看出,隨著粒子密度的增加,R e(f2) 也隨之增加.當增加邊界層厚度時,R e(f2) 在ρm/ρ0>1時先增加后趨于穩(wěn)定,而在ρm/ρ0<1 時先減小后趨于穩(wěn)定.與之有所不同的是,無論粒子的密度如何(ρm/ρ0=1 除外),I m(f2) 均隨著邊界層厚度的增大先增加而后減小,在≈0.5 時達到最大值.此外,粒子密度與流體密度越接近,I m(f2) 越小.

圖14 均勻球形粒子的偶極散射系數f2 隨 的變化 (a)Re(f2)/Re(f20) ;(b) Im(f2)Fig.14.The dipole scattering coefficient f2 plots for a homogeneous sphere versus (a)Re(f2)/Re(f20);(b) Im(f2).

接下來考慮對聲輻射力的修正.具體地,只需要用 R e(f2(r)) 代 替f2(r) 即可.于是理想流體中任意聲場入射下的體積微元受到的聲輻射力(10)式變?yōu)?/p>

對(56)式作體積分即得到黏性流體中球形粒子的聲輻射力的一般表達式.

考慮零階Bessel 駐波場中心的均勻球形粒子,此時需要將理想流體中的結果(30)式修正為

圖15(a)計算了不同歸一化邊界層厚度下水中均勻球形粒子的歸一化聲輻射力隨kR的變化曲線,其中β=π/6,k z h=π/4,粒子與流體的聲參數之比為ρm/ρ0=1.2,c m/c0=1.1.為了進一步凸顯流體黏滯效應對聲輻射力的影響,圖15(b)給出了不同歸一化邊界層厚度下黏性流體中聲輻射力和理想流體中歸一化聲輻射力的差值曲線,其參數設置與圖15(a)完全相同.計算結果顯示,由于流體黏滯效應的存在,粒子的正向與負向聲輻射力均會得到一定的增強,且黏滯效應的影響隨聲輻射力幅值的增大而增大.當邊界層厚度達到數倍波長時,由于 R e(f2(r)) 趨于穩(wěn)定,聲輻射力也不再隨邊界層厚度的增加而變化.整體而言,流體黏滯效應對聲輻射力的影響較為微弱,因而理想流體中的結果不失為很好的近似.必須指出,這一結論是在粒子與流體密度差異不大的前提下得出的,若兩者的密度相差明顯,黏滯效應的影響可能較為顯著,但此時的模型已經不再滿足Born 近似的適用條件,故而這里不再詳細討論.

圖15 零階Bessel 駐波場中心均勻球形粒子受到的歸一化聲輻射力隨kR 的變化(β=π/6,kzh=π/4,ρm/ρ0=1.2,cm/c0=1.1) (a) 歸一化聲輻射力;(b) 黏性流體與理想流體中歸一化聲輻射力的差值Fig.15.The dimensionless acoustic radiation force plots for a homogeneous sphere versus kR in a zero-order standing Bessel beam (β=π/6,kzh=π/4,ρm/ρ0=1.1,cm/c0=1.1):(a) Dimensionless acoustic radiation force;(b) difference of dimensionless acoustic radiation force in a viscous fluid and in an ideal fluid.

7 結論

本文將Born 近似方法運用到低頻零階Bessel駐波場中,推導了位于聲場中心的任意粒子受到的聲輻射力和力矩的積分表達式,并以球形粒子、橢球形粒子和柱形粒子為例進行大量的仿真計算.主要得到如下結論:

1)在kR<1 的低頻范圍內,Born 近似的結果與精確解符合得很好.隨著kR的增加,Born近似解的誤差逐漸增大.此外,粒子與周圍流體聲參數的差異越大,Born 近似解的誤差也越大;

2)當橢球形粒子或柱形粒子的對稱軸不與聲軸平行或垂直時,會產生由于非對稱性導致的聲輻射力矩作用.橢球形粒子、柱形粒子的聲輻射力和力矩關于θs=π/2 分別呈偶對稱和奇對稱;

3)當粒子的尺寸遠小于聲波波長時,粒子的具體形狀對其聲輻射力特性幾乎沒有影響,而聲輻射力矩特性則受到縱橫比的顯著影響;

4)周圍流體的黏滯效應主要通過影響偶極散射系數進而影響聲輻射力,當粒子與周圍流體密度差異不大時,黏滯效應的影響較為微弱.

有必要再次強調,利用Born 近似的方法求解聲輻射力和力矩時必須滿足兩個前提條件:1)粒子與所處流體介質的聲參數相差不能太大;2)粒子所處的聲場必須低頻駐波場.本文的研究結果提供了一種求解駐波場聲輻射力的新思路,可以為實現粒子的精準聲操控和無容器處理技術提供必要的理論基礎.