張 潔，劉春澤，王 歡，許小芳，周紅生

（1. 中國科學院聲學研究所東海研究站，上海 201815；2. 中國科學院大學，北京 100049）

0 引 言

超聲操控技術是對微小物體進行非接觸操控的一種新型重要手段[1]。在生物醫(yī)學領域，該技術可用于靶向給藥的治療方式[2]，將藥物載體通過超聲場操控移動到指定病灶部位后釋放，實現(xiàn)高效、無創(chuàng)、精準治療。其原理是利用高強、高頻聲波的非線性效應產(chǎn)生聲輻射力來實現(xiàn)對微粒的捕獲和操控[3]。當一束聲波入射到物體表面時，物體會對聲波產(chǎn)生散射、折射和吸收作用，與聲波發(fā)生能量和動量的交換，導致物體受到一個不為0的平均壓力，稱為聲輻射力[4-5]。相比較光學、電場和磁場操控技術，超聲操控具有明顯的優(yōu)勢，能夠以更具生物相容性的方式操控更多的材料[6-7]，擁有廣闊的應用前景。

Gor’kov[8]將物體表面的聲輻射力表示為物體表面動量通量的時間平均，推導出微小球形物體在任意聲場中勢能的解析解，依據(jù)聲場分布便可確定該物體的聲勢能，而聲輻射力大小與聲勢能梯度有關，方向指向聲勢能極小值處，故而可以構造特定聲場將物體限制在指定的空間位置。駐波產(chǎn)生的聲輻射壓力遠大于行波，所以常用駐波場完成顆粒束縛。Park等[9]、Hong等[10]和Xic等[11]基于該理論設計了多種利用駐波聲場的操控裝置，分別實現(xiàn)了在空氣中對固體顆粒、液滴以及昆蟲（活體生物）的懸浮。Wu[12]和 Hcrtz[13]推導了球形顆粒在駐波聲場受到的聲輻射力的解析式，并且在液體環(huán)境下通過兩個聚焦超聲換能器完成了對微米級物體的捕獲實驗。

但是以上理論及模型均忽略了傳播介質的熱黏性作用，并且實際上對顆粒的操控不僅僅依靠聲輻射力，也受到由聲流場產(chǎn)生的曳力的影響。這里的聲流場也是聲波的非線性效應產(chǎn)生的，它是由介質的黏性衰減導致聲場的空間梯度發(fā)生變化，介質振元的凈位移不為零引起的全局流動[14-15]。文獻[16-17]在現(xiàn)有理論的基礎上利用無粘聲波散射理論的遠場解表示粒子上的輻射力，將此解擴展到靠近粒子的近場區(qū)域，然后將其與粒子聲邊界層中不可壓縮的黏性流動問題的解進行匹配，重新推導了黏性流體中的聲輻射力公式。

本文基于 Bruus的粘性介質中的聲輻射力理論，建立了超聲駐波場驅動下的水下顆粒操控模型，利用物理場仿真軟件對模型中的聲場、聲流場進行仿真，分析顆粒小球受到的聲輻射力及聲流場引起的Stokc曳力大小，以及兩種作用力的相對關系，動態(tài)模擬顆粒小球操控過程。最后設計對不同尺寸顆粒小球的水下駐波場操控實驗，驗證了顆粒操控的可行性以及聲流場對顆粒操控會產(chǎn)生一定的影響。

1 黏性介質聲輻射力理論

2 計算機仿真

本文采用基于有限元方法的多物理場仿真軟件COMSOL Multiphysics進行聲場、聲流場和顆粒操控動態(tài)過程的仿真，模型的具體仿真參數(shù)如表1所示[20]。

表1 計算機仿真的模型參數(shù)Table 1 Model parameters of computer simulation

水下顆粒操控系統(tǒng)采用三維建模，模型示意圖如圖1所示，圖中灰色區(qū)域是高為20 mm，半徑為10 mm的圓柱體水域，紫色為兩個凹球面超聲聚焦換能器，孔徑半徑為 4.53 mm，曲率半徑（焦距）為10 mm。除換能器表面外其他壁面貼附吸聲材料，在仿真中將邊界設置完美匹配層，表示無聲波反射。設定發(fā)射超聲的頻率為1 MHz，紅色箭頭為換能器聲場輻射方向。以下方換能器孔徑平面中心為原點O，通過原點O在孔徑平面內任選兩條垂直軸線定為x軸和y軸，上下?lián)Q能器中心點的連接線為z軸，建立三維直角坐標系。

圖1 水下顆粒操控系統(tǒng)示意圖（單位：mm）Fig.1 Schematic diagram of underwater particle control system (unit: mm)

2.1 聲場仿真

聲場仿真采用熱黏性聲學模塊進行建模，兩個換能器孔徑平面的中心間距為20 mm，將換能器的曲面設置為速度振動條件，通過調整速度幅值可以改變聲場聲壓值的大小，振動方向與曲面切線垂直，向空間輻射連續(xù)正弦聲波。計算得到穩(wěn)態(tài)聲場的聲壓如圖2所示。由于聲波的干涉效應，在三維聲場中z軸線上聚焦并形成局部駐波聲場，聲場內最大聲壓值0.295 MPa，其中圖2（a）為三維聲場聲壓圖，圖2（b）為xOz切平面的聲壓圖。

圖2 聚焦聲場聲壓圖Fig.2 Focused acoustic pressure diagram

2.2 聲流場仿真

將上文計算得到的穩(wěn)態(tài)聲場作為源場，利用層流模塊對Navicr-Stokcs方程的時均二階擾動方程進行求解，在模型中添加體積力和質量源，最終得到時均二階速度場v2即所謂的聲流速度場。為減少計算機運算量且該三維模型具有軸對稱性，故采用xOz切平面聲場解進行流場仿真。設定初始條件為域內流體速度為0，根據(jù)式（3）～（4）計算得到的穩(wěn)態(tài)流場如圖3（a）所示，圖中背景為流速大小，白線為流線，與線上箭頭共同表示流動方向，靜態(tài)水域在聲輻射力的作用下從兩個換能器輻射表面開始加速，并在中心位置匯聚減速，隨后向左右兩側進行分流，碰到兩側的壁后形成回流。圖3（b）為z軸上的流速值，以軸線中點為中心鏡像對稱分布，軸線上最大流速為0.09 mm·s-1。

圖3 穩(wěn)態(tài)聲流場仿真結果Fig.3 Simulation results of steady acoustic flow field

2.3 顆粒操控仿真

操控的顆粒小球采用密度和水接近的聚苯乙烯材料，仿真中忽略顆粒重力和浮力影響，同時忽略顆粒之間的相互作用，使其僅受到聲輻射力和聲流曳力的作用，在聲場和流場穩(wěn)態(tài)解的基礎上，采用粒子追蹤模塊與瞬態(tài)研究，根據(jù)式（5）～（6）添加聲輻射力與聲流曳力，初始時刻將小球均勻分布在xOz切平面聲場內，計算6 s中的小球運動變化軌跡。

圖4為顆粒聚集動態(tài)過程圖，選取了以坐標（-1.5 mm， 8.6 mm）、（-1.5 mm， 11.4 mm）、（1.5 mm，8.6 mm）和（1.5 mm， 11.4 mm）為頂點的局部駐波場區(qū)域進行展示，背景為式（8）計算的該區(qū)域內時間平均勢能Urad。0 時刻顆粒小球均勻分布在平面內，在合力的作用下發(fā)生移動，顆粒集群運動趨勢為快速的縱向聚集，然后以相對緩慢的速度完成橫向聚集，最終在6 s時達到穩(wěn)定狀態(tài)，聚集在Urad極小值附近。移動過程中，初始位置距離勢能最低點最遠，即處于駐波場波峰或波谷位置的顆粒，獲得最大速度可達到 0.028 m·s-1，但該位置并不是獲得最大初始加速度的位置，下文予以說明。

圖4 顆粒聚集動態(tài)過程Fig.4 Dynamic process of particles aggregation

3 仿真結果分析

在顆粒聚集動態(tài)圖中，取圖4中第6 s圖中的結果作為顆粒操控的穩(wěn)態(tài)結果，選取z軸數(shù)據(jù)，畫出該軸線上的聲壓、輻射力和時間平均勢的曲線，并進行歸一化處理，結果如圖5所示。圖5中紅點標記處為顆粒聚集點，該點位于駐波場聲波波節(jié)位置。從受力角度分析，聲輻射力曲線的正負代表力的方向，該聚焦點正處于聲輻射力的正值變化為負值的臨界點處；從勢能角度分析，該點對應于平均時間勢能曲線的最小值，與時間平均勢理論相符合。顆粒在聲輻射力極值點處獲得最大初始加速度，該位置對應于駐波場中波節(jié)與波峰（或波谷）的中心處，如圖5中藍點位置。

圖5 z軸線聲壓、輻射力和時間平均勢的歸一化曲線Fig.5 Normalized curves of sound pressure, radiation force and time average potential alongz-axis

3.1 作用力分析

在xOz平面內的駐波聲場中，聲輻射力的方向并不完全與聲傳播方向（z軸）平行，可以分解到z與x兩個方向，分別稱為縱向聲輻射力與橫向聲輻射力。圖6選取了z軸上1/2波長寬度的聲場繪制了聲壓和聲輻射力圖。圖6中綠色箭頭代表聲輻射力，箭頭的方向代表箭頭中心位置處聲輻射力的方向，箭頭的長度與聲輻射力的大小成正比。

圖6 1/2波長區(qū)域內的聲壓和聲輻射力圖Fig.6 Sound pressure and radiation force map in the 1/2 wavelength region

圖7為縱、橫聲輻射力圖，由于縱、橫聲輻射力數(shù)值分布差別較大，所以對箭頭長度進行歸一化處理，綠色箭頭僅表示力的方向，圖7（a）、7（b）兩圖解釋了顆粒小球在縱向聲輻射力的推動下向著波節(jié)位置移動，并在橫向聲輻射力作用下向著z軸進行靠攏，完成顆粒聚集過程。

圖7 1/2波長區(qū)域內的聲壓和聲輻射力圖（分解后）Fig.7 Sound pressure and radiation force map in the 1/2 wavelength region after decomposition

接下來對作用力進行具體的數(shù)值分析，添加了顆粒靜止狀態(tài)（u=0）時的聲流曳力，同樣將其分解到z與x方向，稱為縱向曳力與橫向曳力。選取z軸和同x軸平行且z=10 mm的軸線上繪制一維作用力圖，力場分析結果如圖8所示。藍線為聲輻射力曲線，綠線為聲流曳力曲線，圖8（a）為縱向作用力曲線，最大聲輻射力為 3.6×10-8N，最大曳力為1.07×10-10N，右下角為曳力放大后的波形圖，前者約為后者數(shù)值的300倍，因此在縱向上聲輻射力作為顆?？刂频闹鲗ё饔昧?，將顆粒束縛在駐波的每一個波節(jié)處。圖8（b）為橫向的作用力曲線，最大聲輻射力為7.51×10-10N，約為縱向聲輻射力的1/50，最大曳力為3.9×10-11N，在該方向上兩個力的大小相對差值較小，顆粒的操控力由兩個力的合力決定?？傊?，粒子在駐波場內受到的縱向作用力要遠大于橫向作用力，故而圖4中顆粒會在極短時間內沿著縱向聚集，然后緩慢在橫向進行移動，最終匯聚為一團。

圖8 縱、橫方向上的作用力曲線圖Fig.8 Action force curves in the vertical and horizontal directions

3.2 顆粒小球半徑對操控的影響

顆粒小球在水中的操控過程，可以認為是其受到的聲輻射力與曳力的“博弈”。根據(jù)式（5）、（6）可知，顆粒小球受到的聲輻射力和聲流曳力分別與顆粒半徑rP3和rP成正比。利用超聲駐波場束縛顆粒通常要求物體尺寸小于1/3～1/4波長[21]，為研究顆粒小球半徑對操控的影響，設置半徑為 74、37、18和10 μm，計算橫縱軸線上的作用力，分析兩種力的關系，結果如圖9所示。

圖9 不同顆粒半徑下的橫縱向作用力的曲線圖Fig.9 Curves of transverse and longitudinal forces under different particle radii

圖9中藍色虛線表示聲輻射力，綠色虛線表示聲流曳力，紅色實線代表顆粒小球受到的合力，其中作用力值的正負代表力的方向，合力值為0的點表示該位置處聲輻射力與曳力平衡，即為顆粒小球的束縛點。圖9（a）中綠線幅值較藍線小很多，綠紅兩線基本重合，表明合力與聲輻射力相近，此時顆粒小球主要靠聲輻射力進行操控。圖9（b）、9（c）表示隨著顆粒尺寸的減小，聲輻射力與曳力的差值也在逐漸減小，合力發(fā)生了一些變化，此時合力的零值點有三處，一處位于中心，兩處分別偏離中心約1.5mm位置，其他區(qū)域內顆粒隨著聲流被推出駐波聲場。最后當顆粒半徑減小到10 μm時，如圖9（d）所示，僅有中心一處合力的零值點，此時顆粒放置在該平衡位置外的其他區(qū)域內，均由聲流曳力控制其移動。圖9（c）、圖9（f）分別為縱軸上半徑 74和10 μm顆粒受到的作用力，相比較橫向合力的變化，半徑減小到10 μm時仍以聲輻射力為主要作用力，但是此時由于橫向作用力以曳力為主，顆粒小球會隨著聲流運動，無法形成團聚。

3.3 聲場強度對操控的影響

有研究表明隨著聲場強度持續(xù)增加，聲流場從層流過渡到湍流時，流場分布會發(fā)生變化且更加復雜[22]，所以本文研究的聲流場僅處于層流的范圍內。

改變陣面的振速幅值可以調整聲場強度，令聲場最大聲壓值分別為0.074、0.147和0.295 MPa，大小比例為1:2:4，計算得到三種不同強度的聲場對應的聲輻射力比值為 1:4:16，聲流曳力比值同為1:4:16。在式（5）、（6）忽略式中系數(shù)，主要變量為，由于振速與聲壓存在關系：，p1和v1是聲壓和振速一階量，v2為振速的二階量，故式（5）中聲輻射力大小與p2成正比，式（6）在顆粒靜止（u=0）時，聲流曳力也與p2成正比關系，表明在聲場強度增大的過程中，聲輻射力與聲流曳力的相對大小是一致的，合力在空間內的方向分布沒有發(fā)生變化，顆粒在聲場中的束縛狀態(tài)不變。這里提到的狀態(tài)不變是因為忽略了顆粒的重力與浮力，在實際應用中需要通過調節(jié)聲場強度的大小來補償其他作用力，調節(jié)操控位置。

4 水下顆粒操控實驗

4.1 實驗平臺

本文搭建的水下顆粒操控實驗平臺如圖10所示，包括以下設備：（1） 機械夾持裝置，可以調節(jié)夾持換能器的距離及角度；（2） 凹球面超聲聚焦換能器，中心頻率為2 MHz，曲率半徑（焦距）為20 mm，孔徑半徑為9.07 mm，共軸平行相向固定；（3） 多功能信號發(fā)生器（NF-WF1974）；（4） 功率放大器（NFHSA4101）；（5）數(shù)字示波器（RIGOL-MOS5102）；（6） 有機玻璃材質透明水槽。

圖10 實驗器件放置實物圖Fig.10 Picture of experimental device placement

4.2 實驗過程及結果

實驗采用聚苯乙烯顆粒作為操控目標，選取了半徑為150、75、37.5和7.5 μm的四種顆粒。由于聚苯乙烯顆粒密度略大于水，將顆粒倒入小燒杯中加水，靜置會沉入杯底，實驗時需要對顆粒進行攪拌，使其處于懸停狀態(tài)后吸入滴管。

在實驗前向水槽中加入蒸餾水浸沒裝置，調整兩個換能器孔徑平面距離為1.5倍焦距（30 mm），連接好設備后，控制信號發(fā)生器發(fā)射頻率為2 MHz的正弦波連續(xù)信號，經(jīng)過功率放大器后，示波器顯示輸入電壓的峰-峰值為13.6 V。換能器開啟后，使用滴管吸取包含顆粒的水溶液，重復多次向聲波輻射區(qū)域添加顆粒，待顆粒的束縛狀態(tài)穩(wěn)定后進行實驗記錄。

不同尺寸顆粒的水下聚集狀態(tài)如圖11所示，半徑為150、75 μm顆粒在換能器中軸線區(qū)域被束縛，對應于聲場的局部駐波場區(qū)域，部分顆粒處于中軸上，其余顆粒稍微偏離中軸1～2 mm，位置近似于圖9（b）的三處合力零值點。隨著顆粒尺寸的減小，顆粒的聚集程度減弱，分散地分布在輻射聲場范圍內，當顆粒半徑減小到7.5 μm時，顆粒已經(jīng)無法完成束縛，添加的顆粒會隨著聲流向兩側緩慢移動。

圖11 不同尺寸顆粒的水下聚集狀態(tài)Fig.11 Underwater aggregation of particles with different sizes

同時顆粒操控的穩(wěn)定程度也隨半徑減小而減弱。其中對尺寸為 150 μm 的顆粒操控最為穩(wěn)定，可以跟隨夾持設備一同緩慢移動，但顆粒尺寸減小到37.5 μm后，受到水域的微擾動便會有粒子離開輻射區(qū)域，顆粒尺寸越小穩(wěn)定性越差。

5 結 論

本文通過計算機仿真水下顆粒操控系統(tǒng)的聲場、聲流場和顆粒操控動態(tài)過程，驗證了熱黏性流體中的聲輻射力理論和時間平均勢能理論，進行了顆粒的作用力分析與水下顆粒操控實驗，得到的結論如下：

（1） 仿真僅考慮了顆粒在水中操控過程中受到的聲輻射力與聲流曳力。顆粒在駐波聲場中會聚集并束縛在聲波波節(jié)位置，該位置對應于平均勢能場的最小值。在顆粒的縱向聚集過程中，局部駐波場的縱向聲輻射力為主導力，其數(shù)值遠大于同向的聲流曳力；橫向聚集取決于橫向聲輻射力與曳力的相對大小，在顆粒尺度變小的過程中，聲流曳力漸漸成為主導力，此時依靠聲輻射力將無法對顆粒進行操控。

（2） 仿真中流場分布沒有發(fā)生變化的前提下，聲場強度的改變，聲場內任意位置處的兩種力的相對大小不變，即合力的方向不變，顆粒的操控模式不變，合力大小與聲場聲壓值的平方成正比。但實際情況中，考慮顆粒受到的重力、浮力及環(huán)境擾動后，聲場強度的增大可以增強顆粒的操控與抗擾動能力。

（3） 水下顆粒操控實驗的結果驗證了第一條結論，換能器輻射形成的局部駐波場可以依靠聲輻射力穩(wěn)定地聚集顆粒，但隨著顆粒的尺寸減小，顆粒分布逐漸分散，聲輻射力不再作為主導力，若想對該尺寸下的顆粒進行聚集，需要依靠聲學微渦流等方式實現(xiàn)目標。