高斯束展開法的注記之三：輻射阻抗的簡化計算*

袁迎春，丁德勝

（1.南京信息職業技術學院電子信息學院，南京 210023；2.東南大學電子科學與工程學院，南京 210096）

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

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

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

1 輻射阻抗及其經典結果［12］

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

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

Greenspan推廣了瑞利的經典結果，研究了活塞振動速度分布具有形式

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

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

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

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

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

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

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

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

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

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

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

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

零階貝塞爾函數展開公式

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

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

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

3 計算過程與結果

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

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

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

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

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

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

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

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

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

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

采用這些表達式編程時，可以少加一些校正項。直接利用（6a）、（6b）式編程時，校正項可能要多一些，應當得出相同結果（或許消去的誤差項更多，精度會更高一些）。讀者不妨一試。

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

4 結語

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

最后我們指出，本文提供的這些例子，意味著現在的計算方法，可以作為特殊函數計算（精度要求不是太高）的補充。計算精度主要取決于方程（8）式的展開系數。作者希望我們的應用數學家給出具有更高精度的這類函數的高斯展開系數。此外，讀者也許注意到，Bessel函數的高斯展開，如式（9），是單重級數形式，而Struve函數的展開是雙重的形式。直觀上，后者也應當有類似于式（9）的單重展開形式。有關Struve函數這方面的研究結果，我們將在不久的將來發表。

［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］丁德勝，劉曉峻，黃錦煌.高斯束展開法在計算菲涅爾場積分中的應用，中國聲學進展［M］.程建春，田靜.北京：科學出版社，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］.聲學學報，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］.物理學報，2015，64（12）：124201.

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

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