高斯束展開(kāi)法的注記之三：輻射阻抗的簡(jiǎn)化計(jì)算*

袁迎春，丁德勝

（1.南京信息職業(yè)技術(shù)學(xué)院電子信息學(xué)院，南京 210023；2.東南大學(xué)電子科學(xué)與工程學(xué)院，南京 210096）

高斯束展開(kāi)方法已廣泛應(yīng)用于快速計(jì)算菲涅爾場(chǎng)積分，即聲場(chǎng)分布的計(jì)算。這種方法的實(shí)質(zhì)是將菲涅爾場(chǎng)積分展開(kāi)為一系列簡(jiǎn)單的基本函數(shù)的疊加，因而復(fù)雜的數(shù)值積分簡(jiǎn)化為一些簡(jiǎn)單函數(shù)（如Gaussian-Laguerre、Guassian-Hermite、高斯函數(shù)等）且項(xiàng)數(shù)不多的計(jì)算，計(jì)算量大為降低［1-6］。有關(guān)這一方法的詳細(xì)描述可以參考綜述文獻(xiàn)［6］。

最近，我們給出高斯函數(shù)展開(kāi)法的一種推廣［7-11］。將文獻(xiàn)中的圓形或矩形函數(shù)的高斯函數(shù)展開(kāi)式或系數(shù)作為已知的結(jié)果，通過(guò)簡(jiǎn)單的數(shù)學(xué)變換，將貝塞爾函數(shù)和一階Struve函數(shù)等特殊函數(shù)表示成高斯函數(shù)或其他簡(jiǎn)單函數(shù)的疊加。利用這一方法，給出了聲學(xué)中一類圓形活塞聲源的指向性函數(shù)以及均勻活塞的輻射阻抗函數(shù)的計(jì)算結(jié)果，與文獻(xiàn)所給的相符合［9］。

本文將這種方法直接應(yīng)用于一類圓形活塞聲源輻射阻抗的計(jì)算，包括均勻分布、邊緣簡(jiǎn)單支撐和邊緣鉗定的情形。所得結(jié)果與直接從特殊函數(shù)計(jì)算所得到的值相比，符合很好。本文是文獻(xiàn)［9］研究結(jié)果的一種直接推廣，作為我們一系列研究報(bào)告的第三部分。

1 輻射阻抗及其經(jīng)典結(jié)果［12］

眾所周知，輻射阻抗是聲學(xué)中一個(gè)重要的參量，即由于聲源振動(dòng)，聲輻射引起的附加于聲源的力阻抗。在電聲器件的設(shè)計(jì)中，除了要知道電聲器件振動(dòng)系統(tǒng)的力學(xué)參數(shù)如質(zhì)量、彈性系數(shù)和力阻外，還必須知道由輻射聲場(chǎng)對(duì)聲源的反作用而產(chǎn)生的附加輻射阻和同振質(zhì)量。

求聲源的輻射阻抗，實(shí)際上即求聲源振動(dòng)時(shí)，媒質(zhì)中的輻射聲場(chǎng)對(duì)聲源的反作用力。除了少數(shù)幾種非常簡(jiǎn)單的聲源，一般情況下，輻射阻抗的計(jì)算涉及到雙重面積積分，通常是振蕩型的積分。瑞利（Lord Rayleigh）采用十分巧妙的方法，將無(wú)窮剛性平面障板中圓形均勻振動(dòng)活塞的輻射阻抗雙重面積積分，簡(jiǎn)化成1階貝塞爾（Bessel）函數(shù)和1階斯特魯夫（Struve）函數(shù)等特殊函數(shù)，已成為聲學(xué)中的經(jīng)典結(jié)果。

Greenspan推廣了瑞利的經(jīng)典結(jié)果，研究了活塞振動(dòng)速度分布具有形式

的聲場(chǎng)輻射問(wèn)題［12］。為便于讀者參考比較，我們盡量采用Greenspan的表述方式和記號(hào)。這里，r0為聲源所在平面z=0上的徑向坐標(biāo)，即源點(diǎn)坐標(biāo)，a是聲源或輻射器的半徑，V為聲源的體積速度，即

n為一整數(shù)，H表示亥維賽階躍函數(shù)。顯然，n=0對(duì)應(yīng)于均勻圓形活塞聲源情形；n=1和n=2分別對(duì)應(yīng)于最低階的邊緣（圓周r0=a處）簡(jiǎn)單支撐和邊緣鉗定活塞情形。當(dāng)r0＞a時(shí)，由式（1）總是有v(r0)=0，這意味著活塞聲源帶有無(wú)限大剛性平面障板。

Greenspan研究的基本出發(fā)點(diǎn)乃是根據(jù)圓形軸對(duì)稱時(shí)諧聲場(chǎng)（隨時(shí)間的變化關(guān)系為eiωt，ω是角頻率）的King速度勢(shì)公式。一旦求出聲場(chǎng)的速度勢(shì)，則相應(yīng)的聲場(chǎng)量，諸如聲壓、質(zhì)點(diǎn)振速分布等皆可隨之而得出。對(duì)于具有式（1）振速分布的聲源，Greenspan得出聲場(chǎng)在其表面的反作用力，表示為

式中，時(shí)間因子eiωt已省略。ρ0是媒質(zhì)（一般是空氣）的密度，m=(u2-k2)12，波數(shù)k=ω/c，c為媒質(zhì)的聲速。公式（3）中以及下文出現(xiàn)的Jn(z)和Hn(z)分別代表第一類n階貝塞爾（Bessel）函數(shù)和n階斯特魯夫（Struve）函數(shù)。Greenspan經(jīng)過(guò)十分仔細(xì)而繁復(fù)的推導(dǎo)，得出了n=1，2，3時(shí)式（3）的顯式表示，也就是，將上面的積分簡(jiǎn)化成幾項(xiàng)Bessel函數(shù)和Struve函數(shù)。而將比值F/ρ0cV=R+iX，即歸一化阻抗函數(shù)，化成這些特殊函數(shù)的計(jì)算。以下我們列出這些阻抗函數(shù)公式。

均勻圓形活塞（n=0）的歸一化的阻函數(shù)為

這些近似公式可以作為本文計(jì)算正確與否以及計(jì)算精度的參照。

通常情況下，高于一階的Struve函數(shù)一般不列表。Greenspan考慮到這一情況，將阻抗函數(shù)中出現(xiàn)的一階以上的貝塞爾（Bessel）函數(shù)和斯特魯夫（Struve）函數(shù)，全部化成零階和一階形式。實(shí)際上直到如今21世紀(jì)，一些常用的數(shù)學(xué)軟件Matlab或Mathcad，以及FORTRAN和C等計(jì)算機(jī)語(yǔ)言軟件包中，仍然不包括斯特魯夫（Struve）函數(shù)的計(jì)算程序或庫(kù)函數(shù)。

由于貝塞爾函數(shù)和斯特魯夫函數(shù)在衍射理論中十分重要，而文獻(xiàn)中有關(guān)這兩種函數(shù)（特別是后者）的計(jì)算方法也不很多，所以我們給出貝塞爾函數(shù)Jn(z)和斯特魯夫函數(shù)Hn(z)的一些近似高斯展開(kāi)表達(dá)式，并用于阻抗函數(shù)的計(jì)算。

2 貝塞爾函數(shù)和斯特魯夫函數(shù)的高斯展開(kāi)式

Wen和Breazeale的一篇文章研究了圓形軸對(duì)稱活塞聲場(chǎng)分布的簡(jiǎn)化計(jì)算方法，即高斯束展開(kāi)法［4-5］。這篇文章的一個(gè)附帶結(jié)果表明：數(shù)學(xué)上，圓形函數(shù)

其中Ak和Bk稱之為展開(kāi)系數(shù)和高斯系數(shù)，可以用計(jì)算機(jī)最優(yōu)化方法來(lái)求得。Wen和Breazeale給出了兩組數(shù)據(jù)，其中一組10項(xiàng)的數(shù)據(jù)列在文獻(xiàn)［4］的表1中；文獻(xiàn)［13-14］給出了另外兩組數(shù)據(jù)（項(xiàng)數(shù)稍多，精度稍高，一組25項(xiàng)［13］，一組 15項(xiàng)［14］）。這個(gè)近似展開(kāi)式（8）在衍射理論中十分重要，可以簡(jiǎn)化許多復(fù)雜問(wèn)題的分析和計(jì)算［15-20］。

以下我們將根據(jù)展開(kāi)式（8），列出貝塞爾函數(shù)和斯特魯夫（Struve）函數(shù)的高斯展開(kāi)式。

零階貝塞爾函數(shù)展開(kāi)公式

公式（10）的導(dǎo)出，階數(shù)ν要滿足 Re(ν)＞1/2 的條件，因而直接從公式（10），取ν=0，得出 H0(x)的展開(kāi)式，可能會(huì)有問(wèn)題。利用遞推公式

這里要注意的是，上列各公式中系數(shù)Ak和Bk與方程（8）為同一組數(shù)據(jù)。有關(guān)這些公式右邊的展開(kāi)精度，可以參考文獻(xiàn)［7］。文獻(xiàn)［7］中圖2給出了一個(gè)比較［n=0 ，公式（9a）］。用10項(xiàng)高斯展開(kāi)系數(shù)［4］計(jì)算，在大約0-20區(qū)間上，與J0符合很好，相對(duì)誤差大約1%～2%。而另外一組15項(xiàng)的展開(kāi)系數(shù)［5］在大約0～30的范圍內(nèi)，可以更好地?cái)M合J0。上面這些公式的推導(dǎo)細(xì)節(jié)和和精度等的詳細(xì)分析，這里就不給出了。

以下我們用這些特殊函數(shù)的近似展開(kāi)式來(lái)計(jì)算圓形活塞阻抗函數(shù)。

3 計(jì)算過(guò)程與結(jié)果

這里我們給出均勻活塞（n=0）、邊緣簡(jiǎn)單支撐（n=1）和邊緣鉗定這三種聲源分布的輻射阻抗函數(shù)。最近楊某［14］給出了另外一組15項(xiàng)高斯展開(kāi)系數(shù)，精度更高一些。以下計(jì)算中，主要用Wen和Breazeale的10項(xiàng)高斯展開(kāi)系數(shù)［4］（第1組）和楊某［14］的這組數(shù)據(jù)（第2組），順便比較一下精度。所用計(jì)算軟件為Matlab 6.1。

對(duì)于均勻活塞情形，直接將貝塞爾函數(shù)和斯特魯夫函數(shù)的展開(kāi)式（10）等代入阻函數(shù)（4a）和抗函數(shù)公式（4b）中計(jì)算即可，圖1給出了均勻活塞源歸一化阻函數(shù)和抗函數(shù)的計(jì)算結(jié)果。

圖1 均勻活塞（n=0）的輻射阻抗函數(shù)

圖1中R0(2ka)分別采用近似展開(kāi)（9b）和直接計(jì)算貝塞爾函數(shù)［Matlab中函數(shù)名為besselj（mu，x）］所得，兩種方法所得結(jié)果頗相符合［第1組數(shù)據(jù)所得結(jié)果的最大絕對(duì)誤差約9×10-3，第2組的最大誤差約為5×10-4］。關(guān)于抗函數(shù)X0(2ka)，我們把根據(jù)（10a）式的計(jì)算結(jié)果，與課本中x=0(0.5)20的數(shù)值相比較，相對(duì)誤差大約1%～2%［9］。

對(duì)于n=1，2的情況，編程計(jì)算時(shí)要有些技巧。將公式（10a）、（12）代入（5），即得簡(jiǎn)單支撐活塞的抗函數(shù)近似計(jì)算公式

直接根據(jù)式（13）計(jì)算得到歸一化抗函數(shù)曲線，如圖2所示。對(duì)比Greenspan文獻(xiàn)中圖5，可見(jiàn)符合很好。但是，在計(jì)算阻函數(shù)R1過(guò)程中，出現(xiàn)了意外。我們把第1組數(shù)據(jù)［4］（Wen和Breazeale的10項(xiàng)系數(shù)）代入到式（5a）中，計(jì)算時(shí)發(fā)現(xiàn)：當(dāng)x=2ka較大時(shí)（x約大于4），所得結(jié)果，與直接調(diào)用Matlab中的Bessel函數(shù)的計(jì)算結(jié)果，符合很好；而當(dāng)x=2ka較小時(shí)，卻完全不符。實(shí)際上當(dāng)x很小時(shí)（接近于0，如取0.01），計(jì)算值是發(fā)散的。用第2組數(shù)據(jù)計(jì)算，同樣有此錯(cuò)誤。而根據(jù)阻抗函數(shù)的近似公式（5c），當(dāng)x較小時(shí)，R1(x)≈x2/8。將公式（5b）改為如下形式后

情況變得十分明朗。原來(lái)問(wèn)題出在上面式（14）的最后一項(xiàng)中。當(dāng)J0(x)小近似展開(kāi)后，最后一項(xiàng)中分子中的常數(shù)項(xiàng)應(yīng)當(dāng)為0。而將近似式（9a）代入（14）后，所得常數(shù)項(xiàng)δ=-1實(shí)際上并不等于零（這個(gè)值大致取決于（8）式的近似程度，第1組數(shù)據(jù)給出值約為-0.0127，第2組數(shù)約為0.011 5），那么式（14）主要計(jì)算誤差由最后一項(xiàng)產(chǎn)生，也就是δ x2。若x=0.01，此時(shí)的誤差就會(huì)被放大大約104倍。難怪計(jì)算結(jié)果發(fā)散!編程時(shí)，我們只需加上一些校正項(xiàng)，可以簡(jiǎn)單地消除這些誤差，計(jì)算結(jié)果如圖2所示。

圖2 簡(jiǎn)單支撐活塞（n=1）的輻射阻抗函數(shù)

Greenspan實(shí)際上已經(jīng)注意到類似的問(wèn)題。也許他那個(gè)時(shí)代（1980年代），Bessel函數(shù)等特殊函數(shù)的計(jì)算精度不是太高，他的文章［12］中沒(méi)有給出x＜1的那部分曲線。

采用類似的步驟，我們計(jì)算了邊緣鉗定活塞聲源的輻射阻抗函數(shù)，如圖3所示。計(jì)算過(guò)程中，我們利用貝塞爾函數(shù)和斯特魯夫函數(shù)的遞推公式，分別將式（6a）和（6b）改寫成下面形式

采用這些表達(dá)式編程時(shí)，可以少加一些校正項(xiàng)。直接利用（6a）、（6b）式編程時(shí)，校正項(xiàng)可能要多一些，應(yīng)當(dāng)?shù)贸鱿嗤Y(jié)果（或許消去的誤差項(xiàng)更多，精度會(huì)更高一些）。讀者不妨一試。

圖3 邊緣鉗定活塞（n=2）的輻射阻抗函數(shù)

4 結(jié)語(yǔ)

我們給出了高斯函數(shù)展開(kāi)法的推廣，將圓形函數(shù)的高斯展開(kāi)作為已知結(jié)果，通過(guò)簡(jiǎn)單的數(shù)學(xué)變換，將貝塞爾函數(shù)和Struve函數(shù)表示成高斯函數(shù)的近似和。計(jì)算了聲學(xué)中一類活塞聲源的輻射阻抗函數(shù)，與直接計(jì)算特殊函數(shù)的結(jié)果相比，我們的方法給出了相當(dāng)一致的結(jié)果。這種方法還可以直接應(yīng)用于活塞聲源的輻射功率的計(jì)算。

最后我們指出，本文提供的這些例子，意味著現(xiàn)在的計(jì)算方法，可以作為特殊函數(shù)計(jì)算（精度要求不是太高）的補(bǔ)充。計(jì)算精度主要取決于方程（8）式的展開(kāi)系數(shù)。作者希望我們的應(yīng)用數(shù)學(xué)家給出具有更高精度的這類函數(shù)的高斯展開(kāi)系數(shù)。此外，讀者也許注意到，Bessel函數(shù)的高斯展開(kāi)，如式（9），是單重級(jí)數(shù)形式，而Struve函數(shù)的展開(kāi)是雙重的形式。直觀上，后者也應(yīng)當(dāng)有類似于式（9）的單重展開(kāi)形式。有關(guān)Struve函數(shù)這方面的研究結(jié)果，我們將在不久的將來(lái)發(fā)表。

［1］Cook B D，Arnoult III W J.Gaussian-Laguerre/Hermite Formula?tion for the Nearfield of an Ultrasonic Transducer［J］.J Acoust Soc Amer，1976，59：9.

［2］Cavanagh E，Cook B D.Gaussian-Laguerre Description of Ultra?sonic Fields-Numerical Example：Circular Piston［J］.J Acoust Soc Amer，1980，67：1136.

［3］Thompson R B，Gray T A，J，et al.The radiation of Elliptical and Bicylindrically Focused Piston Transducers［J］.J Acoust Soc Am?er，1987，82：1818.

［4］Wen J J，Breazeale M A.A Diffraction Beam Field Expressed as the Superposition of Gaussian Beams［J］.J Acoust Soc Amer，1988，83：1752.

［5］Huang D，Breazeale M A.A Gaussian Finite-Element Method for Description of Sound Diffraction［J］.J Acoust Soc Am，1999，106：1771-1781.

［6］丁德勝，劉曉峻，黃錦煌.高斯束展開(kāi)法在計(jì)算菲涅爾場(chǎng)積分中的應(yīng)用，中國(guó)聲學(xué)進(jìn)展［M］.程建春，田靜.北京：科學(xué)出版社，2008：55-68.

［7］Ding D S，Tong X J，He P Z.Supplementary Notes on the Gauss?ian Beam Expansion［J］.J Acoust Soc Am，2005，118：608.

［8］Dai Yurong，Ding Desheng.Further Notes on the Gaussian Beam Expansion［J］.Chinese Physics Letters，2012，29：024301.

［9］章力軍，丁德勝.高斯束展開(kāi)法的注記：指向性和輻射阻抗的簡(jiǎn)化計(jì)算［J］.電子器件，2013，36（6）：789-792.

［10］Ding D，Lu H，Shen C.A Novel Algorithm for the Sound Field of El?liptically Shaped Transducers［J］.Chin Phys Lett，2014，31：64301.

［11］鮑冬禺，丁德勝.高斯束展開(kāi)法的注記之二：簡(jiǎn)單支撐矩形換能器聲場(chǎng)［J］.電子器件，2015，38（2）：278-282.

［12］Greenspan M.Piston radiator：Some Extensions of the Theory［J］.J Acoust Soc Am，1979，65：608-621.

［13］Kim H J，Schmerr L W，Sedov A.Generation of the Basis Sets for Multi-Gaussian Ultrasonic Beam Models—An Overview［J］.J Acoust Soc Amer，2006，119：1971-1978.

［14］Liu W，Yang J.A Simple and Accurate Method for Calculating the Gaussian Beam Expansion Coefficients［J］.Chin Phys Lett，2010，27：124301.

［15］丁德勝，陸祖宏.活塞類聲場(chǎng)的簡(jiǎn)化算法［J］.聲學(xué)學(xué)報(bào)，1996，21（增刊）：421.

［16］Ding D S，Lu Z H.A Simplified Method to Calculate the Sound Field of Pistonlike Source［J］.Chin J Acoust，1996，15：213.

［17］Ding D S，Liu X J.Approximate Description for Bessel，Bessel-Gauss and Gaussian Beams with Finite Aperture［J］.J Opt Soc Amer，1999，A 16：1286.

［18］Zhang Y，Liu J Q，Ding D S.Sound Field Calculations of Ellipti?cal Pistons by the Superposition of Two-Dimensional Gaussian Beams［J］.Chin Phys Lett，2002，19：1825.

［19］Ding D S，Zhang Y.A Simple Calculation Approach for the Sec?ond-Harmonic Sound Beam Generated by an Arbitrary Distribu?tion Source［J］.Chin Phys Lett，2004，21：503.

［20］謝曉霞，王碩琛，吳逢鐵.Bessel光束經(jīng)橢圓環(huán)形孔徑后的衍射光場(chǎng)［J］.物理學(xué)報(bào)，2015，64（12）：124201.

袁迎春（1981-），女，碩士，講師，畢業(yè)于東南大學(xué)，現(xiàn)任教于南京信息職業(yè)技術(shù)學(xué)院，電子信息學(xué)院。主要從事微波理論與技術(shù)研究；

丁德勝（1965-），男，江蘇人。理學(xué)博士，東南大學(xué)電子科學(xué)與工程學(xué)院教授、博士生導(dǎo)師。主要從事聲學(xué)和微波技術(shù)的研究及教學(xué)，dds@seu.edu.cn。