童韞哲，范 軍，王 斌

(上海交通大學(xué)高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海200240)

0 引 言

頻散曲線在無(wú)損檢測(cè)、彈性目標(biāo)聲散射機(jī)理及水聲傳播領(lǐng)域都有著廣泛用途[1-5]。在水下目標(biāo)聲散射的研究當(dāng)中，大量的實(shí)驗(yàn)結(jié)果表明，水下目標(biāo)的回波中除了容易解釋的幾何散射波以外，總是會(huì)出現(xiàn)一些不易解釋的非鏡面散射波，這些散射波實(shí)際上是在彈性體中傳播的波再輻射產(chǎn)生的[6]。所以計(jì)算彈性波的頻散曲線，了解彈性波的傳播規(guī)律，對(duì)水下目標(biāo)聲散射的精細(xì)特征研究具有至關(guān)重要的作用。

對(duì)于頻散曲線的計(jì)算一般都?xì)w結(jié)為特征方程零點(diǎn)的求解：

其中：k 為行波波數(shù)；ω 為角頻率。目前求解公式(1)的方法主要有三種：(1) 譜方法。該方法是一種較傳統(tǒng)的解析方法。譜方法數(shù)值插值的方法直接求解特征方程，具有計(jì)算量小，計(jì)算精度高的優(yōu)點(diǎn)[7-11]。(2) 半解析半有限元(Semi-Analytical Finite Element,SAFE)法。該方法結(jié)合了有限元與解析解表達(dá)式的特點(diǎn)，只對(duì)波導(dǎo)截面進(jìn)行建模，因此降低了計(jì)算模型的維度[12]。(3) 波有限元(Wave Finite Element,WFE)方法。該方法利用波導(dǎo)的周期性，通過(guò)建立一部分波導(dǎo)的有限元模型以及周期邊界條件獲得相應(yīng)的剛度矩陣，并求解特征根[13]。

鑒于目前的商用有限元軟件的求解器能夠非常成熟地應(yīng)用周期邊界條件求解特征方程，本文引入了波有限元(WFE)方法求解頻散曲線[13]。常規(guī)有限元方法無(wú)法實(shí)現(xiàn)無(wú)限長(zhǎng)的波導(dǎo)建模，所以需要引入基于Floquet-Bloch 理論的周期性邊界條件。Floquet 理論是微分方程理論的一個(gè)分支，它是求解形如x (t ) = A (t )x (t )的周期性線性微分方程的一種方法，其中A (t )是周期為T 的分段連續(xù)周期函數(shù)，在固體物理學(xué)中叫做Bloch 理論。

1 Floquet-Bloch 理論在求解彈性體頻散曲線中的應(yīng)用

有限元方法求解彈性體頻散曲線的關(guān)鍵在于利用Floquet-Bloch 理論將長(zhǎng)度為L(zhǎng)、高度為h 的有限長(zhǎng)彈性體控制方程中的位移/振速寫成行波形式。省略時(shí)間因子eiωt，彈性體控制方程可以寫成：

其中：σ 是應(yīng)力張量，u 是位移張量。將位移u 的解寫成平面波的形式：

公式(3)表示某一特定振動(dòng)狀態(tài)f (y )沿著x 方向以行波的形式傳播，振動(dòng)狀態(tài)只是y 的函數(shù)，在傳播過(guò)程中不發(fā)生變化，長(zhǎng)度為L(zhǎng)、高度為h 的計(jì)算單元示意圖如圖1 所示。

圖1 計(jì)算單元示意圖Fig.1 Schematic diagram of calculation cell

長(zhǎng)度為L(zhǎng) 的單元左邊位移和右邊位移存在如下關(guān)系：

由于u 以行波方式沿著x 方向傳播，振動(dòng)狀態(tài)與x 無(wú)關(guān)。二維無(wú)限長(zhǎng)板可以近似認(rèn)為沿x 方向周期為L(zhǎng) 的周期性結(jié)構(gòu)，利用Floquet-Bloch 理論，存在如下關(guān)系：

其中：kFB是Floquet-Bloch 理論的指數(shù)因子，kFB與x 方向的波數(shù)kx有如下關(guān)系：

詳細(xì)理論推導(dǎo)參見(jiàn)文獻(xiàn)[13]。

2 有限元計(jì)算結(jié)果與譜方法的對(duì)比

有限元在求解頻散曲線時(shí)的困難在于無(wú)法建立無(wú)限長(zhǎng)的波導(dǎo)，F(xiàn)loquet-Bloch 理論實(shí)際上是對(duì)如圖1 所示的計(jì)算單元施加邊界條件以模擬無(wú)限長(zhǎng)波導(dǎo)。COMSOL 軟件中Floquet 周期邊界條件模塊能夠?qū)崿F(xiàn)這一功能，其公式為

其中：usrc和rsrc是源邊界的位移和坐標(biāo)；udst和rdst是目標(biāo)邊界的位移和坐標(biāo)，詳細(xì)信息參見(jiàn)文獻(xiàn)[14]。

在COMSOL 軟件中建立如圖2 所示的模型，在模型兩邊添加Floquet 周期性邊界條件，并輸入不同的行波波數(shù)kF，將有限長(zhǎng)區(qū)域模擬成無(wú)線長(zhǎng)波導(dǎo)，再利用特征頻率求解器求解相應(yīng)kF所對(duì)應(yīng)的特征頻率ω。

圖2 有限元計(jì)算示意圖Fig.2 Schematic diagram of finite element calculation

因?yàn)橹豢紤]沿x軸傳播的行波，所以向量kF在y 軸投影為0，令其在x 軸的投影為kx。由式(6)可知，當(dāng)n=1時(shí)，kx的取值范圍是0～2π /L，實(shí)際上kx取0～π /L最合理，具體情況會(huì)在下文討論。使用COMSOL 的搜根求解器可以求得kx滿足計(jì)算單元的特征方程所對(duì)應(yīng)的所有特征頻率f，由特定的kx及與其對(duì) 應(yīng)的特 征 頻 率f 可求得 對(duì)應(yīng) 的 相速度Cp= 2π f/ kx。本文計(jì)算了兩邊真空的彈性平板頻散曲 線，參數(shù)為：ρ=7900 kg·m-3，縱波聲速CL=5940m·s-1，剪切波速CT=3 100 m·s-1，板的厚度h=40mm。

圖3 有限元與譜方法計(jì)算結(jié)果對(duì)比圖Fig.3 Comparison of calculation results of WFE method and spectral method

圖3 為有限元與譜方法[3]計(jì)算結(jié)果對(duì)比，橫坐標(biāo)為頻率與板厚之積，縱坐標(biāo)是相速度Cp與剪切波速之比。兩者吻合得非常好，從而證明了基于Floquet-Bloch 理論的有限元建模方法的正確性。

3 對(duì)于計(jì)算單元長(zhǎng)度與高度之比值R的討論

選取合適的R(R 為計(jì)算單元的長(zhǎng)度L 與高度h之比)對(duì)于計(jì)算結(jié)果至關(guān)重要。本文所提到的方法是用Floquet-Bloch 周期邊界條件模擬無(wú)限長(zhǎng)的情況，所以存在一定的局限性，這一局限性主要體現(xiàn)在R 的大小與計(jì)算結(jié)果的關(guān)系。

圖4 是相同計(jì)算參數(shù)、不同R 值的計(jì)算結(jié)果，紅線分割了計(jì)算結(jié)果的有效區(qū)和無(wú)效區(qū)。由式(4)中 eikxL的周期性和對(duì)稱性可知波數(shù)kx必須滿足：

因此，kx所對(duì)應(yīng)的特征頻率f 和相速度Cp必然滿足：

對(duì)公式(7)乘以厚度h，經(jīng)過(guò)變換可以得到相速度pC 與頻率必須滿足：

圖4 不同R 計(jì)算結(jié)果對(duì)比Fig.4 Comparison of calculation results of WFE method for different R values

由圖4(a)～4(d)可以看出：計(jì)算結(jié)果的有效區(qū)域隨R 的減小而增大，這一規(guī)律完全符合公式(13)。此外，根據(jù)經(jīng)驗(yàn)可知，當(dāng)R=0.1 時(shí)計(jì)算結(jié)果較好，本文中的圖3 即是使用R=0.1 計(jì)算得到的結(jié)果。

4 結(jié) 論

本文簡(jiǎn)要闡述了Floquet-Bloch 理論及其在頻散曲線計(jì)算上的應(yīng)用。在此基礎(chǔ)之上，分析了其計(jì)算結(jié)果的局限性，為將此方法進(jìn)一步應(yīng)用于更加復(fù)雜的計(jì)算奠定了基礎(chǔ)。

本文得到以下主要結(jié)論：

(1) 利用基于Floquet-Bloch 理論的周期性邊界條件的有限元建模方法，能夠準(zhǔn)確計(jì)算彈性體的頻散曲線；

(2) 通過(guò)周期性邊界條件將有限長(zhǎng)彈性體近似成無(wú)限長(zhǎng)波導(dǎo)存在一定的局限性，這一局限性體現(xiàn)在計(jì)算單元長(zhǎng)度與高度之比R 對(duì)于計(jì)算結(jié)果正確區(qū)域的影響；

(3) 使用有限元法計(jì)算頻散曲線較之傳統(tǒng)的迭代法與譜方法有著操作簡(jiǎn)單、不受波導(dǎo)形狀影響等優(yōu)勢(shì)。