非線性聲場(chǎng)下包膜微泡動(dòng)力學(xué)與頻率響應(yīng)分析*

索鼎杰 紀(jì)鎮(zhèn)祥 黃曉雲(yún) 靳杰 閆天翼?

1) (北京理工大學(xué)醫(yī)學(xué)技術(shù)學(xué)院,北京 100081)

2) (北京理工大學(xué)前沿交叉科學(xué)研究院,北京 100081)

包膜微泡在非線性聲場(chǎng)下的響應(yīng)對(duì)高強(qiáng)度聚焦超聲治療具有重要意義.本文通過(guò)耦合Gilmore-Akulichev-Zener 模型與脂質(zhì)包膜的非線性模型,采用KZK 方程仿真非線性聲場(chǎng),分析在不同超聲頻率、不同聲壓以及不同包膜材料的黏彈性下的微泡動(dòng)力學(xué)行為和微泡振蕩頻率,并進(jìn)一步對(duì)比了實(shí)際測(cè)量聲場(chǎng)與KZK 方程仿真聲場(chǎng)下的微泡動(dòng)力學(xué)行為和頻率響應(yīng).研究結(jié)果表明: 非線性聲場(chǎng)會(huì)導(dǎo)致微泡壁的瞬時(shí)運(yùn)動(dòng)速率減小,聲壓和頻率的改變對(duì)微泡動(dòng)力學(xué)的影響與在線性聲場(chǎng)中類似;包膜材料的不同可以使振蕩頻率中的諧波分量發(fā)生改變,其中包膜材料的彈性對(duì)微泡的頻率響應(yīng)影響較小,包膜的初始黏性和初始表面張力對(duì)微泡的振蕩頻率分布影響較大,當(dāng)初始黏性越小時(shí),二次分諧波的峰值越高,當(dāng)初始表面張力越大時(shí),主頻的峰值越高.本研究進(jìn)一步闡明非線性超聲激勵(lì)包膜微泡的微泡動(dòng)力學(xué),為包膜微泡在非線性聲場(chǎng)下的頻率響應(yīng)分析奠定理論基礎(chǔ).

1 引言

聲空化[1-3]是指超聲波與傳播介質(zhì)中的氣核之間的相互作用,在負(fù)聲壓的作用下使介質(zhì)中氣核迅速擴(kuò)張,出現(xiàn)空腔的現(xiàn)象;聲空化分為穩(wěn)態(tài)空化和瞬態(tài)空化.穩(wěn)態(tài)空化[4,5]是指氣泡穩(wěn)定的振蕩,存在的時(shí)間較長(zhǎng),可用于聚焦超聲誘導(dǎo)的血腦屏障開(kāi)放,通過(guò)調(diào)節(jié)穩(wěn)態(tài)空化,可以降低腦損傷的風(fēng)險(xiǎn),提高藥物的靶向遞送率;瞬態(tài)空化[6,7]是指聲壓達(dá)到慣性空化閾值后,氣泡劇烈的振蕩,存在的時(shí)間較短,瞬態(tài)空化可以產(chǎn)生沖擊波和微射流,可用于聚焦超聲組織摧毀[8,9]、血栓消融[10,11]、腫瘤治療[12,13]等.此外,空化效應(yīng)還可能導(dǎo)致光學(xué)效應(yīng)(聲致發(fā)光)[14,15]、化學(xué)效應(yīng)[16,17]和其他熱效應(yīng)[18,19].因此,對(duì)聲空化進(jìn)行研究十分重要.

近年來(lái),許多研究人員對(duì)單個(gè)微泡在超聲激勵(lì)下的動(dòng)力學(xué)模型進(jìn)行了研究和修正.Rayleigh-Plesset (RP)[20]方程通過(guò)將液體視為不可壓縮流體被廣泛運(yùn)用于計(jì)算微泡動(dòng)力學(xué)的基礎(chǔ).Keller和Miksis[21]通過(guò)拓展RP 方程,將液體視為輕微可壓縮流體,使得方程在大氣泡和大振幅下更加適用.Doinikov[22]通過(guò)拓展Keller 和Miksis[21]的研究,將氣泡的平移運(yùn)動(dòng)方程和徑向運(yùn)動(dòng)方程進(jìn)行耦合,研究了氣泡平移運(yùn)動(dòng)對(duì)微泡動(dòng)力學(xué)的影響,隨后通過(guò)研究次級(jí)Bjerknes 力[23]和聲學(xué)流線的特性,表明次級(jí)Bjerknes 力對(duì)氣泡之間的相互作用起到了重要的作用,并且聲學(xué)流線的形成與次級(jí)Bjerknes 力密切相關(guān),目前該模型已被廣泛應(yīng)用.Yang 和Church[24]通過(guò)耦合Keller-Miksis (KM)方程與Kelvin-Voigt 方程,研究了軟組織中的聲空化閾值[25].Zilonova 等[26]通過(guò)耦合Gilmore-Akulichev 方程與Zener 黏彈性模型,使得方程在高馬赫數(shù)的情況下適用,并研究了高強(qiáng)度聚焦超聲在組織中的聲空化閾值.

為了降低聲空化閾值,通常使用微米涂層脂質(zhì)包膜泡作為造影劑來(lái)提高空化治療時(shí)的安全性、穩(wěn)定性和可控性.影響脂質(zhì)包膜微泡聲空化的主要因素包括: 膜的黏性、膜的彈性、膜的初始表面張力等.Marmottant 等[27,28]提出了單個(gè)超聲造影劑氣泡的模型,該模型考慮了氣體微氣泡上脂質(zhì)包膜的物理特性,目前已被廣泛地應(yīng)用.于潔等[29]在此基礎(chǔ)上對(duì)脂質(zhì)包膜材料的黏性參數(shù)進(jìn)行了非線性修正,使得模型的應(yīng)用更加廣泛.秦對(duì)等[30]通過(guò)耦合Gilmore-Akulichev-Zener (GAZ)模型與脂類包膜的非線性黏彈性模型,研究了微泡的聲空化動(dòng)力學(xué)行為以及微泡周?chē)M織內(nèi)機(jī)械應(yīng)力.上述研究基本都是在線性波的驅(qū)動(dòng)下進(jìn)行數(shù)值分析,然而,高強(qiáng)度聚焦超聲在傳播時(shí),會(huì)出現(xiàn)衍射、散射、衰減以及聲吸收的情況,導(dǎo)致波形畸變和諧波,產(chǎn)生非線性.通常用于模擬聚焦換能器發(fā)射超聲波的兩種常用非線性模型分別是KZK (Khokhlov-Zabolotskaya-Kuznetsov)[31]方程和Westervelt[32]方程.2019 年,Bakhtiari-Nejad 和Shahab[33]通過(guò)求解KZK 方程的數(shù)值解耦合KM 方程,研究了氣泡在非線性聲場(chǎng)下的徑向運(yùn)動(dòng)和平移運(yùn)動(dòng).

本文主要采用KZK 方程進(jìn)行高強(qiáng)度聚焦聲場(chǎng)仿真,耦合GAZ 方程與包膜的黏彈性模型,采用步長(zhǎng)控制算法和Runge-Kutta5 (RK5)算法,探討在非線性聲場(chǎng)中不同超聲頻率、不同聲壓以及不同包膜材料的黏彈性下的微泡動(dòng)力學(xué)行為和微泡振蕩頻率響應(yīng),進(jìn)一步分析不同超聲參數(shù)和包膜材料對(duì)微泡振蕩頻率中諧波分量的影響,最后采用實(shí)際測(cè)量的非線性聲場(chǎng)數(shù)據(jù)與KZK 方程仿真的聲場(chǎng)數(shù)據(jù)進(jìn)行對(duì)比,探究實(shí)際聲場(chǎng)和仿真聲場(chǎng)之間的差異.

2 理論模型

2.1 GAZ 模型

在數(shù)值仿真中,由于高強(qiáng)度聚焦超聲會(huì)導(dǎo)致微泡的大幅度振蕩,氣泡壁的運(yùn)動(dòng)速度高于1 Mach(1 Mach=340.3 m/s),GA 模型更適用高馬赫環(huán)境,并且Zener[26]模型相比于Kelvin-Voigt 模型能更好地適用于黏彈性組織,能較好地模擬微泡的徑向應(yīng)力,因此本研究采用GAZ 模型:

其中,c0是流體中的聲速,R0是初始半徑,H是氣泡壁的焓變,μ是介質(zhì)中的黏性,p0是靜態(tài)壓力,σ0是初始表面張力,ρ是組織密度,B和n是常數(shù).B=-p0,p1是微泡內(nèi)部的壓力,p是驅(qū)動(dòng)聲壓,τrr|R是徑向應(yīng)力,φ是組織的松弛時(shí)間,且q=.

2.2 包膜的非線性黏彈模型

包膜材料能提高微泡的穩(wěn)定性,其黏彈特性能夠減小氣泡的壓縮和擴(kuò)張.非線性修正的Kelvin-Voigt[28,29]模型能夠較好地考慮包膜微泡的黏彈性[34,35]:

其中,S1為包膜彈性項(xiàng),S2為黏性項(xiàng),且包膜表面張力σ(R) 的徑向變化如下:

其中,Rb為微泡半徑最小值,且Rb=,當(dāng)R≤Rb時(shí),微泡發(fā)生皺縮,Rr為微泡破裂半徑,且Rr=,且σr=σ0,當(dāng)R≥Rr時(shí),微泡發(fā)生破裂.其中k表示包膜微泡的黏度,λ為包膜微泡的彈性,且k=,k0為殼的初始黏性參數(shù),α為特征時(shí)間常數(shù).

2.3 驅(qū)動(dòng)聲壓和數(shù)值條件

驅(qū)動(dòng)聲壓采用KZK 方程的模擬,通過(guò)HIFUSimulator (Food and Drug Administration,Silver Spring,MD)[36]進(jìn)行仿真:

其中p是在介質(zhì)中的聲壓,z是軸向坐標(biāo)(垂直于換能器表面),x和y為橫向坐標(biāo)(平行于換能器表面),且延遲時(shí)間τ=t-z/c0,β為非線性系數(shù),b為耗散系數(shù).通過(guò)頻域的傅里葉級(jí)數(shù)解法展開(kāi)可得

其中Cn為第n次諧波的振幅[37].

假設(shè)微泡內(nèi)部的壓力服從范德瓦耳斯方程[38],則可以采用下式表示:

其中,h=R0/8.4 是微泡的范德瓦耳斯半徑,γ是多向性系數(shù).

表1 參數(shù)名稱和值Table 1.Name and value of the parameters.

3 仿真結(jié)果分析與討論

3.1 非線性超聲的響應(yīng)

超聲波在傳播過(guò)程中會(huì)出現(xiàn)衍射、散射、衰減現(xiàn)象,從而產(chǎn)生波形畸變,采用KZK 方程模擬超聲波在水中的非線性傳播.圖1 對(duì)比了線性超聲和非線性超聲在相同正壓下的微泡徑向運(yùn)動(dòng)行為和振蕩頻率響應(yīng)(其中線性聲場(chǎng)采用驅(qū)動(dòng)聲壓為:P=-Ptsin(2πft),Pt為振幅,f=1 MHz),初始半徑R0=1 μm,聲壓P=2.0 MPa.對(duì)比線性聲場(chǎng)和非線性聲場(chǎng)下的微泡動(dòng)力學(xué)行為,可以發(fā)現(xiàn)超聲波傳播時(shí)的損耗導(dǎo)致負(fù)壓偏小,從而使得微泡的徑向運(yùn)動(dòng)減弱,氣泡壁的運(yùn)動(dòng)速度也相應(yīng)減小,微泡的振蕩頻率峰值也隨之減小.圖2 對(duì)比了非線性聲場(chǎng)下不同聲壓和不同頻率下的微泡動(dòng)力學(xué)行為和頻率響應(yīng),初始半徑R0=1 μm,當(dāng)頻率固定在f=1.0 MHz,聲壓增大時(shí),微泡的徑向運(yùn)動(dòng)顯著增強(qiáng),從而使得振蕩響應(yīng)頻率的峰值也增大.當(dāng)聲壓固定在2.0 MPa 時(shí),頻率增大時(shí),微泡的徑向運(yùn)動(dòng)卻開(kāi)始減小,這是由于微泡的振蕩主頻峰值減小和遠(yuǎn)離共振頻率導(dǎo)致的.

圖1 線性聲場(chǎng)和非線性聲場(chǎng)下微泡動(dòng)力學(xué)行為頻率響應(yīng) (a) 相同正壓下的線性超聲與非線性超聲波形圖;(b),(c) 微泡動(dòng)力學(xué)行為的對(duì)比;(d) 微泡半徑變化的振蕩頻率響應(yīng),藍(lán)色線條表示驅(qū)動(dòng)聲壓為KZK 方程仿真的非線性聲場(chǎng),紅色線條表示線性聲場(chǎng)Fig.1.Dynamics behavior and frequency response of microbubbles under linear and nonlinear ultrasound fields: (a) Linear and nonlinear ultrasonic waveforms under the same positive pressure;(b),(c) comparison of bubble dynamics behaviors;(d) oscillation frequency response of the microbubble radius change,the blue lines represent the driving sound pressure as the nonlinear ultrasound field simulated by the KZK equation,and the red lines represent the linear ultrasound field.

圖2 (a),(b),(e) 不同聲壓下非線性超聲的微泡動(dòng)力學(xué)行為和頻率響應(yīng)的差異;(c),(d),(f) 不同頻率下非線性超聲的微泡動(dòng)力學(xué)行為和頻率響應(yīng)差異Fig.2.(a),(b),(e) Differences in the kinetic behavior and frequency response of microbubbles in nonlinear ultrasound at different acoustic pressures;(c),(d),(f) differences in microbubble kinetic behavior and frequency response of microbubbles in nonlinear ultrasound at different frequencies.

3.2 脂質(zhì)包膜參數(shù)影響

圖3 對(duì)比了不同脂質(zhì)包膜微泡的初始表面張力和初始黏性參數(shù)對(duì)微泡動(dòng)力學(xué)行為和頻率響應(yīng)的影響,初始條件設(shè)置如下:f=1 MHz,P=2.0 MPa,R0=1 μm.可以發(fā)現(xiàn)當(dāng)膜的初始表面張力較大時(shí),微泡的徑向運(yùn)動(dòng)增強(qiáng),頻率響應(yīng)中諧波的數(shù)量增加;當(dāng)包膜的初始黏性減小時(shí),微泡的徑向運(yùn)動(dòng)更加劇烈,頻率響應(yīng)中二次分諧波的峰值升高.為了探究包膜材料對(duì)頻率響應(yīng)的影響,采用同樣標(biāo)準(zhǔn)來(lái)判斷,即C=(P11max1表示最大響應(yīng)頻率fmax1對(duì)于的頻譜值,P11max2表示第二響應(yīng)頻率fmax2對(duì)于的頻譜值),若|C|越大,則微泡的頻譜峰值更高,若小于C＜0 且|C|越小,則表示分諧波得到了增強(qiáng),若C＞0 且|C|越小,則表示次諧波得到了增強(qiáng).圖4 對(duì)比了不同包膜參數(shù)對(duì)微泡動(dòng)力學(xué)行為和頻率響應(yīng)的影響,仿真條件都是:f=1 MHz,P=2.0 MPa,R0=1 μm.圖4(a),(b)分別表示包膜微泡的初始表面張力σ0在0—0.07 N/m 和彈性λ在0—1 N/m 時(shí)微泡動(dòng)力學(xué)行為和頻率響應(yīng).可以發(fā)現(xiàn)包膜微泡的初始表面張力越大,微泡的徑向運(yùn)動(dòng)越劇烈,分諧波的頻率響應(yīng)越大,但包膜微泡的彈性對(duì)微泡動(dòng)力學(xué)行為和頻率響應(yīng)幾乎沒(méi)有影響.圖4(c),(d)分別表示包膜的初始黏性k0在0.5×10-8—1.0×10-7kg/s 和初始彈性λ在0—1 N/m 時(shí)微泡動(dòng)力學(xué)行為和頻率響應(yīng).可以發(fā)現(xiàn)當(dāng)脂質(zhì)包膜的初始黏性增大時(shí),微泡的徑向運(yùn)動(dòng)減小,頻率響應(yīng)中超諧波的分量增強(qiáng);當(dāng)包膜的初始黏性減小時(shí),微泡的徑向運(yùn)動(dòng)增強(qiáng),頻率響應(yīng)也提高了分諧波的分量.圖4(e),(f)分別表示包膜的初始黏性k0在0.5×10-8—1.0×10-7kg/s和包膜微泡的初始表面張力σ0在0—0.07 mN/m時(shí)的微泡動(dòng)力學(xué)行為和頻率響應(yīng),可以發(fā)現(xiàn)當(dāng)包膜的初始黏性和初始表面張力一起改變時(shí)(包膜的初始黏性小且包膜微泡的初始表面張力越大),微泡振蕩的主頻和次諧波分量的峰值變高,因此在進(jìn)行聚焦超聲消融時(shí),可以通過(guò)改變脂質(zhì)微泡的包膜材料,進(jìn)一步調(diào)節(jié)信號(hào)的頻率分量,從而使分諧波或次諧波靠近共振頻率,提供聲空化效應(yīng).

圖3 不同包膜材料對(duì)微泡動(dòng)力學(xué)行為和頻率響應(yīng)影響 (a),(b) 包膜微泡的初始表面張力;(c),(d) 包膜的黏性參數(shù)Fig.3.Effects of different coating materials on bubble dynamic behavior and frequency response of encapsulated microbubbles:(a),(b) Initial surface tension of encapsulated microbubbles;(c),(d) viscosity parameters of encapsulated microbubbles.

圖4 不同包膜材料組合對(duì)微泡動(dòng)力學(xué)行為和頻率響應(yīng)影響 (a),(c),(e) 顏色表示微泡半徑變化;(b),(d),(f) 顏色表示C 的變化Fig.4.Effects of bubble shell on bubble dynamics and frequency response in the nonlinear ultrasound field: (a),(c),(e) Color bar of is the change of radius;(b),(d),(f) color bar is the change of C.

3.3 聲場(chǎng)測(cè)量結(jié)果與KZK 方程仿真對(duì)比

聲場(chǎng)實(shí)際測(cè)量的步驟為: 首先,信號(hào)發(fā)生器的輸入信號(hào)通過(guò)功率放大器進(jìn)行放大并發(fā)送到匹配單元和聚焦換能器;然后,采用三軸定位系統(tǒng)將水聽(tīng)器對(duì)準(zhǔn)聚焦換能器的焦點(diǎn)位置,最后進(jìn)行聲場(chǎng)測(cè)量.為了防止水中的氣泡發(fā)生空化,產(chǎn)生諧波,影響數(shù)值仿真,聚焦換能器和水聽(tīng)器都在除氣水中進(jìn)行測(cè)量.圖5(a)對(duì)比了真實(shí)的聚焦換能器在焦點(diǎn)處的超聲波形和KZK 方程仿真出的焦點(diǎn)處超聲波形,仿真參數(shù)如下:f=1 MHz,P=2.0 MPa,R0=1 μm.圖5(b),(c)是實(shí)際測(cè)量的聲場(chǎng)和KZK 方程仿真聲場(chǎng)下的微泡徑向運(yùn)動(dòng)對(duì)比,可以發(fā)現(xiàn)差距不大,但仿真下的氣泡壁運(yùn)動(dòng)速度與實(shí)際相差0.5 MHz.圖5(d)是微泡動(dòng)力學(xué)行為的脈沖響應(yīng),可以發(fā)現(xiàn)真實(shí)的聲場(chǎng)在0.1 MHz 的頻譜響應(yīng)較大,而在主頻1 MHz 的響應(yīng)低于KZK 方程仿真的聲場(chǎng).這可能是由于信號(hào)發(fā)生器發(fā)射的是脈沖類型,并且電壓在起始位置會(huì)發(fā)生損耗,導(dǎo)致聲壓的波形不一致,使得整個(gè)脈沖持續(xù)時(shí)間被當(dāng)成是一個(gè)波形,導(dǎo)致0.1 MHz 的頻譜響應(yīng)較大.

圖5 實(shí)際測(cè)量值與KZK 仿真對(duì)比圖 (a) 相同正壓下實(shí)際測(cè)量與KZK 方程仿真波形圖;(b),(c) 微泡動(dòng)力學(xué)行為的對(duì)比;(d) 微泡半徑變化的振蕩頻率響應(yīng)Fig.5.Actual measured values are compared with the KZK simulation: (a) Actual measurement and KZK equation simulation waveform diagram under the same positive pressure;(b),(c) comparison of the dynamic behavior of microbubbles;(d) oscillation frequency response of the microbubble radius change.

4 總結(jié)

微泡動(dòng)力學(xué)分析是聚焦超聲應(yīng)用于聲空化、開(kāi)放血腦屏障、消融以及聲致發(fā)光中的重要理論基礎(chǔ).考慮在高強(qiáng)度聚焦超聲治療中使用聚焦超聲換能器發(fā)射的超聲由于能量較高會(huì)在介質(zhì)中傳播時(shí)產(chǎn)生非線性畸變,本文主要采用KZK 方程進(jìn)行高強(qiáng)度聚焦聲場(chǎng)仿真,耦合Gilmore-Akulichev-Zener與包膜的黏彈性模型.相比之前的研究,我們?cè)诳紤]超聲波傳播非線性效應(yīng)的同時(shí)加入了脂質(zhì)微泡模型,可以進(jìn)行非線性聲場(chǎng)的脂質(zhì)微泡的微泡動(dòng)力學(xué)仿真和微泡振蕩頻率響應(yīng),使微泡動(dòng)力學(xué)仿真和頻率響應(yīng)下的環(huán)境更加趨近于真實(shí)聲場(chǎng).研究結(jié)果表明: 當(dāng)包膜材料不變時(shí),超聲波傳播時(shí)的非線性會(huì)導(dǎo)致微泡的擴(kuò)張減弱,但會(huì)產(chǎn)生更多的諧波分量,包膜微泡的振蕩頻率分量中諧波含量增大.當(dāng)超聲聲壓增大時(shí),微泡的徑向運(yùn)動(dòng)增強(qiáng),各頻率的分量也增高;當(dāng)超聲頻率靠近振蕩頻率時(shí),微泡的徑向運(yùn)動(dòng)增強(qiáng),頻率分量中主頻的峰值增高.

當(dāng)改變包膜材料時(shí),包膜微泡的初始黏性越小,振蕩頻率中分諧波的分量越大,包膜微泡的初始表面張力越大,振蕩頻率中出現(xiàn)的諧波分量越多,分諧波的分量越多.通過(guò)調(diào)節(jié)不同的包膜參數(shù)材料組合,可以發(fā)現(xiàn)包膜微泡的徑向運(yùn)動(dòng)行為與振蕩頻率的改變比較一致,當(dāng)諧波靠近共振頻率時(shí),微泡的聲空化效應(yīng)得以增強(qiáng);振蕩主頻分量越高,微泡的聲空化效應(yīng)也會(huì)相應(yīng)增強(qiáng).本研究進(jìn)一步為非線性超聲消融奠定了理論基礎(chǔ).由于本文以單個(gè)微泡作為對(duì)象,且未將微泡內(nèi)部與外部的熱傳遞和熱效應(yīng)考慮進(jìn)去,后續(xù)將重點(diǎn)考慮多微泡的徑向運(yùn)動(dòng)和平移運(yùn)動(dòng),并結(jié)合微泡內(nèi)部與外部的熱效應(yīng),進(jìn)一步分析多微泡在非線性聲場(chǎng)激勵(lì)下的微泡動(dòng)力學(xué)行為和熱能傳遞.