COMSOL仿真在聲波諧振管實(shí)驗(yàn)中的應(yīng)用

李春熠，胡安琦，張格格，白在橋，俞驍翀

(北京師范大學(xué) 物理學(xué)系，北京 100875)

隨著計(jì)算技術(shù)的進(jìn)步，數(shù)值仿真已經(jīng)成為與理論和實(shí)驗(yàn)并列的第三種科學(xué)研究方法.大學(xué)物理教學(xué)實(shí)驗(yàn)的選題一般比較簡(jiǎn)單，理論分析相對(duì)完善，但有時(shí)也會(huì)遇到涉及復(fù)雜幾何形狀或邊界條件的波動(dòng)或流體方面的問(wèn)題.此時(shí)簡(jiǎn)化的理論往往不能給出精確的結(jié)果或精細(xì)的物理圖像.如果針對(duì)這些實(shí)驗(yàn)把數(shù)值計(jì)算引入教學(xué)，不僅可以加深學(xué)生對(duì)所涉及物理問(wèn)題的理解，還可以為他們以后的學(xué)習(xí)和研究打下一定基礎(chǔ).本文以聲波諧振管為例，介紹數(shù)值仿真方法在基礎(chǔ)物理實(shí)驗(yàn)教學(xué)中的應(yīng)用.

聲波在管中傳播的衰減很小，在低頻時(shí)可以用一維聲波近似，因此實(shí)驗(yàn)教學(xué)中常利用聲波諧振管來(lái)研究波的反射、疊加和諧振等現(xiàn)象[1-3].但實(shí)際空氣中的聲波是由三維波動(dòng)方程:

(1)

描述的.該方程主要的約束是在硬邊界(比如管壁)處聲壓p的法向?qū)?shù)為零.雖然一維波模型可以給出諧振頻率的計(jì)算公式，但公式中的管口修正項(xiàng)并不能從一維波模型得到.此外，一些非理想情況也需要借助三維波動(dòng)方程理解.

有限元法(Finite Element Method，F(xiàn)EM)是一種非常有效的計(jì)算偏微分方程近似解的數(shù)值方法.有限元法起源于20世紀(jì)中期對(duì)彈性力學(xué)問(wèn)題的數(shù)值求解[4,5].后來(lái)擴(kuò)展到計(jì)算流體力學(xué)、熱學(xué)、電磁學(xué)、聲學(xué)等領(lǐng)域，在科學(xué)研究和工程中有十分廣泛的應(yīng)用[6].隨著計(jì)算機(jī)輔助設(shè)計(jì)(CAD)、復(fù)雜幾何形狀的自動(dòng)網(wǎng)格劃分、穩(wěn)健的求解器等方面的進(jìn)步，性能優(yōu)秀、使用方便的有限元計(jì)算軟件包應(yīng)運(yùn)而生，COMSOL Multiphysics就是其中的一個(gè)重要代表[7].

本文將介紹如何利用COMSOL分析聲波諧振管中的物理問(wèn)題.

1 COMSOL使用方法

有限元方法的原理是將待求問(wèn)題的幾何區(qū)域網(wǎng)格化，每個(gè)單元是簡(jiǎn)單的幾何形狀.然后通過(guò)變分法最小化能量泛函或加權(quán)殘差法離散化，將偏微分方程轉(zhuǎn)換為包含所有單元的代數(shù)方程組，最后用數(shù)值方法求解得到微分方程的近似解[7].

本文選擇COMSOL的“壓力聲學(xué)/頻域(acpr)”物理場(chǎng)，通過(guò)求解頻域的三維亥姆霍茲方程對(duì)聲波諧振管實(shí)驗(yàn)進(jìn)行數(shù)值仿真.下面以測(cè)量部分封閉管的諧振頻率為例說(shuō)明仿真計(jì)算的步驟.在這個(gè)實(shí)驗(yàn)中，設(shè)想閉管的端面上有一個(gè)孔，通過(guò)改變孔的面積，可以實(shí)現(xiàn)從閉管到開管的連續(xù)過(guò)渡.設(shè)計(jì)這個(gè)實(shí)驗(yàn)的初衷是想探究閉管密封不嚴(yán)對(duì)諧振頻率的影響.為了驗(yàn)證仿真結(jié)果的可靠性，本文還完成了真實(shí)實(shí)驗(yàn)，實(shí)驗(yàn)裝置參照?qǐng)D1.

諧振管

1.1 構(gòu)建三維幾何模型

本實(shí)驗(yàn)用到的幾何模型包括諧振管及兩端的附件、揚(yáng)聲器、有不同孔的有機(jī)玻璃擋片.除此之外，為了保證聲場(chǎng)方程的可解性，還需要定義“無(wú)窮遠(yuǎn)”處的邊界條件.構(gòu)建方法分述如下：

1) 諧振管管體

在三維直角坐標(biāo)系中，不妨設(shè)圓柱形諧振管的中軸線為z軸.以坐標(biāo)原點(diǎn)(0, 0, 0)為底面中心點(diǎn)，沿z軸正方向構(gòu)建高度900 mm、半徑20 mm的圓柱體(參見(jiàn)圖2).對(duì)圓柱體進(jìn)行“分層”，并設(shè)置“層的厚度”(即管壁厚度)為3 mm.

2) 揚(yáng)聲器

在諧振管左端管口附近[坐標(biāo)為(0, 0, -90)]構(gòu)建圓柱體，半徑與高度均為30 mm.在真實(shí)實(shí)驗(yàn)中，為防止揚(yáng)聲器反射的聲波進(jìn)入諧振管(從而在揚(yáng)聲器與管口之間發(fā)生共振)，揚(yáng)聲器不能正對(duì)管口[2].相應(yīng)地，將圓柱體繞x軸旋轉(zhuǎn)45°角(參見(jiàn)圖2).

圖2 諧振管管體和揚(yáng)聲器模型的三維圖像

3) 擋片

在右端管口和端蓋之間構(gòu)建帶有不同尺寸和形狀的開孔有機(jī)玻璃擋片(參見(jiàn)圖3).構(gòu)建方法為以(0, 0, 900)點(diǎn)為底面中心沿z軸正方向構(gòu)建圓柱體.沒(méi)有孔的擋片選擇不分層.圓形孔的擋片采用構(gòu)建圓柱體再“分層”的方法構(gòu)建，“層的厚度”為圓柱體的半徑減去圓形孔的半徑；正方形孔、長(zhǎng)方形孔的擋片則在不分層圓柱體的基礎(chǔ)上構(gòu)建正方體和長(zhǎng)方體；半圓形孔的擋片在不分層圓柱體的基礎(chǔ)上構(gòu)建圓柱體和長(zhǎng)方體取“交集”，如圖4所示.模型的具體尺寸參數(shù)為：“有機(jī)玻璃片”圓柱體高度2 mm，半徑23.5 mm；“圓形孔”圓柱體高度2 mm，半徑23.5 mm，層的厚度為7.5～21.3 mm，間隔2 mm；“正方形孔”長(zhǎng)方體高度2 mm，寬度、深度相等，為3～21 mm，間隔3 mm；“長(zhǎng)方形孔”長(zhǎng)方體高度2 mm，寬度24 mm，深度為2～24 mm，間隔2 mm；“半圓形孔”圓柱體高度2 mm，半徑為2～16 mm，間隔2 mm .

圖3 右端管口的yOz平面框架圖

閉管

4) 完美匹配層

為了使問(wèn)題有確定的解，COMSOL仿真用充分大的“完美匹配層”(Perfectly Matched Layers, PMLs)包圍求解區(qū)域. PMLs完全吸收聲波，不會(huì)產(chǎn)生反射，作用相當(dāng)于定義了“無(wú)窮遠(yuǎn)”處的邊界條件.構(gòu)建高度2 000 mm、半徑500 mm、層厚度300 mm的圓柱體將諧振管、揚(yáng)聲器等幾何體包含在內(nèi)，并將該圓柱體定義為“完美匹配層” (參見(jiàn)圖5).

圖5 完美匹配層模型的三維圖像

目前所建的幾何模型關(guān)于yOz平面對(duì)稱，建立一個(gè)大的長(zhǎng)方體與前面的幾何模型取“交集”，將中心截面(yOz平面)暴露出來(lái)，方便后續(xù)選擇區(qū)域和邊界，以及呈現(xiàn)內(nèi)部的聲壓分布等情況(參見(jiàn)圖6).

圖6 對(duì)半“切割”后的幾何模型的三維圖像

1.2 定義物理場(chǎng)和驅(qū)動(dòng)

選擇管內(nèi)區(qū)域、擋片開孔區(qū)域以及完美匹配層區(qū)域?yàn)榍蠼獾目諝饨橘|(zhì)區(qū)域(參見(jiàn)圖7),并將以上區(qū)域做“對(duì)稱”操作.空氣介質(zhì)區(qū)域的邊界(包括管體、揚(yáng)聲器、擋片等固體)會(huì)自動(dòng)被定義為“硬聲場(chǎng)邊界”.硬聲場(chǎng)邊界可視為理想的剛體.

左端管口局部圖 右端管口局部圖

將材料庫(kù)中的空氣“Air”添加到物理場(chǎng)中，采用“線彈性”流體模型，不考慮空氣的熱傳導(dǎo)、黏性等耗散因素. 空氣溫度設(shè)置為真實(shí)實(shí)驗(yàn)溫度值296.95 K(23.8 ℃)，絕對(duì)壓強(qiáng)設(shè)為1 atm.

定義驅(qū)動(dòng)的方法是：“添加聲激勵(lì)”，給揚(yáng)聲器圓柱體面朝向左端管口的底面(參見(jiàn)圖8)施加“法向位移”1 mm，使之做簡(jiǎn)諧運(yùn)動(dòng)，模擬揚(yáng)聲器的膜片振動(dòng).

圖8 揚(yáng)聲器的振動(dòng)面

最后根據(jù)物理場(chǎng)劃分網(wǎng)格.網(wǎng)格劃分需要設(shè)置“最小單元大小”“最大單元增長(zhǎng)率”等多個(gè)參數(shù)或選項(xiàng).對(duì)于本文的問(wèn)題，從物理上考慮，管內(nèi)單元的尺寸應(yīng)該小于波長(zhǎng)的1/6(參見(jiàn)圖9).

左端管口局部圖

1.3 求解

對(duì)于聲波諧振管，最基本的測(cè)量對(duì)象是諧振頻率.添加“特征頻率”的“研究”，設(shè)置特征頻率搜索基準(zhǔn)值和所需特征頻率數(shù)，點(diǎn)擊“計(jì)算”，讓求解器在基準(zhǔn)值附近按照“最近的絕對(duì)值”的方式搜索本征頻率.以一端開口、一端封閉的情形(閉管)為例，設(shè)800 Hz為基準(zhǔn)值，計(jì)算結(jié)果見(jiàn)表1.需要注意的是，雖然模型忽略了空氣的耗散，但完美匹配層的吸聲作用使得系統(tǒng)仍存在耗散，所以諧振頻率一般為復(fù)數(shù)，其虛部可理解為共振峰的寬度或耗散的大小(快慢).仿真結(jié)果給出了很多諧振頻率，其中只有虛部很小的才有物理意義(表1中方括號(hào)標(biāo)記).圖10顯示了諧振頻率662.95+0.6137 9i對(duì)應(yīng)諧振管內(nèi)的聲壓分布.不同灰度代表不同的聲壓值，白色為聲壓最大值，黑色為聲壓最小值，黑白交替處的聲壓為零(波節(jié)).圖中可見(jiàn)明顯的3個(gè)波節(jié)，因此模式數(shù)n=4.

圖10 閉管發(fā)生諧振時(shí)的聲壓分布圖

實(shí)驗(yàn)中經(jīng)常會(huì)研究諧振頻率對(duì)某個(gè)變量的依賴關(guān)系.這時(shí)可以將該變量設(shè)置為參數(shù).添加“參數(shù)化掃描”，設(shè)置起始值、終止值和變化步長(zhǎng)，使求解器自動(dòng)在參數(shù)的定義域內(nèi)掃描求解.例如，在改變長(zhǎng)方形孔的深度時(shí)，可將其深度設(shè)置為參數(shù)d，為了使其關(guān)于yOz平面對(duì)稱，y的位置設(shè)為-d/2，使d從2 mm到22 mm變化，步長(zhǎng)為2 mm.這樣就可以得到不同寬度矩形孔對(duì)應(yīng)的諧振頻率.

表1 閉管時(shí)的諧振頻率仿真結(jié)果

續(xù)表

還可以添加“頻域”的“研究”，設(shè)置起始頻率、終止頻率和變化步長(zhǎng)，點(diǎn)擊“計(jì)算”，即可自動(dòng)繪制出各頻率(包括諧振頻率和非諧振頻率)下的聲壓分布圖和聲壓級(jí)分布圖. 可將任何繪圖的數(shù)據(jù)以文本格式或者CSV格式導(dǎo)出，便于進(jìn)一步分析.

下一節(jié)將結(jié)合具體實(shí)驗(yàn)對(duì)仿真結(jié)果進(jìn)行更仔細(xì)的分析.

2 仿真實(shí)驗(yàn)及結(jié)果分析

這一節(jié)介紹幾個(gè)用COMSOL仿真的實(shí)驗(yàn).

2.1 管口附近的聲壓分布

在一維波近似下，兩端都為開口的諧振管(稱為開管)的諧振頻率經(jīng)驗(yàn)公式為[2]

(2)

式中c為聲速，n=1,2,3, …為模式數(shù)，而等效管長(zhǎng)為

L=L0+0.8d

(3)

其中L0為實(shí)際管長(zhǎng)，d為管的內(nèi)直徑.0.8d稱為管口修正(end correction).從表達(dá)式上看，管口修正可理解為：聲波可以從管口延伸出來(lái)，在距離管口0.4d的位置形成聲壓的波節(jié).

為了探究管口修正的真實(shí)物理意義，本文計(jì)算了開管情況下的聲壓分布.

以開管的模式數(shù)n=5，特征頻率為936.47 Hz的情況為例.定義點(diǎn)(0, 0, -500)到點(diǎn)(0, 0, 1 500)的“三維截線”(即中心軸線).以該截線的z軸坐標(biāo)為橫坐標(biāo)，聲壓為縱坐標(biāo)繪制一維聲壓曲線，結(jié)果見(jiàn)圖11.為便于比較，圖中還畫出來(lái)一維聲波模型預(yù)言的聲壓曲線(正弦函數(shù)曲線).由圖像可看出，在諧振管內(nèi)部(0～900 mm)，仿真實(shí)驗(yàn)曲線與理論曲線基本重合.但在管口外，仿真曲線迅速衰減，與一維聲波模型偏離較大.

圖11 中心軸線上聲壓的仿真結(jié)果(實(shí)線)與一維模型所得的理論曲線(虛線)

為了進(jìn)一步探究中心軸線聲壓的實(shí)驗(yàn)曲線與簡(jiǎn)化理論曲線的偏離細(xì)節(jié)，在右管口(z=900 mm)附近選取15個(gè)截面(xOy平面)，繪制截面上的聲壓分布圖.截面的z坐標(biāo)范圍從893 mm 到907 mm，間隔為1 mm.為了便于通過(guò)顏色進(jìn)行定性比較，將灰度圖顏色范圍固定為-0.1～1.5 Pa.15個(gè)截面上的聲壓分布如圖12所示. 可以看出，在管內(nèi)部遠(yuǎn)離開口處，聲壓均勻分布，符合一維波模型.在開口附近，不論是管內(nèi)還是管外，聲壓分布不再均勻，呈現(xiàn)聲壓隨著半徑增大而減小的分布規(guī)律.

還可以繪制管口附近的聲壓等值面(圖13).同樣，在管內(nèi)距離管口較遠(yuǎn)的位置，聲壓等值面基本垂直于中心軸，這意味著一維波模型比較符合實(shí)際.但越靠近管口，等壓面越向外凸出為弧面(管口類似一個(gè)聲源)，此時(shí)一維波模型不再符合實(shí)際情況.

圖12 開管時(shí)管口附近不同位置的徑向聲壓分布圖：共15個(gè)截面，z坐標(biāo)為893～907 mm，間隔1 mm

圖13 開管情況管口附近的聲壓等值面圖

2.2 麥克風(fēng)對(duì)諧振頻率的影響

在聲波諧振管實(shí)驗(yàn)中，使用微型麥克風(fēng)在管內(nèi)移動(dòng)以測(cè)量各處的聲壓.麥克風(fēng)的存在必然會(huì)影響管內(nèi)的聲場(chǎng)分布，從而改變諧振頻率.這里定量評(píng)估麥克風(fēng)對(duì)諧振頻率的影響，考慮的影響因素包括麥克風(fēng)的半徑和在管內(nèi)的位置.為簡(jiǎn)單起見(jiàn)，只考慮開管情況.

構(gòu)建圓柱體模擬麥克風(fēng)，設(shè)置半徑為3.5 mm，高度為10 mm(真實(shí)實(shí)驗(yàn)所用麥克風(fēng)的尺寸).為使麥克風(fēng)緊貼管內(nèi)壁，圓柱體底面中心的坐標(biāo)x=0 mm，y=-13.5 mm.z坐標(biāo)設(shè)為參數(shù)，從0 mm到890 mm改變，步長(zhǎng)為10 mm.圖14畫出了在振動(dòng)模式數(shù)n=4、5時(shí)，諧振頻率的移動(dòng)量與麥克風(fēng)的位置的關(guān)系(管內(nèi)未放置麥克風(fēng)時(shí)的諧振頻率分別為

749.1 Hz和936.47 Hz).仿真結(jié)果顯示，麥克風(fēng)的確會(huì)影響諧振頻率，但總體上諧振頻率的移動(dòng)量都不大(都不超過(guò)1 Hz). 有趣的是，諧振頻率的移動(dòng)也呈正弦周期性變化的規(guī)律：在聲壓波節(jié)附近，諧振頻率向下移動(dòng)；而在聲壓波腹附近，諧振頻率向上移動(dòng).此外，麥克風(fēng)在波節(jié)處(相比于波腹處)對(duì)諧振頻率的影響要更大.

n=4

由于麥克風(fēng)在波腹和波節(jié)處對(duì)諧振頻率的影響最大，進(jìn)一步探究當(dāng)麥克風(fēng)在波腹或波節(jié)處時(shí)，其半徑對(duì)諧振頻率的影響.把麥克風(fēng)放置在管中央[底面中心坐標(biāo)x=0 mm，z=445 mm，y=(-17+r)mm]，其中半徑r設(shè)為參數(shù)，從1 mm到11 mm改變，步長(zhǎng)為2 mm.由于開管情況下聲壓的分布具有對(duì)稱性，當(dāng)模式數(shù)n為偶數(shù)時(shí)，麥克風(fēng)處于波節(jié)處；模式數(shù)n為奇數(shù)時(shí)，麥克風(fēng)處于波腹處.計(jì)算不同半徑的麥克風(fēng)引起的諧振頻率移動(dòng)，結(jié)果見(jiàn)表2和圖15.

表2 存在不同半徑麥克風(fēng)時(shí)的諧振頻率仿真結(jié)果

圖15 不同半徑麥克風(fēng)引起的諧振頻率移動(dòng)

仿真結(jié)果顯示麥克風(fēng)的半徑r越大，諧振頻率的移動(dòng)量越大.當(dāng)r較小時(shí)，頻率移動(dòng)量正比于r2(由于圓柱體高度固定，也正比于麥克風(fēng)的體積).此外，同在波節(jié)(或波腹)處，模式數(shù)n越大，諧振頻率移動(dòng)越多.最后一點(diǎn)容易理解：模式數(shù)n越大，波長(zhǎng)越短，麥克風(fēng)的大小就越不可以忽略.以上現(xiàn)象在此前的教學(xué)中并沒(méi)有注意.雖然實(shí)驗(yàn)中麥克風(fēng)的半徑較小，總體上對(duì)諧振頻率的影響不大，但這仍是一個(gè)值得注意和進(jìn)一步研究的問(wèn)題.

2.3 一端部分封閉管的諧振頻率

這個(gè)實(shí)驗(yàn)的意義前文已有說(shuō)明，下面給出實(shí)驗(yàn)結(jié)果.表3列出了模式數(shù)n=5下的圓形孔、正方形孔、長(zhǎng)方形孔和半圓形孔擋片對(duì)應(yīng)諧振頻率的實(shí)測(cè)值(fa)、仿真值(fs).可以看出，絕大部分仿真值與實(shí)測(cè)值的差在1 Hz以內(nèi)，而實(shí)驗(yàn)測(cè)量的重復(fù)性誤差也在1 Hz左右.這表明，對(duì)確定的物理模型，可以用虛擬仿真實(shí)驗(yàn)代替實(shí)際實(shí)驗(yàn)，得到準(zhǔn)確的結(jié)果.

表3 部分封閉管的諧振頻率實(shí)驗(yàn)結(jié)果與仿真結(jié)果

3 總結(jié)與討論

本文介紹了COMSOL數(shù)值仿真在聲波諧振管實(shí)驗(yàn)中一些應(yīng)用，包括：對(duì)管口處聲壓分布的仿真揭示了管口聲波的衰減形式，并給出一維波模型成立的大致范圍；仿真麥克風(fēng)對(duì)諧振頻率的影響，發(fā)現(xiàn)諧振頻率移動(dòng)的方向和大小與麥克風(fēng)的位置和振動(dòng)模式數(shù)有關(guān)；部分封閉管的仿真結(jié)果與實(shí)測(cè)結(jié)果符合很好，驗(yàn)證了仿真結(jié)果的可靠性.以上結(jié)果表明：數(shù)值仿真作為真實(shí)實(shí)驗(yàn)的補(bǔ)充，可以揭示更多的不便實(shí)際測(cè)量的物理細(xì)節(jié)，為建立正確的物理圖像提供了重要參考.同時(shí)，在一些情況下仿真實(shí)驗(yàn)可以設(shè)計(jì)成獨(dú)立的實(shí)驗(yàn)，也可以代替真實(shí)的物理實(shí)驗(yàn)，得出準(zhǔn)確的實(shí)驗(yàn)結(jié)果，從而大大降低實(shí)驗(yàn)成本并縮短實(shí)驗(yàn)時(shí)間.對(duì)于聲波諧振管實(shí)驗(yàn)，已有利用有限元分析設(shè)計(jì)仿真實(shí)驗(yàn)的教學(xué)實(shí)踐[8]. 此外，落球法測(cè)量液體黏性系數(shù)、聲速測(cè)量[9]、動(dòng)態(tài)法測(cè)量楊氏模量[10,11]、音叉共振[12]、弦振動(dòng)[13]等基礎(chǔ)物理實(shí)驗(yàn)也提供了可以進(jìn)行數(shù)值仿真深入研究的問(wèn)題.把數(shù)值仿真融入基礎(chǔ)物理實(shí)驗(yàn)教學(xué)，可以讓學(xué)生更深入地理解所涉及的物理問(wèn)題，同時(shí)提高學(xué)生解決實(shí)際物理問(wèn)題的能力.