叢聚的含氣泡水對(duì)線性聲傳播的影響?

范雨喆 陳寶偉 李海森 徐超

1)(哈爾濱工程大學(xué),水聲技術(shù)重點(diǎn)實(shí)驗(yàn)室,哈爾濱 150001)2)(海洋信息獲取與安全工信部重點(diǎn)實(shí)驗(yàn)室(哈爾濱工程大學(xué)),工業(yè)和信息化部,哈爾濱 150001)3)(哈爾濱工程大學(xué)水聲工程學(xué)院,哈爾濱 150001)(2018年4月18日收到;2018年5月18日收到修改稿)

1 引 言

含氣泡水是重要的自然及工業(yè)現(xiàn)象,常見于艦船螺旋槳空化、水下爆炸和生物醫(yī)學(xué)應(yīng)用等領(lǐng)域中[1?3].關(guān)于含氣泡水中的線性聲傳播問題已有大量研究,主要基于兩種方法.其一是通過集平均[4]或體積平均[5]建立起描述含氣泡水內(nèi)的聲傳播的基本方程,在結(jié)合氣泡動(dòng)力學(xué)方程[6]后,對(duì)整體做線性化近似.當(dāng)然,這些線性聲傳播理論由于忽略了氣泡間的相互作用,會(huì)在含氣泡水濃度較高時(shí)產(chǎn)生較大誤差.通過考慮氣泡與含氣泡水中的平均量相互作用,可以初步解決這一問題[7?10].另一種方法則是先對(duì)氣泡動(dòng)力學(xué)方程做線性化近似,在此基礎(chǔ)上將含氣泡水中的聲傳播問題看作多體多次散射問題[11,12].同樣,在含氣泡水氣泡體積分?jǐn)?shù)較高時(shí),需引入氣泡間的再輻射對(duì)理論進(jìn)行修正[13?16],這種方法也被稱為多泡散射法(multiple-scattering approach,MST)[17].在不考慮氣泡間相互作用時(shí),上述兩種方法在線性聲傳播框架內(nèi)的結(jié)論是一致的.此外,也可以尋求該問題的數(shù)值求解辦法,但是其過高的計(jì)算量阻礙了精確仿真環(huán)境的建立[18].

然而,這些研究的一個(gè)共同的重要假設(shè)是:氣泡在空間均勻分布且相互間統(tǒng)計(jì)獨(dú)立.這導(dǎo)致含氣泡水中聲傳播理論大多忽略了含氣泡水內(nèi)氣泡的空間分布對(duì)聲傳播的影響.事實(shí)上,含氣泡水中氣泡在空間上分布并均勻是十分嚴(yán)苛的條件,即使在實(shí)驗(yàn)室環(huán)境中也很難達(dá)到.現(xiàn)有的含氣泡水線性聲傳播實(shí)驗(yàn)中,相速度和衰減系數(shù)的峰值往往低于理論預(yù)測(cè)值,且存在低頻偏移[19,20].通過實(shí)驗(yàn)觀測(cè)[21],多分散的含氣泡介質(zhì)內(nèi)氣泡間可能存在小范圍的聚集趨勢(shì)(在后文統(tǒng)一稱作叢聚現(xiàn)象),這為合理解釋預(yù)測(cè)和測(cè)量結(jié)論間的偏差提供了研究方向.本文基于MST框架,將準(zhǔn)晶體近似(quasi-crystalline assumption,QCA)引入自洽方法(self-consistent approach,SCA)[21,22],推導(dǎo)出了考慮空間分布時(shí)多分散介質(zhì)中的等效聲波波數(shù)形式,并在此基礎(chǔ)上引入Neyman-Scott點(diǎn)過程(Neyman-Scott point process,NSP)描述了含氣泡水內(nèi)氣泡小范圍的叢聚現(xiàn)象[23],給出了叢聚氣泡水內(nèi)等效聲波波數(shù).考慮到空間分布的統(tǒng)計(jì)信息提取對(duì)相關(guān)研究的精確與否起到重要作用,本文引入了一種比例無(wú)偏估計(jì),通過該方法獲得了仿真環(huán)境下叢聚含氣泡水模型的相速度及衰減系數(shù),該建模及統(tǒng)計(jì)方法也可作為相關(guān)實(shí)驗(yàn)工作的理論基礎(chǔ).

2 氣泡水中聲傳播的自洽方法

2.1 氣泡的小幅度振動(dòng)

考慮氣泡在聲場(chǎng)作用下的動(dòng)力學(xué)方程,假設(shè)其在振動(dòng)過程中始終保持球形,并忽略氣泡與水間的質(zhì)量交換,氣泡的徑向振動(dòng)可由Prosperetti修正的Keller-Miksis方程[6]表示為

其中pin為

其中p0為環(huán)境大氣壓,pv是泡內(nèi)蒸汽壓,σ是液體表面張力,μ是液體的黏滯系數(shù),c0為聲波在液體中的傳播速度,ρl是液體密度,γ是泡內(nèi)氣體的多方指數(shù),a0是氣泡初始半徑,a=a(t)為氣泡任意時(shí)刻的半徑,p∞為聲場(chǎng)和環(huán)境壓力.(1)式在低Mach數(shù)下,能準(zhǔn)確地模擬氣泡在聲場(chǎng)中振動(dòng)的動(dòng)力學(xué)過程.對(duì)(1)式做一階微擾處理,此時(shí)穩(wěn)態(tài)的氣泡隨聲場(chǎng)做同頻率周期振動(dòng),假設(shè)a(t)=a0(1+X),p∞=p0(1?εe?jωt)和pin=pin,equ(1?ΦX),可得氣泡的線性動(dòng)力學(xué)方程[24]

此處

其中pin,equ=p0+2σ/a0,在多方過程假設(shè)下,Re{Φ}=3γ.為建立氣泡振動(dòng)與MST理論的聯(lián)系,必須求解單個(gè)氣泡的散射聲場(chǎng).在以ω為角頻率的聲場(chǎng)中,對(duì)ri處的氣泡,將其散射聲場(chǎng)寫作ps(r)=fspinc(ri)G(r?ri)exp(?iωt),其中G(r?ri)=exp(ik|r?ri|)/|r?ri|是三維空間的格林函數(shù)[25],pinc(ri)為ri位置處的入射聲場(chǎng),ps(r)為點(diǎn)氣泡在r處氣泡的散射聲場(chǎng),fs為氣泡的聲散射幅度.對(duì)脈動(dòng)氣泡,根據(jù)Eular方程,有

表示了氣泡散射聲場(chǎng)和氣泡徑向振動(dòng)的關(guān)系.考慮到氣泡的徑向運(yùn)動(dòng)可以表示為a(t)=a0[1+X(t)],且其隨聲場(chǎng)做穩(wěn)態(tài)的同頻率周期振動(dòng),于是有ω2X=?¨X,將其代入(3)式并結(jié)合(4)式[24],又考慮到在點(diǎn)源假設(shè)下ka0?1,可以求解單個(gè)氣泡的聲散射幅度f(wàn)s為

2.2 自洽方法

圖1 含氣泡水中的聲傳播Fig.1.Acoustic propagation in bubbly water.

當(dāng)聲波在含氣泡水中傳播時(shí),每個(gè)氣泡都會(huì)對(duì)聲波進(jìn)行再輻射,而其輻射的聲波又會(huì)被其他氣泡吸收和再次輻射.如圖1所示,采用SCA時(shí),忽略時(shí)間因子exp(?iωt),將氣泡的聲散射幅度重新寫作fs(ai),其中ai為氣泡i的初始半徑,其相應(yīng)的散射聲場(chǎng)也重新寫為pis(r,ai) =fs(ai)pinc(ri)G(r?ri).考慮到等效波數(shù)并非入射角度的函數(shù),為簡(jiǎn)化問題,假設(shè)含氣泡水內(nèi)聲波沿z軸傳播,可將其在ri處的聲壓簡(jiǎn)化為

其中κ為含氣泡水內(nèi)等效波數(shù),P為相應(yīng)的聲壓幅值.由于研究對(duì)象為多分散介質(zhì),氣泡半徑ai也同樣作為隨機(jī)變量引入,各個(gè)散射體的散射幅度f(wàn)s(ai)都變?yōu)殡S機(jī)變量的函數(shù).

任取一個(gè)氣泡j作為參考點(diǎn),其空間位置變量為rj,氣泡半徑為aj,二者都為隨機(jī)變量(圖1).參考點(diǎn)與ri處大小為ai的氣泡i的相對(duì)距離為rij=|rj?ri|.在參考點(diǎn)處的聲壓可以獲得如下方程

式中pi,in為氣泡i受到的激勵(lì)聲壓,pex(rj)=εp0ejkzj為背景聲場(chǎng).對(duì)含有N個(gè)氣泡的氣泡水而言,pi,in為2N個(gè)隨機(jī)變量ri和ai的函數(shù).對(duì)參考點(diǎn)氣泡,在(7)式兩邊同時(shí)求取條件期望,可得

其中〈P〉j表示氣泡j的條件期望,是在參考?xì)馀菘臻g位置和氣泡大小一定時(shí),其余N?1個(gè)氣泡的空間位置和大小的條件概率密度,其可以進(jìn)一步寫作氣泡j一定時(shí)氣泡k的條件概率密度與氣泡j,k一定時(shí)其余N?2個(gè)氣泡的條件概率密度的乘積f(?N?1|rj,aj)=f(?N?2|rk,ak,rj,aj)f(rk,ak|rj,aj). 將其代入(8)式,并假設(shè)含氣泡水在空間上是統(tǒng)計(jì)均質(zhì)的(簡(jiǎn)稱均質(zhì)),各個(gè)氣泡的統(tǒng)計(jì)特征無(wú)明顯差異,在此基礎(chǔ)上(8)式化簡(jiǎn)為

其中〈pk,in〉jk表示對(duì)氣泡j,k取條件期望.對(duì)(9)式進(jìn)行求解,首先給出隨機(jī)變量的條件概率密度函數(shù).相對(duì)參考點(diǎn)rj,大小為ak的氣泡k位于drk的概率為[22,23]

其中g(shù)([rk;ak],[rj;aj])表示氣泡k和j在空間上的相互影響,對(duì)均質(zhì)氣泡群時(shí),g([rk;ak],[rj;aj])=g(|rk?rj|),空間影響退化為徑向函數(shù)且與具體位置無(wú)關(guān),稱作徑向關(guān)聯(lián)函數(shù)(pair-correlation function,PCF),(10)式中氣泡間的相互影響僅通過PCF表示;m(rk)為氣泡群的氣泡數(shù)量密度函數(shù),在均質(zhì)時(shí)退化為常數(shù)m0;n(ak)是氣泡孔徑的概率密度函數(shù)(這里將m0和n(ak)這兩個(gè)函數(shù)區(qū)分開表示,更方便表述含氣泡水的空間特性).因此,(9)式簡(jiǎn)化為

此時(shí),在(11)式中引入QCA[26],〈P〉jk=〈P〉j,再進(jìn)一步考慮氣泡與含氣泡水中的平均量相互作用[7?10],將〈pj,in〉j=Pexp(iκzj)代入(11)式,并去掉角標(biāo)k以做化簡(jiǎn),(11)式變換為

2.3 叢聚時(shí)含氣泡水對(duì)聲傳播的影響

對(duì)叢聚的含氣泡水,引入NSP對(duì)這一現(xiàn)象進(jìn)行描述[23].采用這一方法,避免了由于氣泡間的復(fù)雜受力[28]導(dǎo)致的含氣泡水空間結(jié)構(gòu)研究上的困難[29],直接給出了不同情況下含氣泡水的統(tǒng)計(jì)特征,用統(tǒng)計(jì)特征建立起了含氣泡水中聲傳播與空間分布的聯(lián)系,在簡(jiǎn)化問題的基礎(chǔ)上給出了更具一般性的結(jié)論.氣泡數(shù)量密度為m的含氣泡水中分布著iNSP個(gè)小范圍的叢聚子氣泡群,我們假設(shè)這些叢聚的成因無(wú)明顯差異并在短時(shí)間內(nèi)無(wú)明顯變化,統(tǒng)計(jì)上認(rèn)為這些叢聚子氣泡群內(nèi)空間分布相同,并在空間上整體服從Poisson分布,各自的數(shù)量密度為m0.對(duì)任意的從聚子氣泡群而言,其內(nèi)部仍服從Poisson分布,在空間上的一次實(shí)現(xiàn)中含有Nc(iNSP)個(gè)氣泡.將NSP的一次實(shí)現(xiàn)寫作

其中叢聚子氣泡群nc=〈Nc〉,m0=〈iNSP〉, 且m=m0nc.本文直接給出其空間分布特征如下(數(shù)學(xué)證明可參考文獻(xiàn)[23]),考慮到引入了QCA,因此只需要考慮二階統(tǒng)計(jì)量,PCF為

這里

其中,R為氣泡聚集范圍的半徑.將(15)式代入(13)式中,就得到了叢聚的含氣泡水內(nèi)的等效波數(shù)

3 數(shù)值模擬

3.1 相速度和聲衰減的變化分析

在(16)式中,含氣泡水中氣泡空間分布的相關(guān)參數(shù)被引入到了等效波數(shù)中,在本文的數(shù)值模擬中,我們始終保持含氣泡水的體積分?jǐn)?shù)為5×10?4.選取邊長(zhǎng)為2.6×10?3m的觀測(cè)窗,在此基礎(chǔ)上改變nc,觀察氣泡空間分布對(duì)含氣泡水的影響,其一次實(shí)現(xiàn)的x-y面?zhèn)纫晥D如圖2所示,其中,含氣泡水空間分布的仿真方法在3.2節(jié)中說明.含氣泡水的氣泡的尺寸分布服從log-normal分布,其孔徑的概率密度函數(shù)為

其中設(shè)氣泡的孔徑分布函數(shù)中aμ=8μm,σa=0.4,叢聚半徑為20aμ.此時(shí),注意到含氣泡水小范圍內(nèi)的孔隙率已遠(yuǎn)大于5×10?4,那么氣泡間的相互作用會(huì)對(duì)聲傳播產(chǎn)生影響,為解決這一因素的影響,進(jìn)一步通過修正(5)式,考慮與背景聲場(chǎng)的相互作用下的單個(gè)氣泡的聲散射能力[15]

由圖2可以看出,當(dāng)nc較小時(shí),直觀上較難區(qū)分出含氣泡水是否產(chǎn)生了叢聚現(xiàn)象,因此,在現(xiàn)有的實(shí)驗(yàn)研究中,其氣泡群是否存在叢聚,是無(wú)法從直觀上判斷的[19?21];隨著nc增長(zhǎng),叢聚現(xiàn)象逐漸明顯.nc=30時(shí),從一次實(shí)現(xiàn)上直觀觀測(cè)的結(jié)果而言,其更類似于空化產(chǎn)生的氣泡群形成的空間球狀結(jié)構(gòu),如文獻(xiàn)[30]中圖2.13中的測(cè)量結(jié)果所示.此外,觀察文獻(xiàn)[31]中圖117的空化結(jié)構(gòu),如果空間網(wǎng)狀結(jié)構(gòu)中的細(xì)絲不明顯存在,圖2(d)也可嘗試用于描述網(wǎng)狀空化氣泡群的空間結(jié)構(gòu).因此,本文的理論研究也適用于空化氣泡群結(jié)構(gòu)特征的研究.在此基礎(chǔ)上,分析叢聚含氣泡水對(duì)聲速變化及聲衰減的影響,如圖3所示.

可見,叢聚現(xiàn)象對(duì)含氣泡水內(nèi)聲傳播的影響,主要集中在聲衰減及相速度的峰值附近,對(duì)比無(wú)叢聚含氣泡水可知(氣泡體積分?jǐn)?shù)較小且無(wú)叢聚時(shí),本文的理論模型退化為Commander和Prosperetti[5]的結(jié)論),即使在nc=1時(shí),其也對(duì)含氣泡水的聲衰減有明顯抑制作用,通過對(duì)比實(shí)驗(yàn)結(jié)果,由圖3可以看出,(16)式可以從空間信息上對(duì)現(xiàn)有實(shí)驗(yàn)結(jié)論中的偏差進(jìn)行解釋,即相速度和衰減系數(shù)的峰值往往低于理論預(yù)測(cè)值且存在低頻偏移是由于空間不均勻?qū)е碌?此外,隨nc增加,聲衰減及相速度在峰值附近受到的抑制逐漸加劇,隨頻率遠(yuǎn)離峰值段,叢聚現(xiàn)象對(duì)含氣泡水中聲傳播的影響逐漸變?nèi)?

圖2 不同氣泡數(shù)量密度下叢聚含氣泡水的一次實(shí)現(xiàn)的側(cè)視圖 (a)無(wú)叢聚;(b)nc=1;(c)nc=3;(d)nc=30Fig.2.Realization of clustered bubbly water with different sub-bubble number density:(a)No cluster;(b)nc=1;(b)nc=3;(b)nc=30.

圖3 不同叢聚半徑下聲衰減及相速度隨頻率的變化關(guān)系 (a)聲衰減;(b)相速度Fig.3.Attenuation and phase speed versus frequency with different radius of clusters:(a)Attenuation;(b)phase speed.

3.2 數(shù)值建模方法及統(tǒng)計(jì)分析

基于前文的結(jié)論,能否在實(shí)驗(yàn)中對(duì)空間信息進(jìn)行有效提取,對(duì)含氣泡水內(nèi)聲衰減及相速度的理論預(yù)測(cè)與測(cè)量結(jié)論是否擬合起到?jīng)Q定性作用.考慮到現(xiàn)有含氣泡水實(shí)驗(yàn)研究缺乏相應(yīng)的統(tǒng)計(jì)手段,本文引入了一種比例無(wú)偏估計(jì)以對(duì)后續(xù)實(shí)驗(yàn)研究進(jìn)行指導(dǎo).事實(shí)上,比例無(wú)偏估計(jì)既可以對(duì)(16)式中的叢聚參數(shù)測(cè)量[32],也可以對(duì)統(tǒng)計(jì)量自身進(jìn)行直接測(cè)量[21],本節(jié)引入后一種方法,因?yàn)槠涓哂衅毡樾?在此基礎(chǔ)上,對(duì)統(tǒng)計(jì)結(jié)論與理論預(yù)測(cè)的誤差進(jìn)行觀察.

考慮到統(tǒng)計(jì)全部含氣泡水內(nèi)各個(gè)氣泡的空間位置是不現(xiàn)實(shí)的,我們?cè)诜抡姝h(huán)境中也只采取其樣本(后文統(tǒng)一稱為觀測(cè)窗W)進(jìn)行分析.在仿真環(huán)境中,為保證考慮到W外的叢聚氣泡群對(duì)W的影響,首先在擴(kuò)展觀測(cè)窗W⊕b(0,R)中建立叢聚子氣泡群,一次實(shí)現(xiàn)的叢聚個(gè)數(shù)服從Poisson分布;其次,建立各個(gè)叢聚子氣泡群中的氣泡,每個(gè)子氣泡群的氣泡總數(shù)服從期望為nc的Poisson分布;最后按log-normal分布產(chǎn)生各個(gè)氣泡的半徑,就可以獲得含氣泡水的樣本空間W,如圖2所示,選取邊長(zhǎng)2.6 mm的正方形觀測(cè)窗,以保證觀測(cè)窗內(nèi)樣本氣泡充足(期望值約為2000個(gè)).

為在仿真環(huán)境中建立圖2所示含氣泡水對(duì)聲傳播的影響,須對(duì)其空間信息進(jìn)行提取.根據(jù)(16)式的需求,首先對(duì)其氣泡數(shù)量密度進(jìn)行估計(jì),其標(biāo)準(zhǔn)估計(jì)為

其中N(W)是觀測(cè)窗內(nèi)含氣泡總數(shù),V(W)為觀測(cè)窗體積.該標(biāo)準(zhǔn)估計(jì)是無(wú)偏的.在此基礎(chǔ)上,對(duì)含氣泡水的PCF進(jìn)行統(tǒng)計(jì)分析,本文不加證明地引入比例無(wú)偏的PCF估計(jì)方法(具體證明過程可參考文獻(xiàn)[21]):

該方法通過對(duì)氣泡數(shù)量密度估計(jì)進(jìn)行修正保證了統(tǒng)計(jì)的比率無(wú)偏性(?mre只用于估計(jì)?g(r)).其中

最小,將核函數(shù)表示為

其中h為經(jīng)驗(yàn)參數(shù),為保證估計(jì)結(jié)果平滑且有效信息捕捉準(zhǔn)確,本文選取h=10?5.此外,將修正的氣泡數(shù)量密度寫作

其中V2(W∩?b(ri,r))表示中心位于ri半徑為r的球體在觀測(cè)窗內(nèi)的表面積.根據(jù)(18)式及(19)式對(duì)圖2(c)進(jìn)行估計(jì),可得其PCF如圖4所示.

圖4中紅色誤差棒為PCF理論值±1 MSE,將統(tǒng)計(jì)結(jié)果代入(13)式得到叢聚含氣泡水中聲速及聲衰減隨頻率的變化,如圖5所示.

圖4 PCF統(tǒng)計(jì)結(jié)果及誤差分析Fig.4.Estimated value of PCF and error analysis.

圖5 聲衰減及相速度隨頻率變化的數(shù)值結(jié)果及誤差分析(a)聲衰減;(b)相速度Fig.5.Numerical results of attenuation and phase speed versus frequency comparing with theoretical value:(a)Attenuation;(b)phase speed.

以理論值求得的κ2{nc,R}做對(duì)比分析,可見當(dāng)叢聚產(chǎn)生時(shí),忽略空間信息會(huì)導(dǎo)致對(duì)聲衰減及相速度的預(yù)測(cè)出現(xiàn)較大誤差,而考慮空間信息時(shí)預(yù)測(cè)值與理論值符合良好.如果要進(jìn)一步縮小考慮空間信息時(shí)仍存在的微小誤差,可以通過適當(dāng)擴(kuò)大樣本,或?qū)τ^測(cè)窗進(jìn)行多次采樣、估計(jì)取平均值來達(dá)到目的.

4 結(jié) 論

含氣泡水的空間結(jié)構(gòu)是復(fù)雜的,自然現(xiàn)象產(chǎn)生的含氣泡水中的氣泡很難保證均勻分布.在實(shí)驗(yàn)室環(huán)境下,含氣泡介質(zhì)內(nèi)氣泡間也常存在小范圍的聚集趨勢(shì),而現(xiàn)有描述水下氣泡群對(duì)聲傳播影響的理論往往基于空間均勻分布假設(shè).因此,當(dāng)含氣泡水存在聚集時(shí),現(xiàn)有理論描述聲波在含氣泡水中的傳播時(shí)存在較大誤差.為分析非均勻分布的含氣泡液體對(duì)聲波線性傳播的影響,我們將QCA引入自洽方法中,推導(dǎo)出了考慮空間分布時(shí)多分散介質(zhì)中的等效聲波波數(shù),在此基礎(chǔ)上引入NSP描述了含氣泡水內(nèi)氣泡的叢聚現(xiàn)象.分析發(fā)現(xiàn),即使直觀上無(wú)法觀測(cè)到的叢聚現(xiàn)象,也會(huì)對(duì)含氣泡水內(nèi)聲傳播產(chǎn)生較大干擾,因此含氣泡水的空間信息不可忽視.

通過本文給出的數(shù)值仿真方法,可觀測(cè)到叢聚的含氣泡水對(duì)聲衰減及相速度的準(zhǔn)確預(yù)測(cè)影響較大,忽略空間信息會(huì)過高估計(jì)聲衰減及相速度的峰值,因此,在進(jìn)行相應(yīng)的實(shí)驗(yàn)研究中必須對(duì)空間信息進(jìn)行有效提取,而本文的建模及比例無(wú)偏估計(jì)方法可作為相關(guān)實(shí)驗(yàn)工作的理論基礎(chǔ).