袁迎春，丁德勝

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

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

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

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

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

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

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

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

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

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

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

式中，時間因子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)過十分仔細(xì)而繁復(fù)的推導(dǎo)，得出了n=1，2，3時式（3）的顯式表示，也就是，將上面的積分簡化成幾項Bessel函數(shù)和Struve函數(shù)。而將比值F/ρ0cV=R+iX，即歸一化阻抗函數(shù)，化成這些特殊函數(shù)的計算。以下我們列出這些阻抗函數(shù)公式。

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

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

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

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

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

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

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

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

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

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

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

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

3 計算過程與結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)語

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

最后我們指出，本文提供的這些例子，意味著現(xiàn)在的計算方法，可以作為特殊函數(shù)計算（精度要求不是太高）的補充。計算精度主要取決于方程（8）式的展開系數(shù)。作者希望我們的應(yīng)用數(shù)學(xué)家給出具有更高精度的這類函數(shù)的高斯展開系數(shù)。此外，讀者也許注意到，Bessel函數(shù)的高斯展開，如式（9），是單重級數(shù)形式，而Struve函數(shù)的展開是雙重的形式。直觀上，后者也應(yīng)當(dāng)有類似于式（9）的單重展開形式。有關(guān)Struve函數(shù)這方面的研究結(jié)果，我們將在不久的將來發(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］丁德勝，劉曉峻，黃錦煌.高斯束展開法在計算菲涅爾場積分中的應(yīng)用，中國聲學(xué)進展［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］章力軍，丁德勝.高斯束展開法的注記：指向性和輻射阻抗的簡化計算［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］鮑冬禺，丁德勝.高斯束展開法的注記之二：簡單支撐矩形換能器聲場［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］丁德勝，陸祖宏.活塞類聲場的簡化算法［J］.聲學(xué)學(xué)報，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)形孔徑后的衍射光場［J］.物理學(xué)報，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。