孫秀娜, 張小鳳, 張光斌，汪 艷, 常國(guó)棟

(陜西師范大學(xué) 物理學(xué)與信息技術(shù)學(xué)院，陜西省超聲學(xué)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710119)

FDTD法計(jì)算剛性球形粒子離軸聲輻射力

孫秀娜, 張小鳳, 張光斌*，汪 艷, 常國(guó)棟

(陜西師范大學(xué) 物理學(xué)與信息技術(shù)學(xué)院，陜西省超聲學(xué)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710119)

為了更精確、簡(jiǎn)便計(jì)算物體的聲輻射力，應(yīng)用時(shí)域有限差分算法對(duì)聲波波動(dòng)方程在時(shí)間和空間上進(jìn)行差分離散化，通過(guò)設(shè)置有效的邊界條件，完成聲壓場(chǎng)和速度場(chǎng)在時(shí)域上的迭代。將時(shí)域有限差分法與應(yīng)力張量方程相結(jié)合，求出了二維情況下水中剛性球形粒子的聲輻射力，研究了聲輻射力隨頻率和粒子偏離高斯波束中心距離的變化關(guān)系。仿真結(jié)果表明，剛性球形粒子在高斯聲場(chǎng)中的聲輻射力隨著粒子軸向偏離波束中心距離的增大而減小。當(dāng)剛性球形粒子偏離波束軸時(shí)，軸向聲輻射力與軸向、橫向的偏移距離成反比。橫向聲輻射力隨橫向偏移距離的增加先增大后減小，在偏移波束中心兩側(cè)的距離相同時(shí)，輻射力的大小相等，方向相反。

時(shí)域有限差分法；高斯波束；離軸聲輻射力；數(shù)值仿真

PACS: 43.35.+d

聲輻射力因其無(wú)損、非接觸等特性被廣泛地應(yīng)用到聲懸浮、聲學(xué)粒子操控、超聲給藥等各個(gè)領(lǐng)域。近年來(lái),隨著微操控技術(shù)的快速發(fā)展,利用聲輻射力操控微納米顆粒成為物理、生物、工程等領(lǐng)域的熱門(mén)研究方向。

聲輻射力的研究，開(kāi)始于上世紀(jì)30年代。1934年，King首次計(jì)算了駐波和行波聲場(chǎng)中剛性小球的聲輻射力，給出了理想流體中剛性球聲輻射力的解析表達(dá)式，該研究奠定了聲操控技術(shù)理論基礎(chǔ)。隨后，各國(guó)學(xué)者就聲輻射力的理論進(jìn)行了不懈的研究。目前，有關(guān)聲輻射力的理論研究，主要包括米氏散射法、聲線法和時(shí)域有限差分法。米氏散射法是利用解析理論，通過(guò)求解波動(dòng)方程計(jì)算聲場(chǎng)，進(jìn)而求解聲輻射力，但是這種方法只能對(duì)球形粒子或者無(wú)限長(zhǎng)的圓柱等規(guī)則形狀物體進(jìn)行求解。文獻(xiàn)[1-3]利用米氏散射法研究了規(guī)則聲場(chǎng)下粒子的聲輻射力特性。聲線法是聲波的波長(zhǎng)相對(duì)散射體尺寸非常小，衍射被忽略時(shí)，將聲波看成類(lèi)似光線的聲線，用聲線的疏密來(lái)表示聲場(chǎng)的強(qiáng)弱，進(jìn)而來(lái)計(jì)算聲輻射力的方法，只適用物體尺寸遠(yuǎn)大于波長(zhǎng)的情況。文獻(xiàn)[4-5]采用聲線理論計(jì)算了聲場(chǎng)中任意位置流體小球的聲輻射力，給出了聲輻射力與粒子捕獲能力的關(guān)系。

時(shí)域有限差分法(Finite Difference Time Domain Method ,FDTD)被廣泛地應(yīng)用于電磁波傳播和散射的數(shù)值仿真中[6]。它通過(guò)直接離散Maxwell方程計(jì)算電磁場(chǎng)，不需要任何形式的導(dǎo)出方程，故不會(huì)因?yàn)閿?shù)學(xué)模型而限制其應(yīng)用范圍，其差分方程中直接包含介質(zhì)的參量，只須賦予各網(wǎng)格相應(yīng)的參量，就能模擬各種散射體的結(jié)構(gòu)。目前，基于時(shí)域有限差分法與麥克斯韋應(yīng)力張量方程相結(jié)合的方法計(jì)算光俘獲力的準(zhǔn)確性已經(jīng)得到了多次驗(yàn)證。 2010年，蔡飛燕、程欣等首次運(yùn)用時(shí)域有限差分法計(jì)算了聲場(chǎng)中的輻射力[7-8]。本文應(yīng)用時(shí)域有限差分法，結(jié)合應(yīng)力張量方程，計(jì)算了高斯波束作用下水中剛性球形粒子在聲場(chǎng)中的不同位置時(shí)的聲輻射力與軸向、橫向偏移距離的變化關(guān)系。計(jì)算結(jié)果與利用米氏散射法計(jì)算得出的聲輻射力非常一致；并且時(shí)域有限差分法不需要借用數(shù)學(xué)工具來(lái)推導(dǎo)大量公式，對(duì)物體的尺寸、形狀、物理性質(zhì)沒(méi)有特殊的限定。

1 基于FDTD的聲輻射力計(jì)算

1.1 聲波動(dòng)方程的離散化

在均勻的流體介質(zhì)中，二維直角坐標(biāo)下的聲波方程為

(1)

(2)

(3)

式中：c為介質(zhì)聲速，ρ為介質(zhì)密度，vx、vy為質(zhì)點(diǎn)在x和y方向的振動(dòng)速度，p為質(zhì)點(diǎn)聲壓。

用FDTD仿真聲場(chǎng)時(shí)，首先將(1)—(3)式的微分方程轉(zhuǎn)化為離散差分方程，在時(shí)間軸上逐步推進(jìn)求解。 由于計(jì)算區(qū)域是有限的，聲波在計(jì)算區(qū)域的邊界處會(huì)發(fā)生反射，影響聲波的正常傳播。為了減小在邊界處反射波的影響，選擇完全匹配層吸收邊界條件(Perfectly Matched Layer，PML)[9]。為使聲波在PML中快速衰減，需要對(duì)聲波波動(dòng)方程進(jìn)行指數(shù)差分變換。為了便于計(jì)算，將(1)—(3)式中的聲壓p分解為兩個(gè)子分量px和py，引入衰減系數(shù)αx和αy，使聲壓和速度在PML層中以指數(shù)衰減。完全匹配層吸收邊界條件內(nèi)部的指數(shù)差分公式為

(4)

(5)

(6)

(7)

1.2 應(yīng)用FDTD計(jì)算聲輻射力

聲輻射力的求解是對(duì)物體表面的應(yīng)力張量進(jìn)行迭代積分，即

(8)

(9)

在時(shí)域有限差分算法中，聲壓和速度分量在時(shí)間和空間上是交錯(cuò)分布的。求解聲輻射力時(shí)，需要計(jì)算的是物體表面上分布的網(wǎng)格點(diǎn)處的應(yīng)力張量。由于聲場(chǎng)是不斷變化的，計(jì)算聲輻射力時(shí)，對(duì)瞬時(shí)輻射力取一個(gè)周期進(jìn)行時(shí)間平均，即

(10)

式中T為入射激勵(lì)源的周期。

為了方便計(jì)算，利用傅里葉變換將聲壓從時(shí)域變換到頻域，即

(11)

式中ω為入射波的角頻率,N為時(shí)間迭代次數(shù)。當(dāng)對(duì)所有變量做傅里葉變換后，每個(gè)Yee網(wǎng)格的法向和切向力的表達(dá)式[12]可以分別表示為

(12)

(13)

(14)

I=∫〈p2〉/(ρc2)dl。

(15)

其中l(wèi)是繞物體表面的長(zhǎng)度。

2 數(shù)值仿真

2.1 邊界仿真

為了驗(yàn)證加入PML吸收邊界條件的時(shí)域有限差分法的計(jì)算區(qū)域滿足無(wú)限大區(qū)域聲場(chǎng)的仿真要求，在Matlab軟件中對(duì)一個(gè)連續(xù)上升余弦激勵(lì)源在水中的聲場(chǎng)進(jìn)行了仿真。仿真時(shí)的參數(shù)設(shè)置為：二維計(jì)算區(qū)域400×400網(wǎng)格，網(wǎng)格間隔Δx=Δy=λmin/20，其中λmin=c/fmax是最小波長(zhǎng)，最大頻率fmax=3 750 Hz，時(shí)間步長(zhǎng)Δt=Δx/2c，聲速c=1 500 m/s，完全匹配層個(gè)數(shù)h=10，衰減因子m=1.8。圖1和圖2分別為未加邊界和加入PML吸收邊界時(shí)1 000時(shí)間步的聲場(chǎng)仿真結(jié)果。由圖1和圖2的比較可以看出，未加邊界時(shí)，聲波傳到邊界被反射回來(lái)，入射波和反射波在計(jì)算區(qū)域相互疊加。加PML吸收邊界后，聲波被匹配層吸收，幾乎看不到反射波，有效解決了用有限區(qū)域模擬無(wú)限大區(qū)域聲場(chǎng)的問(wèn)題。

圖1 未加邊界條件的聲波傳播Fig.1 Acoustic wave propagation without boundary conditions

圖2 加PML邊界條件的聲波傳播Fig.2 Acoustic wave propagation with PML boundary condition

2.2 調(diào)制高斯脈沖波源的仿真

設(shè)聲源為調(diào)制高斯脈沖波源，其聲壓表達(dá)式為

(17)

其中:x0是高斯波波束中心的坐標(biāo)，w是高斯波束的束腰半徑。ω是調(diào)制高斯脈沖波的角頻率，τ是脈沖寬度。仿真時(shí)的參數(shù)設(shè)置為：x0=200網(wǎng)格處，角頻率ω=24 000，脈沖寬度τ=9π/(2w)，其他參數(shù)設(shè)置同2.1節(jié)。利用FDTD方法，分別對(duì)束腰半徑為w=0.01λ和w=λ時(shí)，200時(shí)間步的高斯聲場(chǎng)進(jìn)行了仿真，結(jié)果如圖3和圖4所示。

圖3 高斯波束束腰半徑w=0.01λ時(shí)聲場(chǎng)的分布Fig.3 The acoustic distribution field of Gaussianbeam with beam width w=0.01λ

圖4 高斯波束束腰半徑w=λ時(shí)聲場(chǎng)的分布Fig.4 The acoustic distribution field of Gaussian beam with beam width w=λ

從圖3和圖4的仿真結(jié)果可以看出，高斯波束的束腰半徑不同，對(duì)應(yīng)的聲場(chǎng)空間分布不同。束腰半徑相對(duì)于波長(zhǎng)較小時(shí)，高斯波的傳播規(guī)律接近于點(diǎn)源的傳播。當(dāng)束腰半徑等于1倍波長(zhǎng)時(shí)，高斯波束相當(dāng)于一條線形波源，聲波強(qiáng)度在空間上呈現(xiàn)高斯分布。

2.3 聲輻射力的數(shù)值仿真

2.3.1 聲輻射力隨頻率的變化 對(duì)水中剛性球形粒子在高斯波束作用下的離軸聲輻射力進(jìn)行了仿真。仿真時(shí)的參數(shù)設(shè)置為：計(jì)算區(qū)域1 000×1 000網(wǎng)格，聲速c=1 500 m/s，水的密度ρ=1×103kg/m3，剛性球形粒子的半徑a=40Δx。調(diào)制高斯脈沖波的束腰半徑取1倍波長(zhǎng)。高斯波束的中心位置位于x0=500,y0=30網(wǎng)格處，球形粒子位于x=500,y=110網(wǎng)格處。根據(jù)(14)式，分別計(jì)算出剛性球形粒子的軸向和橫向聲輻射力隨無(wú)量綱頻率ka的變化關(guān)系曲線,如圖5、圖6所示。

圖5 軸向聲輻射力隨ka的變化關(guān)系Fig.5 The relationship of axial acoustic radiation force and ka

圖6 橫向聲輻射力隨ka的變化關(guān)系Fig.6 The relationship of transverse radiation force and ka

從圖5的仿真結(jié)果可知，在ka接近于0的位置處，軸向聲輻射力是一個(gè)負(fù)值，隨ka的增大逐漸由負(fù)值趨向于0逐漸變化為正輻射力，并隨ka的增大而迅速增大，達(dá)到一個(gè)最大值。在ka<1時(shí)，軸向聲輻射力隨ka的增加有一個(gè)劇烈的波動(dòng)，從負(fù)方向上的最大值變化到正方向上的最大值。不僅力的大小迅速變化，力的方向也發(fā)生改變。在12時(shí)，剛性球形粒子的軸向聲輻射力隨ka的增加變化幅度不大，基本不受ka的影響。仿真結(jié)果與文獻(xiàn)[10-12]中運(yùn)用解析方法計(jì)算的聲輻射力隨頻率的變化規(guī)律一致。圖6是剛性球形粒子橫向聲輻射力隨ka的變化關(guān)系曲線，從圖中的計(jì)算結(jié)果可見(jiàn)，在02時(shí)，聲輻射力是一個(gè)趨于零的值。這是由于高斯波束的對(duì)稱性，當(dāng)球形粒子位于聲波傳播軸上時(shí)，球形粒子橫向聲輻射力的合力為0。由圖5和圖6 的仿真結(jié)果可見(jiàn)，在頻率很小(ka<2)時(shí)，頻率對(duì)聲輻射力的影響很大，隨著頻率的變化，輻射力的大小和方向都發(fā)生很大變化；當(dāng)ka>2時(shí)，聲輻射力受ka的影響相對(duì)較小。

2.3.2 聲輻射力與粒子沿軸向偏離波束中心距離的關(guān)系 利用FDTD方法，對(duì)高斯聲場(chǎng)中剛性球形粒子的聲輻射力隨粒子沿軸向偏移波束中心的距離的變化關(guān)系進(jìn)行了研究。數(shù)值仿真中，取高斯脈沖波的束腰半徑w=λ。剛性球形粒子沿著傳輸軸從90網(wǎng)格移動(dòng)到610網(wǎng)格處。剛性球形粒子所受到的聲輻射力與偏移距離的關(guān)系如圖7所示。

圖7 軸向聲輻射力隨y/a的變化關(guān)系Fig.7 The relationship of axial radiation force and y/a

圖7中橫坐標(biāo)是粒子距離波束中心的軸向距離與粒子半徑之比。縱坐標(biāo)是軸向聲輻射力。從圖7的計(jì)算結(jié)果可見(jiàn)，粒子距離波束中心越近，軸向聲輻射力值越大，當(dāng)粒子在軸向逐漸遠(yuǎn)離波源時(shí)，軸向聲輻射力的值逐漸減小。粒子的聲輻射力與粒子到波源的距離成反比關(guān)系。

2.3.3 聲輻射力與剛性球形粒子橫向離軸距離關(guān)系

為了研究球形粒子位于聲場(chǎng)任意位置處聲輻射力的變化規(guī)律，對(duì)粒子分別在軸向和橫向上偏移波束中心不同的距離時(shí)軸向和橫向聲輻射力的變化進(jìn)行了仿真，計(jì)算結(jié)果如圖8和圖9所示。

圖8 不同偏移距離下軸向聲輻射力與x/a的變化關(guān)系Fig.8 The relationship of axial radiation force and different x/a

從圖8的仿真結(jié)果可見(jiàn)，當(dāng)剛性球形粒子在軸向距離波束中心的距離不變時(shí)，在波束軸上，即x/a=0，軸向的聲輻射力值最大。當(dāng)粒子從波束中心向兩側(cè)偏移時(shí)，軸向聲輻射力隨粒子橫向偏離波束中心距離的增大而逐漸減小，且在兩側(cè)相同橫向位置的軸向聲輻射力近似相等。當(dāng)橫向距離不變時(shí)，軸向聲輻射力與軸向距離成反比。仿真結(jié)果與文獻(xiàn)[13-15]的數(shù)值積分計(jì)算結(jié)果一致。圖8中的軸向聲輻射力不完全對(duì)稱是由于時(shí)域有限差分算法是將計(jì)算區(qū)域離散網(wǎng)格化，當(dāng)向兩側(cè)移動(dòng)球形粒子時(shí)，兩側(cè)所取的球上網(wǎng)格節(jié)點(diǎn)數(shù)略有不同，會(huì)產(chǎn)生兩側(cè)不完全對(duì)稱的情況，從而使計(jì)算結(jié)果稍有不同。

圖9 不同偏移距離下橫向聲輻射力與x/a的變化關(guān)系Fig.9 The relationship of transverse radiation force and different x/a

圖9是橫向聲輻射力與粒子橫向偏移距離的變化關(guān)系圖，圖中的3條曲線分別對(duì)應(yīng)的是軸向偏離波束中心1a、2a、5a的距離時(shí)，橫向聲輻射力隨橫向偏離波束中心距離的變化關(guān)系曲線。從圖9的仿真結(jié)果可見(jiàn)，剛性球形粒子在波束軸上時(shí)，三條曲線的橫向聲輻射力都近似為0。當(dāng)粒子從波束中心位置向兩側(cè)橫向移動(dòng)，球形粒子受到的橫向輻射力大小相等，方向相反。隨著偏離波束中心距離的增大，橫向聲輻射力的值逐漸增大。當(dāng)球形粒子距離波束中心的距離約為粒子半徑的3倍時(shí)，力達(dá)到最大值。隨著偏移的距離繼續(xù)增大，粒子受到的橫向聲輻射力開(kāi)始減小。從圖中可見(jiàn)，兩側(cè)相同位置的橫向聲輻射力是一對(duì)大小相等，方向相反的力，在相同橫向位置處，不同的軸向距離，橫向聲輻射力的值是不同的，粒子在軸向距離波束中心1a時(shí)橫向聲輻射力最大，在軸向距離波束中心5a時(shí)，橫向聲輻射力最小，橫向聲輻射力的大小與軸向距離波束中心的距離成反比。

3 結(jié)論

本文以聲波理論為基礎(chǔ)，應(yīng)用時(shí)域有限差分法，結(jié)合應(yīng)力張量方程，計(jì)算了離軸高斯波束作用下水中剛性球形粒子的聲輻射力。仿真研究結(jié)果表明，剛性球形粒子在高斯聲場(chǎng)中的聲輻射力隨著粒子沿軸向偏離波束中心距離的增大而減小。剛性球形粒子位于聲場(chǎng)任意位置時(shí)，軸向聲輻射力與軸向和橫向偏移距離的變化成反比關(guān)系。橫向聲輻射力隨橫向偏移距離的增加是先增大后減小。在波束中心兩側(cè)相同距離處，橫向輻射力大小相等，方向相反。時(shí)域有限差分法是一種簡(jiǎn)單、直觀、適用性廣的數(shù)值方法，能直觀模擬聲波與物體相互作用的過(guò)程，該研究將為時(shí)域有限差分法在不同形狀、不同材料的多個(gè)粒子間聲輻射力的計(jì)算提供依據(jù)。

[1] Marston L Philip. Negative axial radiation forces on solid spheres and shells in a Bessel beam[J]. Journal of the Acoustic Society America, 2007, 122:3162-3168.

[2] Mitri F G, Silva G T. Off-axial acoustic scattering of high-order Bessel vortex beam by a rigid sphere[J]. Wave Motion, 2011, 48(5):392-400.

[3] Wu Juru, Du Gonghuan. Acoustic radiation force on a small compressible sphere in a focused beam[J]. Journal of the Acoustic Society America, 1990, 87(3):997-1003.

[4] Lee J, Shung K K. Radiation forces exerted on arbitrarily located sphere by acoustic tweezers [J]. Journal of the Acoustic Society America, 2006, 120 (2): 1084-1094.

[5] Lee J, The S Y, Lee A, et al. Single beam acoustic trapping[J]. Applied Physics Letter, 2009, 95: 123-126.

[6] Yee K S. Numerical solution of initial boundary value problems involving Maxwell equations in isotropic media[J]. IEEE Transactions on Antennas and Propagation, 1996, 14(3):302-307.

[7] Cai Feiyan, Meng Long, Jiang Chunxiang, et al. Computation of the acoustic radiation force using the finite-difference time-domain method[J]. Journal of the Acoustic Society America, 2010, 128 (4): 1617-1622.

[8] Cheng Xin, Cai Feiyan, Meng Long, et al. Computation of the acoustic radiation force on a sphere based on the 3-D FDTD Method[C]//2010全國(guó)壓電和聲波理論及器件技術(shù)研討會(huì)論文集.廈門(mén):廈門(mén)大學(xué), 2010: 236-239.

[9] 黃坤朋，趙越喆. 聲學(xué)FDTD法完全匹配層數(shù)和衰減系數(shù)對(duì)衰減性能的影響[J].華南理工大學(xué)學(xué)報(bào):自然科學(xué)版，2011,39(4):135-139.

[10] 陳冬梅，張小鳳，張光斌,等.水中剛性柱在高斯聲場(chǎng)中的聲輻射力特性[J].陜西師范大學(xué)學(xué)報(bào):自然科學(xué)版，2014,42(2)：27-32.

[11] 陳冬梅，張小鳳,張光斌,等.水中球形粒子對(duì)高斯波束的聲散射[J]. 陜西師范大學(xué)學(xué)報(bào):自然科學(xué)版，2013,41(4)：31-34.

[12] 宋智廣,張小鳳,張光斌.高斯波束對(duì)水中柱形粒子的聲輻射力研究[J].壓電與聲光,2013,35(6):792-796.

[13] Wei W, Thiessen D B, Marston P L.Acoustic radiation force on a compressible cylinder in a standing wave[J]. Journal of the Acoustic Society America, 2004,116: 201-208.

[14] Mahdi Azarpeyv, Mohammad Azarpeyvand. Acoustic radiation force on a rigid cylinder in a focused Gaussian beam[J]. Journal of Sound and Vibration, 2013, 332: 2338-2349.

[15] Roumeliotis J A, Agissilaos-Georgios P Ziotopoulos, Kokkorakis G C.Acoustic scattering by a circular cylinder parallel with another of small radius[J]. Journal of the Acoustic Society America, 2001, 109(3): 870-877.

〔責(zé)任編輯 李 博〕

Computation of the off-axial acoustic radiation force on a rigid spherical particles via FDTD

SUN Xiuna, ZHANG Xiaofeng, ZHANG Guangbin*, WANG Yan, CHANG Guodong

(School of Physics and Information Technology, Shaanxi Normal University,Shaanxi Key Laboratory of Ultrasonics, Xi′an 710119, Shaanxi, China)

In order to accurately and simply calculate the acoustic radiation force,the FDTD (Finite Difference Time Domain) method is applied to discrete the acoustic wave equations in both spatial and temporal domain. By setting the appropriate boundary conditions, the pressure and velocity are computed by numerical iteration in time domain. The two-dimensional acoustic radiation force on a rigid spherical particle in water is calculated by the momentum-flux tensor equations and the finite difference domain method. The influence of frequence and shift distance from beam center on acoustic radiation force are obtained for rigid sphere particle.The results show that the axial acoustic radiation force decreases with the axial shift distance. When the particle is moved away from the beam axial, the axial acoustic radiation force is reduced. The transverse radiation force is increasing first then decreasing as the shift distance increased.The directions of the transverse force are opposite and the amplitudes of them are same when the particle shifts from the beam center with same distance.Key words: FDTD; Gaussian beam; off-axial acoustic radiation force; numerical simulation

1672-4291(2015)04-0028-06

10.15983/j.cnki.jsnu.2015.04.241

2015-01-24

中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(GK201302049); 陜西省自然科學(xué)基金(2012JM1013)

孫秀娜，女，碩士研究生，研究方向?yàn)楣β食?。E-mail:819168301@qq.com

*通信作者：張光斌，男，副教授，博士。E-mail：guangbinzhang@snnu.edu.cn

O422.2; TB556

A