二維自由振動問題的自適應有限元分析初探

袁 駟，孫浩涵

(清華大學土木工程系，土木工程安全與耐久教育部重點實驗室，北京 100084)

自振頻率和振型反映了結構動力特性，是抗震分析和結構設計的重要基礎，根據(jù)已有結構的振型信息判斷其受損分布也是當前損傷識別研究的熱點[1]。二維自由振動在數(shù)學上歸結為偏微分方程特征值問題，精確求解缺乏可行性，因而集中質(zhì)量法、能量法和有限元法(FEM)等一系列離散和近似方法得到了充分發(fā)展。憑借處理復雜結構的靈活性和通用性，各類有限元法逐漸成為求解的主要手段[2]。

自適應有限元法(AFEM)是提升求解質(zhì)量和效率的一種有效方法[3―4]，主要思路是反復在前一次網(wǎng)格的求解信息基礎上，利用可靠的誤差估計手段對有限元解精度做出估計，并結合有效的網(wǎng)格細化技術調(diào)整網(wǎng)格，以獲得優(yōu)化的網(wǎng)格，直至滿足用戶事先給定的誤差限。

單元能量投影(Element Energy Projection，簡稱EEP)法是袁駟等[5―6]基于數(shù)學理論和力學概念，提出的新型有限元超收斂后處理算法，基于其超收斂的優(yōu)良特性，一整套有限元自適應分析方法已得到充分的發(fā)展和應用。該算法基于位移超收斂解估計誤差，指導網(wǎng)格細分，可給出按最大模度量的、逐點滿足給定誤差限的解答，在二維[7―9]及三維[10]線性問題中取得一系列成功，并已成功推廣至部分非線性問題[11―12]。

基于EEP法的一維FEM和二維FEMOL特征值問題自適應分析已經(jīng)取得成功[13―14]。本文方法基于類似的線性化思想，合理引入二維線性問題的EEP自適應分析技術，提出了二維自由振動問題的自適應有限元求解策略。文末若干算例表明，本法可靠，冗余度小，效果頗佳。

1 基本模型及有限元求解

本文以彈性薄膜為模型問題，進行相應的公式推導及算法說明。不失一般性，該算法流程可推廣至更加復雜的線彈性體自由振動及失穩(wěn)分析。

1.1 彈性薄膜自由振動問題

彈性薄膜自由振動在數(shù)學上對應于二維Laplace特征值問題，其控制方程可以表示為：

式中：L為線性微分算子；▽2為Laplace算子；λ=ω2ρ/T，ω為自振頻率，ρ為薄膜密度，T為薄膜張力；u為對應的振型函數(shù)。為方便討論，沿用數(shù)學表達，下文也將λ稱為特征值，u稱為特征函數(shù)，將二者合稱為特征對。此外，還應滿足適當給定的邊界條件。

1.2 二維有限元離散

本文采用雙m次四邊形單元，單元試探函數(shù)uh借助雙m次形函數(shù)對結點位移插值得到：

式中：Ni和Nj為m次Lagrange或Hierarchical形函數(shù)；dij為相應的結點位移。單元檢驗函數(shù)vh采用相同的形函數(shù)插值形式。

由于二維EEP超收斂計算要求，本方案的有限元網(wǎng)格須是“擬線法網(wǎng)格”，即先利用FEMOL的離散方式用一組結線對求解區(qū)域進行半離散，然后再沿結線維度進一步離散，得到的網(wǎng)格既是FEMOL常微分方程組(Ordinary Differential Equations，簡稱ODEs)的廣義一維有限元網(wǎng)格，也是原問題的二維有限元網(wǎng)格，如圖1所示。

圖1 二維問題逐維離散示意圖Fig.1 D-by-D discretization of 2D problems

給定網(wǎng)格劃分π，標準的FEM過程將建立如下的矩陣廣義特征值問題：

式中：D為振型向量；K和M分別為傳統(tǒng)FEM中的整體靜力剛度矩陣和一致質(zhì)量矩陣，二者均與λ無關。

1.3 矩陣廣義特征值問題求解

經(jīng)由1.2節(jié)的有限元過程，在當前網(wǎng)格π下，原問題1近似轉(zhuǎn)化為式(3)所示的矩陣廣義特征值問題。在本文算法中，采用基于Sturm序列的計數(shù)法及逆冪(子空間)迭代進行求解。

對于式(3)，設任給一試探值λa并作特征值移位，則可寫作如下形式：

式中，Ka=Kλ(λa)。根據(jù)Sturm序列的性質(zhì)，低于λa的特征值總數(shù)可通過J計數(shù)公式求得，即：

式中，s{Ka}代表對Ka符號的計數(shù)，即以Gauss消元法將Ka消成上三角陣后主對角線上負元素的個數(shù)。

以J計數(shù)為基礎，可通過二分法確定待求階特征值的區(qū)間范圍(λl,λu)，從而獲得較好的移位值λa∈(λl,λu)(一般可取為λa=(λl+λu)/2)，而后對問題3進行特征值移位后的逆冪(子空間)迭代(即對式(4)進行求解)，并不斷更新上、下界(λl,λu)，既保證了特征值的精度，也大幅提高了收斂速度。

2 誤差估計及自適應求解

2.1 自適應求解目標

設要求解前n階特征對(λk,uk) (k=1,2,…,n)，用戶給定特征值及特征函數(shù)的誤差限均為tol。本文算法的最終目標是：在精確解(λk,uk)(k=1,2,…,n)未知的情況下，得到充分好的FEM網(wǎng)格πk，使得該網(wǎng)格下的FEM解(k=1,2,…,n)同時滿足。

在實際計算時，由于精確解λk、uk未知，上述目標不能作為停機準則，因而采用以下誤差控制準則：

2.2 誤差估計及自適應網(wǎng)格劃分

二維有限元誤差估計是本文算法的核心，主要思想是通過等價線性問題[11]，將二維線性有限元EEP超收斂計算方法引入非線性的特征值問題求解中。與傳統(tǒng)的基于超收斂分片恢復SPR等方法的誤差估計[15―17]不同，在獲取超收斂的位移場后，我們采用最大模而非能量模進行誤差控制，從而直接控制特征函數(shù)求解的誤差，十分自然合理。為簡化表示，此處以(λh,uh)表示第k階特征對在當前FEM網(wǎng)格下的解。此時若將式(1)的右端項λu代換為λhuh，則其轉(zhuǎn)化為如下線性問題：

由于逆冪(子空間)迭代已經(jīng)收斂，若無網(wǎng)格調(diào)整，式(10)的FEM解仍為(λhuh)。又由于此時原非線性問題已轉(zhuǎn)化為式(10)的線性問題，從而可引入二維線性EEP超收斂求解公式進行后處理計算，獲得合理的誤差估計。

限于篇幅，這里不討論二維線性有限元EEP超收斂公式的具體推導過程，僅簡單介紹思路。與“逐維離散”方式相呼應，二維EEP法超收斂計算采用的是“逐維修復”方案。如圖2所示，借助有限元線法作為橋梁，分兩步利用單元能量投影定理導出單元上任一點的超收斂計算公式。其主要過程為：

圖2 二維問題逐維恢復示意圖Fig.2 D-by-D recovery of 2D problems

1) 二維有限元解：在當前網(wǎng)格πk下，經(jīng)由第1節(jié)所述的有限元過程獲得第k階的有限元解(λh,uh)。

2) 擬線法解：采用式(11)對線法常微分方程組的廣義一維有限元解進行超收斂計算得到“擬線法解”，用之代替線法解，其計算公式為：

3) 全域超收斂解：將問題視為有限元線法的超收斂計算，再次利用單元能量投影定理，通過式(12)獲得全域的二維有限元超收斂解，其計算公式為：

式(11)和式(12)中各參數(shù)的含義及計算請詳見文獻[7―8]，在此不再贅述。通過以上過程獲得EEP超收斂解后，即可定義單元上任意一點的估計誤差：

在本文中，基于該誤差估計，采用均差法[9]進行網(wǎng)格細分，其高效穩(wěn)健特性已被大量數(shù)值算例驗證。該法生成的網(wǎng)格分布合理、誤差冗余度小。

2.3 自適應求解步驟

與基于FEMOL的特征值問題自適應求解類似[13]，對每階待求的特征對，基于以上分析本文算法可歸結為：

1) 劃界階段：在當前二維FEM網(wǎng)格上(首階網(wǎng)格π0人為給定，第k+1階初始網(wǎng)格取為第k階最終網(wǎng)格)，基于J計數(shù)采用二分法得到待求特征值的上、下界。

2) 定解階段：采用移位的逆冪(子空間)迭代法得到該網(wǎng)格上特征對的有限元解，利用二維有限元的EEP超收斂公式計算超收斂解，檢驗誤差并進行必要的網(wǎng)格細分；在獲得的新網(wǎng)格上重復該迭代過程，不斷更新待求特征值的上、下界及相應特征函數(shù)，直至滿足誤差要求。

該算法高效、精確、通用、可靠，實現(xiàn)了精確的頻率和振型求解。相對于課題組之前提出的特征值問題FEMOL自適應分析算法，本文算法在二維的兩個方向進行離散，進一步降低了冗余度，提高了求解效率，增強了靈活性。

3 數(shù)值算例

本文算法已編制成Fortran90程序。若無特殊說明，下文算例均采用三次元進行計算，誤差限給定為tol=10-3。為保證超收斂計算的準確性，在進行逆冪(子空間)迭代時，收斂指標設置為更加苛刻的10-8，即當前后兩次迭代符合

認為該次逆冪(子空間)迭代已經(jīng)收斂。下文中，第k階最終特征值取為估計上、下界的中值，即誤差比較采取特征值混合誤差

和振型絕對最大誤差

并相應定義特征值誤差比ελ/tol及振型誤差比εu/tol，誤差比小于1即表明自適應求解圓滿完成。下面通過若干典型的彈性薄膜自由振動算例表明算法的有效性。

例1.固支方膜

考慮圖3所示的四邊固支方膜的自由振動問題，其精確解如下：

其中，m(m=1,2,3,…)及n(n=1,2,3,…)分別表示在x及y方向的半波數(shù)目，當二者交換取值時將出現(xiàn)重頻。采用本文方法求解前10階特征對，并將結果匯總至表1。

圖3 固支方膜計算簡圖Fig.3 Clamped square membrane model

表1 固支方膜自由振動求解結果Table 1 Results for clamped square membrane vibration

從表1結果可以看到，該算法能夠很好地進行誤差控制，各階特征對求解均滿足誤差要求，且從振型誤差比可以看到，本方法冗余度小。重頻情況可以有無窮多種振型組合，因此其振型誤差比不在此列出。

圖4為固支方膜的第1階振型圖，圖5和圖6分別為最終網(wǎng)格下的第1階振型真實誤差分布及估計誤差分布，可以看到，針對有限元解進行的EEP后處理過程在全域給出了相當準確的誤差估計，從而可以有效地指導網(wǎng)格劃分。本文方法對于重頻求解依然有效，以第5階及第6階的重頻為例，得到的振型分別如圖7、圖8所示，對應于m=1,n=3及m=3,n=1的情況。

圖4 固支方膜第1階FEM解振型圖Fig.4 1st order vibration mode of clamped square membrane by FEM

圖5 固支方膜第1階振型域內(nèi)誤差分布Fig.5 Error distribution of 1st order vibration mode of clamped square membrane

圖6 固支方膜第1階振型域內(nèi)估計誤差分布Fig.6 Error estimation of 1st order vibration mode of clamped square membrane

圖7 固支方膜第5階FEM解振型圖Fig.7 5th order vibration mode of clamped square membrane by FEM

圖8 固支方膜第6階FEM解振型圖Fig.8 6th order vibration mode of clamped square membrane by FEM

例2.固支L形膜

考慮L形固支薄膜自由振動問題。計算簡圖如圖9所示，初始網(wǎng)格采用2個四邊形單元。該例中求解域的凹角造成了一定的奇異性，從而增大了求解的難度。作為經(jīng)典的具有奇異性的問題，已有很多學者對該問題進行過研究，文獻[18]給出了其較為準確的頻率上、下界估計。

圖9 固支L形膜計算簡圖(左)及初始網(wǎng)格(右)Fig.9 Clamped L-shape membrane model (left) and initial mesh (right)

圖10為首階特征對自適應求解的最終網(wǎng)格。不同于文獻[13]中采用特殊的射線網(wǎng)格以適應問題的奇異性，我們看到隨著自適應過程的進行，網(wǎng)格自動地向奇異點附近加密以適應問題性質(zhì)。圖11和圖12給出了具有典型性的第1階及第5階振型圖。表2給出了該問題前10階求解結果，均嚴格滿足求解目標。

圖10 固支L形膜第1階特征對自適應求解最終網(wǎng)格Fig.10 Final mesh of adaptive procedure for 1st order eigenpair of clamped L-shaped membrane

圖11 固支L形膜第1階FEM解振型圖Fig.11 1st order vibration mode of clamped L-shaped membrane by FEM

圖12 固支L形膜第5階FEM解振型圖Fig.12 5th order vibration mode of clamped L-shaped membrane by FEM

表2 固支L形膜自由振動求解結果Table 2 Results for clamped L-shaped membrane vibration

例3.帶裂縫方膜

考慮具有奇異性的帶裂縫方膜問題，其計算模型及初始網(wǎng)格如圖13所示。根據(jù)對稱性，取半結構進行計算，只計算對稱情況的自由振動。該問題可定義為：

圖13 帶裂縫方膜Fig.13 Square membrane with an inside crack

在文獻[19]中，作者采用Trefftz方法對Laplace算子的特征值問題進行求解，并對該帶裂縫方膜問題給出了前兩階數(shù)值精確的特征值，可用來檢驗本文算法的可靠性。圖14為自適應求解第1階特征對的最終網(wǎng)格，可以看到無須人為干涉，自適應過程便自動識別了該問題在O點的強奇異性。

圖14 帶裂縫方膜第1階特征對自適應求解最終網(wǎng)格Fig.14 Final mesh of adaptive procedure for 1st order eigenpair of crack membrane

圖15和圖16分別為第1階和第3階振型圖。表3列舉了前5階的求解結果，與文獻中前2階的結果相比較，本算法獲得了令人滿意的效果。

圖15 帶裂縫方膜第1階FEM解振型圖Fig.15 1st order vibration mode of crack membrane by FEM

圖16 帶裂縫方膜第3階FEM解振型圖Fig.16 3rd order vibration mode of crack membrane by FEM

表3 帶裂縫方膜自由振動求解結果Table 3 Results for crack membrane vibration

例4.固支圓膜

考慮定義在半徑為1的圓形區(qū)域的Laplace特征值問題，計算簡圖及初始網(wǎng)格如圖17所示。其極坐標下解析的特征值及特征函數(shù)[20]為：

式中：Jm為第一類m階Bessel函數(shù)；jm,n為其第n階根，即：

圖17 固支圓膜計算簡圖(左)及初始網(wǎng)格(右)Fig.17 Clamped circular membrane model (left) and initial mesh (right)

對于該圓膜自由振動問題，文獻[21]也通過理論分析和數(shù)值方法給出了自振頻率。我們采用四次精確幾何單元進行計算，將前10階的計算結果匯總至表4。由于圓膜的對稱性，在有限元求解過程中我們發(fā)現(xiàn)對于其中若干階次，將求得重頻率的情況。但遵循文獻[20―21]的表示，表4中不再對其單獨列舉。由于該問題并不具有奇異性，從圖18所示的自適應最終網(wǎng)格我們也可以看到其性質(zhì)相對均勻，再次說明基于EEP誤差估計的網(wǎng)格劃分具有良好的合理性。圖19～圖22分別為有限元求得的固支圓膜第3階、4階、5階、7階振型。

表4 固支圓膜自由振動求解結果Table 4 Results for clamped circular membrane vibration

圖18 固支圓膜第1階、3階、5階自適應最終網(wǎng)格Fig.18 Final mesh of 1st, 3rd and 5th order with adaptive procedure for the clamped circular membrane

圖19 固支圓膜第3階FEM解振型圖Fig.19 3rd order vibration mode of clamped circular membrane by FEM

圖20 固支圓膜第4階FEM解振型圖Fig.20 4th order vibration mode of clamped circular membrane by FEM

圖21 固支圓膜第5階FEM解振型圖Fig.21 5th order vibration mode of clamped circular membrane by FEM

圖22 固支圓膜第7階FEM解振型圖Fig.22 7th order vibration mode of clamped circular membrane by FEM

圖23及圖24分別給出了第4階振型自適應最終網(wǎng)格下的實際誤差及估計誤差，全域上逐點誤差比均滿足要求，且估計誤差分布與實際誤差分布接近。

圖23 固支圓膜第4階振型域內(nèi)誤差分布Fig.23 Error distribution of 4th order vibration mode of clamped circular membrane

圖24 固支圓膜第4階振型域內(nèi)估計誤差分布Fig.24 Error estimation of 4th order vibration mode of clamped circular membrane

4 結論

本文提出的基于EEP法的二維自由振動問題有限元自適應分析方法，通過合理的線性化步驟，將業(yè)已成熟的二維線性有限元EEP后處理技術引入，思路清晰，結果可靠，可得到滿足精度要求的自振頻率和按最大模度量滿足用戶給定誤差限的振型，展示出穩(wěn)定的優(yōu)良特性。本文算法可拓展性強，未來結合本課題組在網(wǎng)格局部加密方面的最新研究成果[22]，有望進一步提高求解效率。