增材制造技術(shù)中的流固耦合度影響因素

摘? ?要： 對流固界面的耦合程度為增材制造技術(shù)帶來的挑戰(zhàn)進行報道。以流體柱體外圍繞無限遠各向同性固體模型為基準，在角頻率—波數(shù)域建立流體滿足的聲波方程、固體滿足的彈性波方程和流固耦合邊界條件。引入流固耦合度這一概念，導(dǎo)出法向位移斷續(xù)、法向應(yīng)力斷續(xù)、切向應(yīng)力為零的流固界面不完全耦合邊界條件，考察不同中心頻率的單極子聲源在模型中激發(fā)的頻散特征和波形特征。流體柱體與硬性固體相耦合，降低了各階模式的截止頻率，壓制了各階模式的傳播速度，法向位移耦合度或法向應(yīng)力耦合度的增大，分別使得各階模式更加遠離或靠近流體柱體外圍不存在硬性固體的情形;流體柱體與軟性固體相耦合，頻散波僅存在一個具有較低截止頻率的0階模式，法向位移耦合度或法向應(yīng)力耦合度的減小，均使得0階模式更加遠離流體柱體外圍不存在軟性固體的情形，且截止頻率持續(xù)降低直至消失，甚至可能導(dǎo)致0階模式發(fā)生反轉(zhuǎn)。法向應(yīng)力耦合度為零時的復(fù)雜波串并非頻散曲線中的某階模式，而是某種其他非頻散波。留存于流體柱體內(nèi)的“滯能縱波”和近乎于流體柱體外圍繞橫向各向同性固體時的傳播特性，是流體柱體與軟性固體相耦合時產(chǎn)生的特殊數(shù)學(xué)物理現(xiàn)象。

關(guān)鍵詞： 增材制造;流固耦合度;角頻率—波數(shù)域;流固界面不完全耦合邊界條件;單極子聲源;頻散特征;截止頻率;0階模式;反轉(zhuǎn);非頻散波;“滯能縱波”;橫向各向同性

引言

聲波、彈性波的傳播特性，是評價增材制造工藝效果的基礎(chǔ)指標，是檢驗增材制造技術(shù)制備所得物品性能的得力工具[1]。將聲波、彈性波的傳播特性研究貫穿于增材制造全生命周期始終，將有力推動增材制造行業(yè)精益生產(chǎn)，促進增材制造產(chǎn)品質(zhì)量和競爭力提升[2]。

2017年，我國提出行動計劃，明確要求提升增材制造裝備、核心器件及軟件質(zhì)量[3]，其中數(shù)學(xué)物理方法和算法是支撐增材制造工藝智能規(guī)劃、在線檢測、過程控制等軟件功能的核心[4]。聲波、彈性波傳播的數(shù)學(xué)物理研究方法涵蓋理論法和數(shù)值法兩種，理論法包括基于并矢格林函數(shù)的求解方法、角頻率—波數(shù)域分析法等，數(shù)值法包括各時空域的有限差分法、有限元法等。文獻[5]采用融合了完全匹配層全吸收邊界的、適用于斜向各向同性（屬于一種各向異性）介質(zhì)的時域交錯網(wǎng)格有限差分法，對彈性材料的各向異性為增材制造技術(shù)帶來的挑戰(zhàn)進行了報道，呼吁科研人員重視各向異性介質(zhì)的Thomsen參數(shù)大小以及介質(zhì)在各個方位的傳播速度差異程度對彈性波傳播特性的影響;文獻[6]采用彈性波傳播的并矢格林函數(shù)解，以增材制造工藝中常見的負泊松比超材料為研究對象，初步探索了超材料彈性參數(shù)線性化的普適方法，為增材制造領(lǐng)域共性技術(shù)創(chuàng)新提供了驅(qū)動力;文獻[7]采用頻域有限元方法，以1.5維聲波、彈性波方程為例，對內(nèi)在機理相對復(fù)雜、細節(jié)把控要求較高、求解具有一定技術(shù)門檻的半整數(shù)維波動方程進行了考察，尤其是對完全匹配層拉伸尺度選取、等效積分弱形式推導(dǎo)、流固耦合邊界條件設(shè)置、極點規(guī)避等技術(shù)細節(jié)進行了注記，為更高維空間波傳播特性的研究打開了局面;文獻[8-9]基于聲波和彈性波在流固界面的雙向轉(zhuǎn)換機制，推導(dǎo)了應(yīng)用于時域有限元方法的，帶有完全匹配層的聲波—彈性波耦合（即流固耦合）方程的等效積分弱形式，并立足于增材制造技術(shù)應(yīng)用需求，構(gòu)思了若干聲波—彈性波耦合模型，研討了遠場、近場條件下聲波—彈性波耦合對波場傳播的影響，得到了大量具有實際價值的結(jié)論和啟示。文獻[10]為流固耦合問題提出了一種來源于經(jīng)驗的傳遞矩陣校正方法，為流固耦合模型頻散曲線的智能求取提供了便利。上述文獻雖然覆蓋了各維度、各時空域下基于各種數(shù)學(xué)物理方法和算法的聲波、彈性波及二者之間相互耦合機制的研究，但并未考慮到聲波、彈性波之間的耦合程度對波場傳播的影響，而增材制造工藝的固有特點又是導(dǎo)致這種耦合程度的不確定性的根源所在。據(jù)此，本文引入流固耦合度這一概念，在角頻率—波數(shù)域考察不同流固耦合度下波場傳播的頻散特征和波形特征，以期為增材制造工藝效果評價和物品制備性能檢驗提供更加嚴密、嚴謹?shù)睦碚摶A(chǔ)。

1? 流體柱體中聲波的傳播

式（40）稱作流固界面不完全耦合邊界條件。其中，、分別稱作法向位移耦合度、法向應(yīng)力耦合度，處于區(qū)間[0， 1];稱作切向應(yīng)力耦合度，僅取0或1。當、、均為0時，式（40）退化為第1章所討論的流體柱體中聲波的傳播情形;當三者均為1時，式（40）退化為第2章所討論的流體柱體外圍繞無限遠固體，且流固界面完全耦合時，流體柱體內(nèi)聲波的傳播情形;對于、（二者可以不相等）且的情形，本文將邊界條件描述為：法向位移斷續(xù)、法向應(yīng)力斷續(xù)、切向應(yīng)力為零。

4? 應(yīng)用與討論

4.1? 模型與參數(shù)

將第1章所述流體柱體的情形記作模型1，將第2章所述流體柱體外圍繞無限遠固體的情形記作模型2。限于篇幅，本文僅對單極子聲源激發(fā)的波場進行考察。多極子聲源情形的研究結(jié)論必然有所不同，但研究思路類似。

模型基本參數(shù)如表1所示，其中硬性固體指橫波速度大于相鄰流體中聲波傳播速度的固體，軟性固體反之。

4.2? 頻散特征計算與討論

所有滿足式（46）的點集構(gòu)成模型2的頻散曲線。對于式（46），可根據(jù)文獻[10]，首先在給定的波源頻率和軸向相速度范圍內(nèi)對二維譜進行大體觀察，判斷二維譜的極性翻轉(zhuǎn)特性，然后針對不同表現(xiàn)形式的極性翻轉(zhuǎn)特性，采用校正因子實施“反翻轉(zhuǎn)”，最后根據(jù)“反翻轉(zhuǎn)”過程整理得到傳遞矩陣校正圖譜，從而通過牛頓迭代法[11]求解。牛頓迭代法詳見附錄B。

4.2.1? 流固界面完全耦合

首先，在模型2的流體柱體外圍選用硬性固體，且將流固界面設(shè)定為完全耦合（即），繪制頻散曲線，并與模型1相對比，如圖1a所示。對比得知，在0～25 kHz頻段，兩個模型均具有0階、1階、2階模式，其中0階模式無截止頻率，其他階模式具有較高的截止頻率。此外，流體柱體與硬性固體相耦合，不僅降低了各階模式的截止頻率，而且壓制了各階模式的傳播速度。在模型2中，1階以上（含）模式在截止頻率下傳播速度最大，等于硬性固體的橫波速度，而在模型1中，1階以上（含）模式在截止頻率下傳播速度無限大（事實上，模型1的流體柱體模型是一個典型的“剛性壁”波導(dǎo)問題，其等效于流體柱體外圍繞無限遠“波速無限大的超硬性固體”的情形）。

其次，在模型2的流體柱體外圍選用軟性固體，仍將流固界面設(shè)定為完全耦合，繪制頻散曲線，并與模型1相對比，如圖1b所示。與圖1a截然不同的是，模型2的頻散曲線似乎與模型1毫無相關(guān)性。在0～25 kHz頻段對模型2的頻散曲線進行求解，僅可得到0階模式，且該模式具有一個較低的截止頻率（<1 kHz）。0階模式在截止頻率下傳播速度最大，等于軟性固體的橫波速度。

4.2.2? 流固界面不完全耦合

令人出乎意料的是，對于任意滿足式（47）的流固耦合度，無論模型2采用硬性固體還是軟性固體，計算得到的頻散曲線均與流固界面完全耦合的情形重合。

進一步，考察法向位移耦合度、法向應(yīng)力耦合度不同的情形。不失一般性，以下僅討論0階、1階模式，并且為便于對比，模型1的上述模式也將在頻散特征分析圖中用虛線繪出。

（1）硬性固體

① 法向位移耦合度不變，法向應(yīng)力耦合度變化

令不變，將分別取值為0.25、0.5、0.75、1，計算結(jié)果如圖2a所示。在的不同取值下，計算結(jié)果完全相同，其中模型1和模型2的0階模式發(fā)生了重合，1階模式發(fā)生了交叉。與圖1a對比得知，頻散曲線交叉現(xiàn)象似乎是由模型2的1階以上（含）模式之間相互交融而引起的。

令不變，將分別取值為0、0.25、0.5、0.75、1，計算結(jié)果如圖2b所示。當時，模型2的0階模式不存在，1階模式在截止頻率下傳播速度無限大。隨著的增大，模型2的0階模式開始顯現(xiàn)，并不斷趨近于模型1的0階模式。在的不同取值下，模型2的1階模式近乎巧合地交于一點，當實際頻率一定，且大于這一巧合交點所在的頻率時，傳播速度隨著的增大而增大，反之亦然。初步推斷這一巧合交點所在的頻率為模型1的1階模式的截止頻率，該截止頻率有理論解[12]，即

其中，為截止頻率，為第一類變形Bessel函數(shù)的第一個虛部大于0的根[13]。將相關(guān)參數(shù)代入式（48），求得。采用牛頓迭代法在該截止頻率附近對的不同取值進行精算，發(fā)現(xiàn)確鑿如此。

分別令、、不變，取值同上，計算結(jié)果分別如圖2c、2d、2e所示。除數(shù)值必然不同外，上述結(jié)果與時的規(guī)律相似，其中一個顯著的區(qū)別是，當一定時， 的增大，使得模型2的0階模式更加遠離模型1。

② 法向應(yīng)力耦合度不變，法向位移耦合度變化

限于篇幅，僅討論的情形。將分別取值為0、0.25、0.5、0.75、1，計算結(jié)果如圖2f所示。顯而易見，當一定時，的增大，除了使得模型2的0階模式更加遠離模型1以外，也使得模型2的1階模式更加稍微遠離模型1。

（2）軟性固體

章節(jié)4.2.1揭示了流體柱體外圍繞無限遠軟性固體時頻散特性的捉摸不定。因此，以下從流固耦合度由大至小開展試探性分析。

① 法向位移耦合度不變，法向應(yīng)力耦合度變化

令不變，將分別取值為1、0.75、0.5、0.25、0，計算結(jié)果如圖3a所示。剛減小到0.75，模型2的0階模式的截止頻率即消失;當減小到0.25時，模型2的0階模式甚至發(fā)生了反轉(zhuǎn);當時，模型2的0階模式不復(fù)存在。

令不變，取值同上，計算結(jié)果如圖3b所示。與時不同，當減小到0.5時，模型2的0階模式的截止頻率方消失;和時的頻散特性與時相似。

令不變，取值同上，計算結(jié)果如圖3c所示。與以上情形不同，當減小到0.25時，模型2的0階模式的截止頻率方消失，且模型2的0階模式未再發(fā)生反轉(zhuǎn)。

令不變，取值同上，計算結(jié)果如圖3d所示。在的不同取值下，模型2的0階模式的截止頻率均未消失，且模型2的0階模式均未發(fā)生反轉(zhuǎn)。此外，模型2的0階模式對的敏感度已變得非常弱。

令不變，將分別取值為1、0.75、0.5、0.25，計算結(jié)果如圖3e所示。在的不同取值下，計算結(jié)果完全相同，模型2的0階模式頻散極其輕微，且其截止頻率在以上情形中最大（≈3 kHz）。

② 法向應(yīng)力耦合度不變，法向位移耦合度變化

限于篇幅，僅討論的情形。將分別取值為0.25、0.5、0.75、1，計算結(jié)果如圖3f所示。顯而易見，當一定時，的增大，使得模型2的0階模式更加接近流固界面完全耦合的情形。

4.3? 波形特征計算與討論

結(jié)合章節(jié)4.2計算得到的頻散特征，考察代表性中心頻率3 kHz和10 kHz的單極子聲源在模型1和模型2中激發(fā)的波形。聲源與接收器之間的距離取1.0～2.0 m，接收器之間的間隔取0.2 m。

4.3.1? 流固界面完全耦合

（1）硬性固體

繪制中心頻率為3 kHz的單極子聲源在模型1中激發(fā)的波形，并與其中直達波的貢獻相對比，如圖4a所示。在模型1中，直達波的貢獻較為微弱，且衰減較快;反射波占據(jù)主要地位，為純粹的0階模式。進一步，繪制該聲源在模型2中激發(fā)的波形，并與模型1相對比，如圖4b所示。流體柱體外圍硬性固體的存在，使得流體柱體內(nèi)的模式波發(fā)生了遲滯，即傳播速度降低，這與圖1a所示的頻散曲線相符合。

同理，繪制中心頻率為10 kHz的單極子聲源在模型1中激發(fā)的波形，并與其中直達波的貢獻相對比，如圖5a所示。在模型1中，直達波的貢獻較為微弱，且衰減更快，但直達波與反射波的幅度比大于聲源中心頻率為3 kHz的情形。盡管如此，反射波依舊占據(jù)主要地位，且在一個0階模式脈沖后有復(fù)雜波串緊隨，這些波串很容易被認作1階模式。進一步，繪制該聲源在模型2中激發(fā)的波形，如圖5b所示。如同章節(jié)4.2.1所提到的，流體柱體外圍硬性固體的存在，壓制了各階模式的傳播速度，充其量不會超過硬性固體的橫波速度，這使得模式波特征得到了很大程度的簡化。事實上，圖5b的波串中也蘊含著多種模式，且以非頻散波為主，但這已超出本文的討論范疇，感興趣的讀者可與筆者私下討論。

（2）軟性固體

繪制中心頻率為3 kHz的單極子聲源在模型2中激發(fā)的波形，并與模型1相對比，如圖6a所示。流體柱體外圍軟性固體的存在，使得流體柱體內(nèi)的模式波受到了對沖，幅度嚴重衰減，其他模式波的幅度也遠遠不及流體柱體外圍不存在軟性固體的情形。事實上，模型2的波串中蘊含著兩種模式，一是傳播速度較快的，以固體縱波速度傳播的非頻散波，未在圖1b中展示;二是傳播速度較慢的，以低于固體橫波速度傳播的頻散波，即圖1b中的0階模式。

同理，繪制中心頻率為10 kHz的單極子聲源在模型2中激發(fā)的波形，如圖6b所示。此時，波串中仍蘊含著兩種模式，一是幅度較大、傳播速度較快的，以固體縱波速度傳播的非頻散波;二是衰減較嚴重的，與上述非頻散波相交融的，以更加稍微低于固體橫波速度傳播的頻散波。此外，將模型2中的波形與模型1（即圖5a中實線）相對比，發(fā)現(xiàn)流體柱體外圍繞軟性固體時產(chǎn)生的非頻散波，似乎是將流體柱體外圍不存在軟性固體時產(chǎn)生的反射波的一部分釋放出去的結(jié)果，這使得留存在流體柱體內(nèi)的這一部分持續(xù)了約1 ms左右即趨于平穩(wěn)。部分文獻將這一非頻散波稱作“漏能縱波”（Leaky P-wave），但按照上述分析，將留存在流體柱體內(nèi)的這一部分稱作“滯能縱波”（Hysteresis P-wave）似乎更為貼切。

4.3.2? 流固界面不完全耦合

首先考察法向位移耦合度、法向應(yīng)力耦合度相同的情形。與章節(jié)4.2.2的分析相一致，對于任意滿足式（47）的流固耦合度，無論模型2采用硬性固體還是軟性固體，計算得到的波形均與流固界面完全耦合的情形重合。進一步，考察法向位移耦合度、法向應(yīng)力耦合度不同的情形。

（1）硬性固體

令不變，將分別取值為0.5、1，繪制中心頻率為3 kHz的單極子聲源激發(fā)的波形，計算結(jié)果與模型1完全相同。進一步，將單極子聲源中心頻率調(diào)整至10 kHz，依舊如是。

令不變，將分別取值為0、0.5、1，繪制中心頻率為3 kHz的單極子聲源激發(fā)的波形，計算結(jié)果對比圖如圖7a和7b所示。當時，復(fù)雜波串再次顯現(xiàn)，根據(jù)圖2c，3 kHz的中心頻率遠遠小于模型2的1階模式的截止頻率，排除了這一復(fù)雜波串為1階模式的可能，它可能是0階模式或某種其他非頻散波。隨著的增大，上述復(fù)雜波串消失，此外流體柱體內(nèi)的模式波傳播速度相應(yīng)提高，這與圖2c所示的頻散曲線相符合。進一步，將單極子聲源中心頻率調(diào)整至10 kHz，計算結(jié)果對比圖如圖7c和7d所示。當時，復(fù)雜波串再次顯現(xiàn)，但由于模型2的0階模式并不存在，進一步排除了這一復(fù)雜波串為0階模式的可能，它只能是某種其他非頻散波。與聲源中心頻率為3 kHz的情形相似，隨著的增大，上述復(fù)雜波串消失，此外流體柱體內(nèi)的模式波傳播速度相應(yīng)提高，這也與圖2c所示的頻散曲線相符合。但可意外發(fā)現(xiàn)，當時，波串持續(xù)時間相對更長。

令不變，將分別取值為0、0.5、1，繪制中心頻率為3 kHz、10 kHz的單極子聲源激發(fā)的波形，分析結(jié)論與時基本相同。限于篇幅，計算結(jié)果不再展示。

（2）軟性固體

令不變，將分別取值為0.5、1，繪制中心頻率為3 kHz的單極子聲源激發(fā)的波形，計算結(jié)果與模型1完全相同。進一步，將單極子聲源中心頻率調(diào)整至10 kHz，依舊如是。

令不變，將分別取值為0、0.5、1，繪制中心頻率為3 kHz的單極子聲源激發(fā)的波形，計算結(jié)果對比圖如圖8a和8b所示。當時，復(fù)雜波串再次顯現(xiàn)，只是比硬性固體的情形更規(guī)則有序一些，同樣根據(jù)圖3c，0階模式并不存在，因此這一復(fù)雜波串只能是某種其他非頻散波。隨著的增大，上述復(fù)雜波串消失，此外流體柱體內(nèi)的模式波傳播速度相應(yīng)提高，這與圖3c所示的頻散曲線相符合，但可以驚奇地發(fā)現(xiàn)，、時的模式波竟至于體現(xiàn)為近乎于流體柱體外圍繞無限遠橫向各向同性（屬于一種各向異性）固體時流體柱體內(nèi)聲波的傳播特性！這極可能造成對流體柱體外圍固體各向異性程度的誤判。進一步，將單極子聲源中心頻率調(diào)整至10 kHz，計算結(jié)果對比圖如圖8c和8d所示。當時，復(fù)雜波串再次顯現(xiàn)，與前文分析相同，其只能是某種其他非頻散波。與聲源中心頻率為3 kHz的情形相似，隨著的增大，上述復(fù)雜波串消失，此外流體柱體內(nèi)的模式波傳播速度相應(yīng)提高，這也與圖3c所示的頻散曲線相符合，且仍可發(fā)現(xiàn)，當時，波串持續(xù)時間相對更長。

令不變，將分別取值為0、0.5、1，繪制中心頻率為3 kHz、10 kHz的單極子聲源激發(fā)的波形，分析結(jié)論與時基本相同。限于篇幅，計算結(jié)果不再展示。

5? 結(jié)論與局限

本文以流體柱體外圍繞無限遠固體模型為基準，引入流固耦合度這一概念，在角頻率—波數(shù)域考察了不同流固耦合度下波場傳播的頻散特征和波形特征，為增材制造工藝效果評價和物品制備性能檢驗提供了更加嚴密、嚴謹?shù)睦碚摶A(chǔ)。主要結(jié)論有：

（1） 流體柱體與硬性固體相耦合，不僅降低了各階模式的截止頻率，而且壓制了各階模式的傳播速度。法向位移耦合度的增大，使得各階模式更加遠離流體柱體外圍不存在硬性固體的情形;法向應(yīng)力耦合度的增大，使得各階模式更加靠近流體柱體外圍不存在硬性固體的情形。

（2） 流體柱體與軟性固體相耦合，頻散波僅存在0階模式，且該模式具有一個較低的截止頻率。法向位移耦合度或法向應(yīng)力耦合度的減小，均使得0階模式更加遠離流體柱體外圍不存在軟性固體的情形，且截止頻率持續(xù)降低直至消失，甚至可能導(dǎo)致0階模式發(fā)生反轉(zhuǎn)。

（3） 流體柱體無論是與硬性固體還是軟性固體相耦合，法向應(yīng)力耦合度為0時的復(fù)雜波串均并非頻散曲線中的某階模式，而是某種其他非頻散波。

（4） 留存于流體柱體內(nèi)的“滯能縱波”和近乎于流體柱體外圍繞橫向各向同性固體時的傳播特性，是流體柱體與軟性固體相耦合時產(chǎn)生的特殊數(shù)學(xué)物理現(xiàn)象。

本研究盡管力求全面，但限于篇幅，仍存在一定的局限，包括但不限于：

（1） 當法向位移耦合度和法向應(yīng)力耦合度相同且非零時，頻散特征和波形特征并未顯現(xiàn)出差別，表明流固耦合度的完備性有待加強;

（2） 多極子聲源，如偶極子、四極子聲源在流固耦合模型中的傳播特性有待分析;

（3） 流體與各向異性固體相耦合時的波場傳播特性有待考察;

（4） 固體—固體界面的彈性耦合度勢必比流固耦合度更為復(fù)雜，有待深入;

（5） 為了有效檢驗由增材制造技術(shù)制備所得的更為復(fù)雜的物品的性能，流固耦合度、彈性耦合度在數(shù)值法，如各時空域的有限差分法、有限元法中的數(shù)學(xué)物理表達，有待探究。

參考文獻

[1] 邵長金， 楊振清， 周廣剛， 等. 場與波[M]. 東營： 中國石油大學(xué)出版社， 2015.

[2] 彼得 S. 潘迪， 羅伯特 P. 紐曼. 六西格瑪管理法： 世界頂級企業(yè)追求卓越之道： 第2版[M]. 北京： 機械工業(yè)出版社， 2017.

[3] 增材制造產(chǎn)業(yè)發(fā)展行動計劃（2017-2020年）[R].

[4] 顧樵. 數(shù)學(xué)物理方法[M]. 北京： 科學(xué)出版社， 2019.

[5] 張闊， 劉鶴. 增材制造技術(shù)中的彈性各向異性影響因素[J]. 工業(yè)技術(shù)創(chuàng)新， 2017， 4（4）： 53-58.

[6] 張闊. 各向同性負泊松比超材料彈性參數(shù)線性化初探[J]. 工業(yè)技術(shù)創(chuàng)新， 2018， 5（4）： 53-60.

[7] 張闊. 關(guān)于半整數(shù)維波動方程求解的幾點注記[J]. 新一代信息技術(shù)， 2019， 2（21）： 13-21.

[8] 張闊. 增材制造技術(shù)中基于時域有限元方法的聲波—彈性波耦合（一）： 理論[J]. 工業(yè)技術(shù)創(chuàng)新， 2019， 6（5）： 74-85， 90.

[9] 張闊. 增材制造技術(shù)中基于時域有限元方法的聲波—彈性波耦合（二）： 應(yīng)用[J]. 工業(yè)技術(shù)創(chuàng)新， 2019， 6（6）： 75-86.

[10] 張闊. 柱狀模型傳遞矩陣校正圖譜[J]. 新一代信息技術(shù)， 2020， 3（2）： 1-9.

[11] 李信真， 車剛明， 歐陽潔， 等. 計算方法： 第2版[M]. 西安： 西北工業(yè)大學(xué)出版社， 2010.

[12] 張海瀾， 王秀明， 張碧星. 井孔的聲場和波[M]. 北京： 科學(xué)出版社， 2004.

[13] 王元明. 數(shù)學(xué)物理方程與特殊函數(shù)[M]. 北京： 高等教育出版社， 2004.